# R/varlog.lam.R In SightabilityModel: Wildlife Sightability Modeling

#### Documented in varlog.lam

#' Calculates the variance of the log rate of change between 2 population
#' estimates that rely on the same sightability model.
#'
#' Calculates the variance of the log rate of change between 2 population
#' estimates that rely on the same sightability model.
#'
#' This function uses the delta method to calculate an approximate variance for
#' the log rate of change, log(tau^[t+1])-log(tau^[t]), while accounting for
#' the positive covariance between the two estimates (as a result of using the
#' same sightability model to correct for detection).
#'
#' @param sight1 Sightability model object for the first population estimate
#' (formed by calling Sight.Est function)
#' @param sight2 Sightability model object for the second population estimate
#' (formed by calling Sight.Est function)
#' @return \item{loglambda}{log rate of change = log(tau^[t+1]/tau^[t])}
#' \item{varloglamda}{approximate variance of loglambda}
#' @author John Fieberg
#' @keywords methods
#' @examples
#'
#' # Example using moose survey data
#'   data(obs.m) # observational moose survey data
#'   data(exp.m) # experimental moose survey data
#'   data(sampinfo.m) # information on sampling rates
#'
#' # Estimate population size in 2006 and 2007
#'   sampinfo <- sampinfo.m[sampinfo.m$year==2007, ] #' tau.2007 <- Sight.Est(observed ~ voc, odat = obs.m[obs.m$year==2007, ],
#'                           sdat = exp.m, sampinfo.m[sampinfo.m$year == 2007, ], #' method = "Wong", logCI = TRUE, alpha = 0.05, Vm.boot = FALSE) #' tau.2006 <- Sight.Est(observed ~ voc, odat = obs.m[obs.m$year==2006, ],
#'                           sdat = exp.m, sampinfo.m[sampinfo.m$year == 2006, ], #' method = "Wong", logCI = TRUE, alpha = 0.05, Vm.boot = FALSE) #' #' # Log rate of change #' varlog.lam(tau.2006, tau.2007) #' #' @export varlog.lam varlog.lam <- function(sight1,sight2){ if(sight1$call$form != sight2$call$form){ stop("Need same sightability model to calculate the covariance") } if(is.null(sight1$sight$note)){ varbeta <- vcov(sight1$sight)
beta <- coef(sight1$sight) }else{ varbeta <- sight1$sight$varbet beta <- sight1$sight$bet } y1 <- sight1$odat$total y2 <- sight2$odat$total n1 <- nrow(sight1$odat)
n2 <- nrow(sight2$odat) inv.srate1 <- 1/sight1$odat$samp.rates inv.srate2 <- 1/sight2$odat$samp.rates fo <- sight1$call$form class(fo) <- "formula" tempnm1 <- terms(fo, data = sight1$odat)
tempnm2 <- attr(tempnm1, "term.labels")
covars1 <- sight1$odat[, tempnm2] covars2 <- sight2$odat[, tempnm2]

# Do as much of the matrix multiplication outside of loop as possible
xdat1 <- as.matrix(cbind(rep(1, n1), covars1))
xdat2 <- as.matrix(cbind(rep(1, n2), covars2))
xb1 <- xdat1%*%beta  # X'beta
xb2 <- xdat2%*%beta  # X'beta
xbb <- kronecker(xb1, t(xb2), FUN = "+")  # X1+x2
xvarbeta <- xdat1%*%varbeta%*%t(xdat2)  # X Sig X
smat <- matrix(0, n1, n2) # holder for some of the terms in the expression for smat
for(i in 1:n1){
for(j in 1:n2){
xtemp1 <- as.vector(xdat1[i, ], mode = "numeric")
xtemp2 <- as.vector(xdat2[j, ], mode = "numeric")
xtot <- t(xtemp1+xtemp2)
smat[i,j] <- (xtot%*%varbeta%*%t(xtot))/2
}
}

smat.cov <- exp(-xbb-smat)*(exp(xvarbeta)-1)
y.p1 <- as.matrix(y1*inv.srate1, n1, 1)
y.p2 <- as.matrix(y2*inv.srate2, n2, 1)

cov.total <- t(y.p1)%*%smat.cov%*%y.p2

var.tau1.tau2 <- matrix(c(sight2$est[2], cov.total, cov.total, sight1$est[2]), ncol = 2, byrow = TRUE)
dfs <- matrix(c(1/sight2$est[1], -1/sight1$est[1]), 1, 2)
varloglam <- list(loglamda = log(sight2$est[1]/sight1$est[1]), varloglamda = dfs%*%var.tau1.tau2%*%t(dfs))
return(varloglam = varloglam)
}


## Try the SightabilityModel package in your browser

Any scripts or data that you put into this service are public.

SightabilityModel documentation built on May 2, 2023, 1:12 a.m.