Nothing
#' Mark-Recapture Distance Sampling (MRDS) IO - FI
#'
#' Mark-Recapture Analysis of Independent Observer Configuration with Full
#' Independence
#'
#' The mark-recapture data derived from an independent observer distance
#' sampling survey can be used to derive conditional detection functions
#' (p_j(y)) for both observers (j=1,2). They are conditional detection
#' functions because detection probability for observer j is based on seeing or
#' not seeing observations made by observer 3-j. Thus, p_1(y) is estimated by
#' p_1|2(y).
#'
#' If detections by the observers are independent (full
#' independence) then p_1(y)=p_1|2(y),p_2(y)=p_2|1(y) and for the union, full
#' independence means that p(y)=p_1(y) + p_2(y) - p_1(y)*p_2(y) for each
#' distance y. In fitting the detection functions the likelihood given by eq
#' 6.8 and 6.16 in Laake and Borchers (2004) is used. That analysis does not
#' require the usual distance sampling assumption that perpendicular distances
#' are uniformly distributed based on line placement that is random relative to
#' animal distribution. However, that assumption is used in computing
#' predicted detection probability which is averaged based on a uniform
#' distribution (see eq 6.11 of Laake and Borchers 2004).
#'
#' For a complete description of each of the calling arguments, see
#' \code{\link{ddf}}. The argument \code{model} in this function is the same
#' as \code{mrmodel} in \code{ddf}. The argument \code{dataname} is the name
#' of the dataframe specified by the argument \code{data} in \code{ddf}. The
#' arguments \code{control},\code{meta.data},and \code{method} are defined the
#' same as in \code{ddf}.
#'
#' @export
#' @method ddf io.fi
#' @param dsmodel not used
#' @param mrmodel mark-recapture model specification
#' @param data analysis dataframe
#' @param meta.data list containing settings controlling data structure
#' @param control list containing settings controlling model fitting
#' @param call original function call used to call \code{ddf}
#' @param method analysis method; only needed if this function called from
#' \code{ddf.io}
#' @return result: an io.fi model object
#' @author Jeff Laake
#' @seealso
#' \code{\link{ddf.io}},\code{\link{summary.io.fi}},\code{\link{coef.io.fi}},
#' \code{\link{plot.io.fi}},\code{\link{gof.io.fi}},\code{\link{io.glm}}
#' @references Laake, J.L. and D.L. Borchers. 2004. Methods for incomplete
#' detection at distance zero. In: Advanced Distance Sampling, eds. S.T.
#' Buckland, D.R.Anderson, K.P. Burnham, J.L. Laake, D.L. Borchers, and L.
#' Thomas. Oxford University Press.
#' @keywords Statistical Models
#' @importFrom stats optimHess
#' @importFrom methods is
#'
ddf.io.fi <- function(dsmodel = NULL, mrmodel, data,
method, meta.data=list(), control=list(),
call=""){
#Rename arguments to match generic ddf
model <- mrmodel
# Functions used: assign.default.values, process.data, create.model.frame
# ioglm, predict(predict.io.fi), NCovered (NCovered.io.fi)
# NOTE: gams are only partially implemented
# The following are dummy glm and gam functions that are defined here to
# provide the list of arguments for use in the real glm/gam functions.
# These dummy functions are removed after they are used so the real ones can
# be used in the model fitting.
glm <- function(formula, link="logit"){
formula <- fixformula(formula)
if(is(link, "function")){
link <- substitute(link)
}
link <- match.arg(link, c("logit"))
return(list(fct="glm", formula=formula, link=substitute(link)))
}
gam <- function(formula, link="logit"){
formula <- fixformula(formula)
if(is(link, "function")){
link <- substitute(link)
}
link <- match.arg(link, c("logit"))
return(list(fct="gam", formula=formula, link=substitute(link)))
}
# Save current user options and then set design contrasts to treatment style
save.options<-options()
options(contrasts=c("contr.treatment", "contr.poly"))
# Set up meta data values
meta.data <- assign.default.values(meta.data, left=0, width=NA, binned=FALSE,
int.range=NA, point=FALSE)
# Set up control values
control <- assign.default.values(control, showit=0,
estimate=TRUE, refit=TRUE, nrefits=25,
initial=NA, lowerbounds=NA,
upperbounds=NA, mono.points=20)
# Assign model values; this uses temporarily defined functions glm and gam
modpaste <- paste(model)
modelvalues <- try(eval(parse(text=modpaste[2:length(modpaste)])))
if(inherits(modelvalues, "try-error")){
stop("Invalid model specification: ",model)
}
rm(glm,gam)
# Process data if needed
if(is.data.frame(data)){
data.list <- process.data(data, meta.data)
meta.data <- data.list$meta.data
xmat <- data.list$xmat
}else{
xmat <- data
}
# Setup default breaks
if(meta.data$binned){
meta.data$breaks <- c(max(0,
min(as.numeric(levels(as.factor(xmat$distbegin))))
),
as.numeric(levels(as.factor(xmat$distend))))
}
# Create result list with some arguments
result <- list(call=call, data=data, model=model, meta.data=meta.data,
control=control, method="io.fi")
class(result) <- c("io.fi", "ddf")
# Create formula and model frame (if not GAM)
xmat$offsetvalue <- rep(0, dim(xmat)[1])
model.formula <- paste("detected", modelvalues$formula)
p.formula <- as.formula(model.formula)
xmat2 <- xmat[xmat$observer==2,]
xmat1 <- xmat[xmat$observer==1,]
GAM <- FALSE
xmat <- create.model.frame(xmat, as.formula(model.formula), meta.data)
model.formula <- as.formula(paste(model.formula, "+offset(offsetvalue)"))
# Fit the conditional detection functions using io.glm
suppressWarnings(result$mr <- io.glm(xmat, model.formula, GAM=GAM))
# if(GAM)result$mr$data=xmat
if(result$mr$converged){
# if the glm did converge, then do one quick round of BFGS to
# compute the hessian
xx <- rbind(xmat1,xmat2)
dmrows <- nrow(xmat1)
xmatH <- model.matrix(p.formula, xx)
xmat1H <- xmatH[1:dmrows, , drop=FALSE]
xmat2H <- xmatH[(dmrows+1):(2*dmrows), , drop=FALSE]
result$hessian<- optimHess(result$mr$coefficients, lnl.io,
x1=xmat1H, x2=xmat2H,
x1_detected=xmat1$detected,
x2_detected=xmat2$detected,
models=list(p.formula=p.formula))
}else{
# only resort to optim if there was not convergence in the glm!
# Now use optimx with starting values perturbed by 5%
fit <-try(optimx(result$mr$coefficients*1.05,
lnl.io, method="L-BFGS-B", hessian=TRUE,
x1=xmat1H, x2=xmat2H,
x1o=xmat1, x2o=xmat2,
models=list(p.formula=p.formula)))
# did this model converge?
topfit.par <- coef(fit, order="value")[1, ]
details <- attr(fit, "details")[1,]
fit <- as.list(summary(fit, order="value")[1, ])
fit$par <- topfit.par
fit$message <- ""
names(fit)[names(fit)=="convcode"] <- "conv"
fit$hessian<-details$nhatend
if(fit$conv!=0 | inherits(fit, "try-error")){
# first try the old way of just setting the starting values to zero,
# this seems (with the crabbie data) to converge back to the values in
# result$mr$coefficients but without having the convergence issues
# NEED TO THINK ABOUT THIS MORE...
fit <- try(optimx(result$mr$coefficients*0, lnl.io, method="L-BFGS-B",
hessian=TRUE,
x1=xmat1H, x2=xmat2H,
x1o=xmat1, x2o=xmat2,
models=list(p.formula=p.formula)))
topfit.par <- coef(fit, order="value")[1, ]
details <- attr(fit,"details")[1,]
fit <- as.list(summary(fit, order="value")[1, ])
fit$par <- topfit.par
fit$message <- ""
names(fit)[names(fit)=="convcode"] <- "conv"
fit$hessian <- details$nhatend
}
# if nothing worked...
if(fit$conv!=0 | inherits(fit, "try-error")){
stop("No convergence in ddf.io.fi()")
}else{
result$mr$mr$coefficients <- fit$par
result$hessian <- fit$hessian
}
}
# Compute the L_omega portion of the likelihood value, AIC and hessian
cond.det <- predict(result)
result$par <- coef(result$mr)
npar <- length(result$par)
p1 <- cond.det$p1
p2 <- cond.det$p2
p.c.omega <- p1^result$mr$data$detected[result$mr$data$observer==1]*
(1-p1)^(1-result$mr$data$detected[result$mr$data$observer==1])*
p2^result$mr$data$detected[result$mr$data$observer==2]*
(1-p2)^(1-result$mr$data$detected[result$mr$data$observer==2])
result$lnl <- sum(log(p.c.omega)) - sum(log(cond.det$fitted))
#if(GAM)
# result$hessian <- result$mr$Vp
#else
# result$hessian <- solve(summary(result$mr)$cov.unscaled)
# If this is method=io.fi then compute lnl for L_y and add to L_omega before
# computing AIC. Note the code in predict is currently specific to logit link
if(method=="io.fi"){
# Compute fitted values - integral of p(x)/width
result$fitted <- predict(result,integrate=TRUE)$fitted
result$Nhat <- NCovered(result$par,result)
distances <- result$data$distance[result$data$observer==1]
if(!meta.data$binned){
if(meta.data$point){
result$lnl <- result$lnl +
sum(log(cond.det$fitted*2*distances/meta.data$width^2)) -
sum(log(result$fitted))
}else{
result$lnl<- result$lnl +
sum(log(cond.det$fitted/meta.data$width)) -
sum(log(result$fitted))
}
}else{
for(i in 1:(nrow(xmat)/2)){
int.val <- predict(result, newdata=xmat[(2*(i-1)+1):(2*i),],
int.range=as.vector(as.matrix(xmat[(2*(i-1)+1),
c("distbegin", "distend")])),
integrate=TRUE)$fitted
result$lnl <- result$lnl + log(int.val)
}
result$lnl<- result$lnl- sum(log(result$fitted))
}
}
result$criterion <- -2*result$lnl + 2*npar
# Restore user options
options(save.options)
# Return result
return(result)
}
io.p01 <- function(xmat1, xmat2, beta){
p1 <- plogis(xmat1%*%beta)
p2 <- plogis(xmat2%*%beta)
return(p2*(1-p1))
}
io.p10 <- function(xmat1, xmat2, beta){
p1 <- plogis(xmat1%*%beta)
p2 <- plogis(xmat2%*%beta)
return(p1*(1-p2))
}
io.p11 <- function(xmat1, xmat2, beta){
p1 <- plogis(xmat1%*%beta)
p2 <- plogis(xmat2%*%beta)
return(p1*p2)
}
io.pdot<-function(xmat1, xmat2, beta){
p1 <- plogis(xmat1%*%beta)
p2 <- plogis(xmat2%*%beta)
return(p1+p2-p1*p2)
}
p.io <- function(par, x1, x2, models){
# # create design matrix for p1,p2 and delta
# x <- rbind(x1,x2)
# dmrows <- nrow(x1)
# xmat <- model.matrix(models$p.formula,x)
# xmat1 <- xmat[1:dmrows,,drop=FALSE]
# xmat2 <- xmat[(dmrows+1):(2*dmrows),,drop=FALSE]
xmat1 <- x1
xmat2 <- x2
p01 <- io.p01(xmat1, xmat2, beta=par)
p10 <- io.p10(xmat1, xmat2, beta=par)
p11 <- io.p11(xmat1, xmat2, beta=par)
pdot <- io.pdot(xmat1, xmat2, beta=par)
return(list(p11=p11,
p01=p01,
p10=p10,
pdot=pdot))
}
lnl.io <- function(par, x1, x2, models, x1_detected, x2_detected){
# Compute probabilities
p.list <- p.io(par,x1,x2,models)
p11 <- p.list$p11
p01 <- p.list$p01
p10 <- p.list$p10
p01[p01==0] <- 1e-6
p10[p10==0] <- 1e-6
# Compute log-likelihood value
lnl <- sum((1-x1_detected)*x2_detected*log(p01)) +
sum((1-x2_detected)*x1_detected*log(p10))+
sum(x1_detected*x2_detected*log(p11)) -
sum(log(p.list$pdot))
return(-lnl)
}
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.