In all the examples of Unit 10 we based the estimation procedure on only one presence/absence vector.
This has played a fundamental role into the construction of the R code for maximum likelihood estimation.
Indeed, the likelihood function is defined as
LL<- function(x,y) {<br>
R = dbinom(o,1,inc(1:length(o),x,y),log = TRUE)<br>
-sum(R)<br>
}
where “o” is a presence/absence vector, like “p” or “p1” in our reference metapopulation REFMETAPOP.
The function LL uses the function dbinom, requiring a vector of integers as first argument (see Remark 10.1). Then, this function cannot be directly used as it is if the observed mean occupancy time is used, as in Eq. (5.1) and Example 5.2.
However, only a slight modification is required: we only need to consider a longer presence/absence vector, the vector obtained by collecting all the vectors together. Indeed, as mentioned in Unit 5, the temporal order of the data is not taken into account, and the likelihood becomes
| (11.1) |
|
\(J_1^{\delta_{11} + \delta_{21}}(1-J_1)^{1-(\delta_{11} + \delta_{21})}\) x ... x \(J_n^{\delta_{1n} + \delta_{2n}}(1-J_n)^{1-(\delta_{1n} + \delta_{2n})}\) |
in the case of two observations per patch (and analogously for more than 2). This requires some simple code gimmick.
Example 11.1. Basic IFM, with more than one snapshot.
We refer to REFMETAPOP, as usual.
pm<-(p+p1)/2
S<-rowSums(edis[,pm>0])
p2<-c(p,p1)
length(p2)
# for each presence/absence vector inc should be called, in order to
# match each patch observation with the right area and connectivity.
# Be aware to follow the same order in collecting the vectors as in p2
LL<- function(x,y) {<br>
R = dbinom(p2,1,c(inc(1:length(p),x,y),inc(1:length(p1),x,y)),log = TRUE)<br>
-sum(R)<br>
}
library(stats4)
fit <- mle(LL, start = list(x=0.03, y=0.01), method = "L-BFGS-B", lower = rep(0, 2))
summary(fit)
Example 11.2. The boundary problem.
As explicitly mentioned in Example 10.3, the optimization problem we are considering is a constrained problem,
as parameters x and y should be positive.
In this kind of problems, it can occur that the optimum belongs to the boundary that is, either \(\hat{x}\) or \(\hat{y}\).
The R instruction mle takes constraint into account by the option lower = rep(0, 2)
(equivalent to lower = c(0,0)).
However, the behavior of the log-likelihood function on the boundary domain can be explicitly checked, by fixing x = 0
and searching for the minimizing value of y,
fit1 <- mle(LL, start = list(x=0.03, y=0.01), fixed=list(x=0),
method = "L-BFGS-B", lower = rep(0, 2))
summary(fit1)
and viceversa:
fit2 <- mle(LL, start = list(x=0.03, y=0.01), fixed=list(y=0),
method = "L-BFGS-B", lower = rep(0, 2))
summary(fit2)
It can be easily seen that both the minus 2 x log-likelihood values are higher than the optimum value of Example 11.1, and therefore the minimum is reached out of the boundaries, as stated in Example 11.1.
Example 11.3. The basic IFM with rescue effect: more than one snapshot.
pm<-(p+p1)/2
S<-rowSums(edis[,pm>0])
# As in Example 10.4:
incRE<-function(i, x, y){<br>
(A[i]^x)*S[i]^2/((y^2*A0^x)+(S[i]^2*A[i]^x))<br>
}
# As in Example 11.1:
p2<-c(p,p1)
LL<- function(x,y) {<br>
R = dbinom(p2,1,c(incRE(1:length(p),x,y),incRE(1:length(p1),x,y)),log = TRUE)<br>
-sum(R)<br>
}
fit1 <- mle(LL, start = list(x=0.3, y=0.005), method = "L-BFGS-B", lower = rep(0, 2))
summary(fit1)