#'This function asseses some frequentist properties when naively fitting Mb to capture histories that are missing data. Models
#'are fit in Mark via Rmark
#'@param N an integer that is the true population size
#'@param p a vector of size 2 containing the capture and recapture probabilities
#'@param occ an integer that is the number of capture occasions to simulate
#'@param K an integer that is the number of traps to simulate
#'@param traptype a character string specifying the trap type. "single"
#'allows individuals to be caught in at most 1 trap per occasion and "multi" allows
#'individuals to be caught in multiple traps per occasion.
#'@param btype a character string specifying the type of behavioral response to capture
#'Current options are "global" for a global trap response and "trap" for a trap-specific
#'response. "trap" currently only works with traptype=multi.
#'@param lambda_h a positive value specifying the parameter for the zero-truncated Poisson
#'hair deposition process (S_ijk|W_ijk=1)
#'@param lambda_c a positive value specifying the parameter for the zero-truncated Poisson
#'cluster deposition process if cluster=TRUE
#'@param kappa_h an integer specifying the number of hair samples to retain per occasion-trap if cluster=FALSE
#'or per occasion-trap-cluster if cluster=TRUE. Used in sub-hair and sub-cluster methods
#'@param kappa_t an integer specifying the number of traps per occasion at which to retain kappa_h hair samples.
#'If no value is specified no trap-level subsampling occurs. Only works when cluster=FALSE. Used in the
#'sub-trap method
#'@param delta a numeric value between 0 and 1 specifying the probability a hair sample will be retained
#'in the subsample. This is used if subsampling is done by pooling samples on each occasion and taking
#'a simple random sample.
#'@param cluster a logical indicating whether or not to simulate from the cluster model (multiple clusters per
#'individual-occasion-trap)
#'@param kappa_c an integer specifying the number of clusters per trap-occasion at which to retain kappa_h hair samples.
#'Only works when cluster=TRUE. Sub-cluster method.
#'@param alpha a numeric value between 0 and 1 specifying the probability that a hair sample produces
#'an individual identification
#'@param sims the number of simulated data sets to generate and fit Mb to
#'@param cores the number of cores to do the simulation on. A typical for loop is used if cores=1 and
#'the foreach package is used if cores>1
#'@return a list with the simulation results (expand later)
#'@export
assess=function(N,p,occ,K,traptype,btype="global",lambda_h,lambda_c=NULL,delta=NULL,kappa_h=NULL,
kappa_t=NULL,cluster=FALSE,kappa_c=NULL,alpha=NULL,sims=100,cores=1){
if(cores==1){
#preallocate storage structures
storeN=matrix(NA,nrow=sims,ncol=3)
storep=rep(NA,sims)
storec=rep(NA,sims)
storeFail=rep(NA,sims)
storeCapsTrue=matrix(NA,nrow=sims,ncol=occ)
storeCapsObs=matrix(NA,nrow=sims,ncol=occ)
storeFirstcapstrue=matrix(NA,nrow=sims,ncol=occ)
storeMiss1caps=matrix(NA,nrow=sims,ncol=occ)
storeMdottrue=rep(NA,sims)
storeMdotobs=rep(NA,sims)
storeMt1true=rep(NA,sims)
storeMt1obs=rep(NA,sims)
storeS=matrix(NA,nrow=sims,ncol=occ)
storeU=matrix(NA,nrow=sims,ncol=occ)
storeR=matrix(NA,nrow=sims,ncol=occ)
storepstats=array(NA,dim=c(5,occ,sims))
for(i in 1:sims){
data=simMBmisscap(N,p,occ,K,traptype,btype,lambda_h,lambda_c,delta,kappa_h,
kappa_t,cluster,kappa_c,alpha)
Mdata=as.data.frame(apply(data$Wobs2D,1,paste,collapse=""))
colnames(Mdata)="ch"
Mdata$ch=as.character(Mdata$ch)
out=RMark::mark(Mdata,model="Closed",output=FALSE,delete=TRUE) #Mb
storeN[i,]=unlist(out$results$derived[c(1,3,4)])
storep[i]=RMark::get.real(out,"p")[1]
storec[i]=RMark::get.real(out,"c")[1]
storeFail[i]=data$Wstats$fail
storeCapsTrue[i,]=colSums(data$W2D)
storeCapsObs[i,]=colSums(data$Wobs2D)
storeFirstcapstrue[i,]=data$Wstats$firstcapsTrue
storeMiss1caps[i,]=data$Wstats$miss1caps
storeMdottrue[i]=data$Wstats$Mdot_true
storeMdotobs[i]= data$Wstats$Mdot_obs
storeMt1true[i]=data$Wstats$Mt1_true
storeMt1obs[i]=data$Wstats$Mt1_obs
storeS[i,]=data$HairStats$Sdotj
storeU[i,]=data$HairStats$Udotj
storeR[i,]=data$HairStats$Rdotj
storepstats[,,i]=data$pstats
}
den=sims
}else{
#load multicore packages
require(foreach)
require(snow)
require(doSNOW)
require(RMark)
cl.tmp = makeCluster(rep("localhost",cores), type="SOCK")
registerDoSNOW(cl.tmp)
#clusterExport(cl.tmp, varlist=c("alpha","btype","cluster","delta"), envir = environment())
clusterExport(cl.tmp, list=ls(environment()), envir = environment())
out2=foreach(i=1:sims,.packages=c("VGAM","RMark"),.export="simMBmisscap",.errorhandling = "remove") %dopar% {
data=simMBmisscap(N,p,occ,K,traptype,btype,lambda_h,lambda_c,delta,kappa_h,
kappa_t,cluster,kappa_c,alpha)
Mdata=as.data.frame(apply(data$Wobs2D,1,paste,collapse=""))
colnames(Mdata)="ch"
Mdata$ch=as.character(Mdata$ch)
out=mark(Mdata,model="Closed",output=FALSE,delete=FALSE,threads=1) #Mb
f0=get.real(out,"f0",se=TRUE)[3:4]
c=exp(1.96*sqrt((log(1+f0[2]^2/(f0[1]^2)))))
Nhat=unlist(c(data$n+f0[1],data$n+f0[1]/c,data$n+f0[1]*c))
#Nhat=out$results$derived[c(1,3,4)]
params=list(Nhat,
get.real(out,"p")[1],
get.real(out,"c")[1])
return(list(params,
data$Wstats$fail,
colSums(data$W2D),
colSums(data$Wobs2D),
data$Wstats$firstcapsTrue,
data$Wstats$miss1caps,
data$Wstats$Mdot_true,
data$Wstats$Mdot_obs,
data$Wstats$Mt1_true,
data$Wstats$Mt1_obs,
data$HairStats$Sdotj,
data$HairStats$Udotj,
data$HairStats$Rdotj,
data$pstats))
}
stopCluster(cl.tmp)
cleanup(ask=FALSE,lx=ls(.GlobalEnv))
den=length(out2) #successful fits RMark didn't loose track of
#preallocate storage structures
storeN=matrix(NA,nrow=den,ncol=3)
storep=rep(NA,den)
storec=rep(NA,den)
storeFail=rep(NA,den)
storeCapsTrue=matrix(NA,nrow=den,ncol=occ)
storeCapsObs=matrix(NA,nrow=den,ncol=occ)
storeFirstcapstrue=matrix(NA,nrow=den,ncol=occ)
storeMiss1caps=matrix(NA,nrow=den,ncol=occ)
storeMdottrue=rep(NA,den)
storeMdotobs=rep(NA,den)
storeMt1true=rep(NA,den)
storeMt1obs=rep(NA,den)
storeS=matrix(NA,nrow=den,ncol=occ)
storeU=matrix(NA,nrow=den,ncol=occ)
storeR=matrix(NA,nrow=den,ncol=occ)
storepstats=array(NA,dim=c(5,occ,den))
for(i in 1:den){
storeN[i,]=out2[[i]][[1]][[1]]
storep[i]=out2[[i]][[1]][[2]]
storec[i]=out2[[i]][[1]][[3]]
storeFail[i]=out2[[i]][[2]]
storeCapsTrue[i,]=out2[[i]][[3]]
storeCapsObs[i,]=out2[[i]][[4]]
storeFirstcapstrue[i,]=out2[[i]][[5]]
storeMiss1caps[i,]=out2[[i]][[6]]
storeMdottrue[i]=out2[[i]][[7]]
storeMdotobs[i]=out2[[i]][[8]]
storeMt1true[i]=out2[[i]][[9]]
storeMt1obs[i]=out2[[i]][[10]]
storeS[i,]=out2[[i]][[11]]
storeU[i,]=out2[[i]][[12]]
storeR[i,]=out2[[i]][[13]]
storepstats[,,i]=out2[[i]][[14]]
}
}
#Calcuate N stats
Nmean=mean(storeN[,1],na.rm=TRUE)
dir=Nmean[1]-N
Nbias=ifelse(dir>0,paste("+",round(100*(Nmean[1]/N-1),2),"%"),paste("-",round(100*(1-Nmean[1]/N),2),"%"))
Ncover=round(sum(storeN[,2]<N&storeN[,3]>N,na.rm=TRUE)/den,2)
Nwidth=round(mean(storeN[,3]-storeN[,2],na.rm=TRUE),2)
Nfail=sims-den
LikFails0=sum(storeFail<=0)
LikFails100=sum(storeFail<=100)
Nstats=data.frame(MeanNhat=round(Nmean,2),Bias=Nbias,Coverage=Ncover,CIwidth=Nwidth,
simFails=Nfail,LikFails0=LikFails0,LikFails100=LikFails100)
#Calculate capstats stats
CapStats=rbind(colMeans(storeCapsTrue),colMeans(storeCapsObs),colMeans(storeCapsTrue)-colMeans(storeCapsObs))
rownames(CapStats)=c("True Capture Events","Observed Capture Events","NMissing")
CapStats=cbind(CapStats,rowSums(CapStats))
colnames(CapStats)=c(paste("T",1:occ,sep=""),"Total")
CapStats=rbind(CapStats,CapStats[2,]/CapStats[1,])
rownames(CapStats)[4]="%remain"
CapStats=rbind(CapStats,c(colMeans(storeFirstcapstrue),sum(colMeans(storeFirstcapstrue))))
CapStats=rbind(CapStats,c(colMeans(storeFirstcapstrue-storeMiss1caps),sum(colMeans(storeFirstcapstrue-storeMiss1caps))))
CapStats=rbind(CapStats,CapStats[5,]-CapStats[6,])
CapStats=rbind(CapStats,CapStats[6,]/CapStats[5,])
CapStats=round(CapStats,2)
rownames(CapStats)[5:8]=c("True First Captures","Observed First Captures","NMissing","%remain")
#Calculate hair stats
HairStats=rbind(colMeans(storeS),colMeans(storeU),colMeans(storeR))
HairStats=cbind(HairStats,rowSums(HairStats))
colnames(HairStats)=c(paste("T",1:occ,sep=""),"Total")
HairStats=rbind(HairStats,HairStats[2,]/HairStats[1,],HairStats[3,]/HairStats[1,])
rownames(HairStats)=c("Sdotj","Udotj","Rdotj","%remainU","%remainR")
HairStats=round(HairStats,2)
#calculate p stats
pstats=apply(storepstats,1:2,mean)
rownames(pstats)=c("p_true","p_obs","p_obs_true","p_obs_false","percent_false")
colnames(pstats)=paste("T",1:occ,sep="")
pstats=round(pstats,2)
#Calculate Mb stats
MbStats=matrix(c(mean(storeMt1true),mean(storeMt1obs),mean(storeMdottrue),mean(storeMdotobs)),nrow=2,byrow=TRUE)
rownames(MbStats)=c("M_t+1","M.")
colnames(MbStats)=c("True","Observed")
return(list(Nstats=Nstats,CapStats=CapStats,HairStats=HairStats,MbStats=MbStats,pstats=pstats,storeN=storeN,storep=storep,storec=storec,storeFail=storeFail,
storeMdottrue=storeMdottrue,storeMdotobs=storeMdotobs,storeMt1true=storeMt1true,
storeMt1obs=storeMt1obs,storeS=storeS,storeU=storeU))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.