# R/covtheta.R In SightabilityModel: Wildlife Sightability Modeling

#### Documented in covtheta

#' Estimates var/cov matrix of inflation factors (1/prob detection) using a
#' non-parametric bootstrap.
#'
#' Estimates var/cov matrix of inflation factors (1/prob detection) using a
#' non-parametric bootstrap.  Called by function Sight.Est if Vm.boot = TRUE.
#'
#'
#' @param total Number of animals in each independently sighted group
#' @param srates Plot sampling probability (associated with the independently
#' observed animal groups)
#' @param stratum Stratum identifiers (associated with the independently
#' observed animal groups)
#' @param subunit Plot ID (associated with the independently observed animal
#' groups)
#' @param covars Matrix of sightability covariates (associated with the
#' independently observed animal groups)
#' @param betas Logistic regression parameter estimates (from fitted
#' sightability model)
#' @param varbetas Estimated variance-covariance matrix for the logistic
#' regression parameter estimates (from fitted sightability model)
#' @param nboots Number of bootstrap resamples.
#' @return \item{smat }{Estimated variance-covariance matrix for the inflation
#' factors theta = (1/probability of detection).  This is an n.animal x
#' n.animal matrix. }
#' @author John Fieberg
#' @keywords methods
#' @export covtheta
covtheta <-
function(total, srates, stratum, subunit, covars, betas, varbetas, nboots){

#  Arguments (the first four should have one observation for each individual in the obs. dataset
#  total = total number of animals in the observed group
#  srates = sampling rate for the stratum (for each observed animal)
#  subunit = subunit identifier
#  covars = matrix of covariates used in the sightability model
#  beta = vector of parameter estimates from logistic model
#  varbeta = var/cov matrix for above parameter estimates
# nboots = number of bootstrap replicates

# Check on length of vector arguments
n <- length(total)
if(length(srates) != n) {stop("Srates vector needs to be same dimension as total vector")}
if(length(stratum) != n) {stop("Stratum vector needs to be same dimension as total vector")}
if(length(srates) != n) {stop("Subunit vector needs to be same dimension as total vector")}

# Make sure strata are numbered 1 to h
stratum <- as.numeric(factor(stratum))

# Form xmatrix
xdata <- as.matrix(cbind(rep(1, n), covars))
ncovars <- dim(xdata)
nxdata <- nrow(xdata)
if(ncovars != n) {stop("Covariate matrix needs to be same length as total vector")}

# Make sure beta is a matrix of dim ncovar x 1
thetas <- matrix(NA, nboots, n)
for(j in 1:nboots){
tempbet <- betas[j, ]
tempvarbeta <- varbetas[j, , ]

#  Estimate theta for each animal = 1/prob(detection)
thetas[j, ] <- 1+exp(-xdata%*%tempbet-diag(xdata%*%tempvarbeta%*%t(xdata))/2)
}

# Use empirical variance for Cov(thetas)
smat <- cov(thetas)
return(smat)
}


## Try the SightabilityModel package in your browser

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

SightabilityModel documentation built on May 20, 2022, 9:06 a.m.