sightability_analysis: Example Sightability Data Analysis
In jlaake/DoubleObs: Analysis of double-observer survey data

An example of sightability data analysis using marked individuals to build a logistic regression model and then estimating abundance from those seen.

A data frame with 129 observations of 7 variables:

ch: a character vector containing the three character encounter history; first occasion indicaes whther the group was marked or not
count: number of animals in the group
het: variable used to model heterogeneity (re-capture probability.
lngs: natural log of count
veg: random vegetation variable
type: "unmarked" or "marked"

Data were generated with the following code using simhet function.

simres=simhet(1,150,25,beta=c(-.5,-.5,-1,2,.5),data=data.frame(veg=runif(175)), gen_formula="~-1+loc+veg+cov+lngs",use_gs=TRUE,gs_param=c(1,10),seed=9234397) example_data=simres[[1]][[1]]$results[[1]]$data$data

The groups sizes were randomly generated from a gamma distribution with scale=10 and shape=1 (gs_param). The true simulated population size is 175 groups and 1732 animals. A uniform random veg variable, and unmodelled heterogeneity variable cov (uniform -1,1) and natural log of group size are used to generate p for each location(loc).

data(example_data)
df=example_data
# Only use "sum" of second and third capture history positions; first position for radio(marked)
df$seen=ifelse((substr(df$ch,2,2)=="1") | (substr(df$ch,3,3)=="1"),1,0)
# For logistic model development only use marked groups.
lrdata=df[df$type=="marked",]
# restrict df to only those seen for prediction
df=df[df$seen==1,]
# Only fit a single model for this example
lrmodel=glm(seen~veg+lngs,family="binomial",data=lrdata)
# get predicted value of p and std error for each survey observation
pred=predict(lrmodel,df,type="response",se.fit=TRUE)
pdot=pred$fit
# construct variance-covariance matrix of predictions p
vcov=summary.glm(lrmodel)$cov.unscaled
dm=model.matrix(~veg+lngs,data=df)
deriv=dm*pdot*(1-pdot)
vc_pdot=deriv%*%vcov%*%t(deriv)
# construct variance
deriv=matrix(-1/pdot^2,nrow=1)
#' # Group abundance
Ngrp=sum(1/pdot)
se_Ngrp=sqrt(sum((1-pdot)/pdot^2)+deriv%*%vc_pdot%*%t(deriv))
Mt1=nrow(df)
Ngrp=Ngrp-Mt1
cv_Ngrp=se_Ngrp/Ngrp
C=exp(1.96*sqrt(log(1+cv_Ngrp^2)))
Ngrp_lcl=Ngrp/C+Mt1
Ngrp_ucl=Ngrp*C+Mt1
Ngrp=Ngrp+Mt1
paste("Ngrp = ",round(Ngrp), " se = ", round(se_Ngrp), " CI = (",round(Ngrp_lcl),",
",round(Ngrp_ucl),")",sep="")

# Animal abundance
deriv=matrix(-df$count/pdot^2,nrow=1)
N=sum(df$count/pdot)
se_N=sqrt(sum(df$count^2*(1-pdot)/pdot^2)+deriv%*%vc_pdot%*%t(deriv))
Mt1=sum(df$count)
N=N-Mt1
cv_N=se_N/N
C=exp(1.96*sqrt(log(1+cv_N^2)))
N_lcl=N/C+Mt1
N_ucl=N*C+Mt1
N=N+Mt1
paste("N = ",round(N), " se = ", round(se_N), " CI = (",round(N_lcl),",",round(N_ucl),")",sep="")