Nothing
svocc.fit <-
function(Y, X, Z, link.sta = "cloglog", link.det = "logit", penalized = FALSE, auc = FALSE,
method = c("optim", "dc"), inits, ...)
{
obs <- Y
occ <- X
det <- Z
## internal funs
`singleocc.MLE` <-
# function(parameters, observations, X, Z, link1, link2)
function(parameters, observations, X, Z, link1, link2)
{
par.occ <- parameters[1:NCOL(X)]
par.det <- parameters[-(1:NCOL(X))]
# linkinvfun.occ <- binomial(link1)$linkinv
# linkinvfun.det <- binomial(link2)$linkinv
phi <- linkinvfun.occ(X %*% par.occ)
detec <- linkinvfun.det(Z %*% par.det)
q <- detec * phi
loglik.bin <- sum(observations * log(q)) + sum((1 - observations) * log(1 - q))
if (is.na(loglik.bin))
loglik.bin <- -1e+05
return(loglik.bin)
}
#-------------------------------------------------------------------------------------------
#Name: LikelihoodPenalizedSingleSurvey.fn
#Function: To calculate the penalized log likelihood for a single survey occupancy model
#Input parameters:
# parameters: a vector with the value of the parameters for the occupancy and detection model.
# observations:vector of size n.sites that contains 1 if the specied was observed,0 otherwise
# X:matrix that contains the value of the covariates used in the occupancy model
# Z:matrix that contains the value of the covariates used in the detection model
# bmax:naive estimated parameters for the occupancy model
# dmax:naive estimated parameters for the detection model
# psi0:mean estimated occupancy obtained from the MLE
# d0:mean estimated detection obtained from the MLE
# link1:Link function for the occupancy model
# link2:Link function for the detection model
# v1:Square root of the sum of the estimated variance for the detection parameters obtained by the MLE
# v2:Square root of the sum of the estimated variance for the occupancy parameters obtained by the MLE
#Output: Value of the penalized log likelihood function evaluated at the value of the parameters
#-------------------------------------------------------------------------------------------
`singleocc.PMLE` <-
function(parameters, observations, X, Z, bmax, dmax, psi0, d0, link1, link2, v1, v2)
{
param.occup <- parameters[1:NCOL(X)]
param.detec <- parameters[-(1:NCOL(X))]
linkinvfun1 <- binomial(link1)$linkinv
linkinvfun2 <- binomial(link2)$linkinv
# occup <- linkinvfun1(X %*% param.occup)
# occup.max <- linkinvfun1(X %*% bmax)
# detec <- linkinvfun2(Z %*% param.detec)
# detec.max <- linkinvfun2(Z %*% dmax)
occup <- linkinvfun.occ(X %*% param.occup)
occup.max <- linkinvfun.occ(X %*% bmax)
detec <- linkinvfun.det(Z %*% param.detec)
detec.max <- linkinvfun.det(Z %*% dmax)
p <- occup * detec
LogLikZIB <- (observations * log(p)) + ((1 - observations) * log(1 - p))
penalty.occupancy <- (1 - mean(occup.max)) * d0 * v1 * sum(abs(param.occup - bmax))
penalty.detection <- (1 - mean(detec.max)) * psi0 * v2 * sum(abs(param.detec - dmax))
penalization <- penalty.occupancy + penalty.detection
LogLik <- sum(LogLikZIB) - penalization
return(LogLik)
}
## BUGS model for 2^3 kinds of link pairs
mcmcSS.all <- c("model {",
"for (i in 1:N.sites) {",
"Y[i] ~ dbin(p[i], k)",
"logit(p.occ[i]) <- inprod(X[i, ], beta)",
"logit(p.det[i]) <- inprod(Z[i, ], theta)",
"probit(p.occ[i]) <- inprod(X[i, ], beta)",
"probit(p.det[i]) <- inprod(Z[i, ], theta)",
"cloglog(p.occ[i]) <- inprod(X[i, ], beta)",
"cloglog(p.det[i]) <- inprod(Z[i, ], theta)",
"p[i] <- p.occ[i] * p.det[i]",
"}",
"for (j in 1:num.cov.occ) {",
"beta[j] ~ dnorm(prior.occ[j,1], prior.occ[j,2])",
"}",
"for (j in 1:num.cov.det) {",
"theta[j] ~ dnorm(prior.det[j,1], prior.det[j,2])",
"}",
"}")
## evaluation starts here
method <- match.arg(method)
# X <- occ
# Z <- det
# Y <- obs
observations <- Y
N.sites <- length(Y)
num.cov.occ <- NCOL(X)
num.cov.det <- NCOL(Z)
linkinvfun.occ <- binomial(link=make.link(link.sta))$linkinv
linkinvfun.det <- binomial(link=make.link(link.det))$linkinv
## GLM estimates
glm.occ <- glm.fit(X, Y, family = binomial(link.sta))
glm.det <- glm.fit(Z, Y, family = binomial(link.det))
coef.occ <- coef(glm.occ)
coef.occ[is.na(coef.occ)] <- 0
coef.det <- coef(glm.det)
coef.det[is.na(coef.det)] <- 0
control.optim <- getOption("detect.optim.control")
## logLik fun must be maximized
control.optim$fnscale <- -1
opmeth <- getOption("detect.optim.method")
if (is.null(opmeth))
opmeth <- "BFGS"
control.dc <- if (method=="dc")
getOption("detect.dc.control") else NULL
control.mcmc <- if (method=="dc")
getOption("detect.mcmc.control") else NULL
Control <- list(optim.control=control.optim,
optim.method=opmeth,
dc.control=control.dc,
mcmc.control=control.mcmc)
pmle.problem <- FALSE
converged <- c(mle=NA, pmle=NA)
if (missing(inits)) {
# inits <- if (method=="dc")
# NULL else c(coef.occ, coef.det)
inits <- c(coef.occ, coef.det)
}
# if (inherits(inits, "svisit"))
# inits <- coef(inits)
inits <- rep(inits, length(c(coef.occ, coef.det)))[1:length(c(coef.occ, coef.det))]
## MLE from MCMC
if (method=="dc") {
prec <- control.dc$prec
n.clones <- control.dc$n.clones
## prior specifications
prior.occ <- cbind(coef.occ, rep(prec, length(coef.occ)))
prior.det <- cbind(coef.det, rep(prec, length(coef.det)))
excl.occ <- switch(link.sta,
"logit"=c(6,8),
"probit"=c(4,8),
"cloglog"=c(4,6))
excl.det <- switch(link.det,
"logit"=c(7,9),
"probit"=c(5,9),
"cloglog"=c(5,7))
model <- dclone::custommodel(mcmcSS.all, c(excl.occ, excl.det))
dat <- list(Y=obs, X=occ, Z=det, k=1,
N.sites=N.sites, num.cov.occ=num.cov.occ, num.cov.det=num.cov.det,
prior.occ=prior.occ, prior.det=prior.det)
# dat <- dclone(dat, n.clones, multiply=c("Y","k"), unchanged=names(dat)[-c(1,4)])
## old specification, depending on options
## new specification
# mle.res <- jags.engine(dat, c("beta", "theta"), model, inits=NULL, cl=cl, ...)
# mle.res <- jags.fit(dat, c("beta", "theta"), model, inits, ...)
## latest preference is that it depends on dcmle
requireNamespace("dcmle")
#requireNamespace("dclone")
dcd <- dcmle::makeDcFit(model=model, data=dat, params=c("beta", "theta"),
multiply=c("Y","k"), unchanged=names(dat)[-c(1,4)])
mle.res <- dcmle::dcmle(dcd, n.clones=n.clones, ...)
}
## MLE from optim
if (method=="optim") {
mle.res <- optim(inits, singleocc.MLE, gr = NULL, method = opmeth, hessian = TRUE,
control = control.optim,
observations = observations, X = X, Z = Z, link1 = link.sta, link2 = link.det)
if (rcond(mle.res$hessian) < 1e-06 && penalized) {
warning("cannot calculate PMLE, condition number too small")
penalized <- FALSE
}
}
## PMLE
if (penalized) {
if (method=="dc") {
#mle.parameters <- coef(mle.res)
#vv <- dclone::dcsd(mle.res)^2
vv <- dcmle::dcsd(mle.res)^2
mle.parameters <- dcmle::coef(mle.res)
}
if (method=="optim") {
mle.parameters <- mle.res$par
vv <- -diag(.solvenear(mle.res$hessian))
}
mle.par.occ <- mle.parameters[1:num.cov.occ]
mle.par.det <- mle.parameters[-(1:num.cov.occ)]
mle.phi <- linkinvfun.occ(X %*% mle.par.occ)
mle.delta <- linkinvfun.det(Z %*% mle.par.det)
## sometimes negative variance values occurr (strange but true)
## in this case MLE is returned, not PMLE, but attr is there to remember this unintended change
if (any(vv < 0)) {
pmle.res <- NULL
result <- mle.res
penalized <- FALSE
pmle.problem <- TRUE
warning("negative variance values in optim, penalization not used")
} else {
vv.occ <- sqrt(sum(vv[1:num.cov.occ]))
vv.det <- sqrt(sum(vv[-(1:num.cov.occ)]))
# parameters, observations, X, Z, bmax, dmax, psi0, d0, link1, link2, v1, v2
# v1:Square root of the sum of the estimated variance for the detection parameters obtained by the MLE
# v2:Square root of the sum of the estimated variance for the occupancy parameters obtained by the MLE
pmle.inits <- mle.parameters
pmle.res <- optim(pmle.inits, singleocc.PMLE, gr = NULL, method = opmeth, hessian = TRUE,
control = control.optim,
observations = observations, X = X, Z = Z,
bmax = coef.occ, dmax = coef.det,
psi0 = mean(mle.phi), d0 = mean(mle.delta),
link1=link.sta, link2=link.det,
v1 = vv.occ, v2 = vv.det)
result <- pmle.res
}
} else {
pmle.res <- NULL
result <- mle.res
}
## final evaluation of results
if (method=="dc" && !penalized) {
#std.error <- dclone::dcsd(mle.res)
#parameters <- coef(mle.res)
std.error <- dcmle::dcsd(mle.res)
parameters <- dcmle::coef(mle.res)
loglik <- singleocc.MLE(parameters, observations, X, Z, link1 = link.sta, link2 = link.det)
} else {
if (rcond(result$hessian) <= 1e-06)
std.error <- rep(NA, num.cov.occ + num.cov.det)
if (rcond(result$hessian) > 1e-06) {
opvar <- -diag(.solvenear(result$hessian))
if (any(opvar < 0)) {
opvar[opvar < 0] <- NA
## if pmle.problem=TRUE, this warning would be repeated
if (!pmle.problem)
warning("negative variance values in optim, NAs produced")
}
std.error <- sqrt(opvar)
}
parameters <- result$par
loglik <- result$value
}
converged[1] <- if (method=="dc")
dcmle::gelman.diag(mle.res)$mpsrf < 1.1 else mle.res$convergence == 0
if (penalized)
converged[2] <- pmle.res$convergence == 0
if (any(converged[!is.na(converged)]==FALSE))
warning("model did not converge\n")
## separating occ and det results
par.occ <- parameters[1:num.cov.occ]
par.det <- parameters[-(1:num.cov.occ)]
phi <- linkinvfun.occ(X %*% par.occ)
delta <- linkinvfun.det(Z %*% par.det)
se.occ <- std.error[1:num.cov.occ]
se.det <- std.error[-(1:num.cov.occ)]
temp <- (phi * (1 - delta))/((1 - phi) + (phi * (1 - delta)))
phi.obs <- ((observations == 1) * 1) + ((observations == 0) * temp)
names(par.occ) <- names(se.occ) <- colnames(X)
names(par.det) <- names(se.det) <- colnames(Z)
## AUC calculation
if (auc) {
zp <- (phi * delta)[observations == 0]
op <- (phi * delta)[observations == 1]
u <- sum(unlist(lapply(op, function(z) sum(zp < z))))
auc.out <- u / (sum(observations) * (length(observations) - sum(observations)))
}
## SE is not correct when using PMLE
if (penalized) {
se.occ[] <- NA
se.det[] <- NA
}
## assembling return object
out <- list(coefficients = list(sta = par.occ, det = par.det),
std.error = list(sta = se.occ, det = se.det),
estimated.probabilities = as.numeric(temp),
fitted.values = as.numeric(phi.obs),
method = method,
link = list(sta=link.sta, det=link.det),
nobs = N.sites,
df.null = N.sites - 2,
df.residual = N.sites - num.cov.occ - num.cov.det,
inits = inits,
# phi = as.numeric(phi),
occurrence.probabilities = as.numeric(phi),
# delta = as.numeric(delta),
detection.probabilities = as.numeric(delta),
loglik = loglik,
# results = list(glm.sta = glm.occ, glm.det = glm.det, mle = mle.res, pmle = pmle.res),
results = list(mle = mle.res, pmle = pmle.res),
converged = converged,
penalized = penalized,
control = Control)
# y = observations)
if (auc)
out$auc <- auc.out
Conv <- if (!out$penalized)
out$converged[1] else out$converged[2]
if (!Conv)
cat("Warning:\n Model did not converge\n\n")
return(out)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.