############################################################################################## # R code to conduct simulations in program CAPWIRE - Pennell et al. 2013, Miller et al. 2008 # ############################################################################################## # Install required package library(capwire) # Conduct abundance estimate for target population (i.e. MR black bears using 3 hair sampling sessions) #Create vector of capture histories MRBear.CapHis <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, 1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3, 4,4,4,5,6,6) MRBear.CapHis # Build CAPWIRE input file MRBear.Table <- buildClassTable(MRBear.CapHis) MRBear.Table # Fit equal capture rates model (ecm) MRBear.ecm <- fitEcm(data=MRBear.Table, max.pop=200) MRBear.ecm # Fit two innate rates model (tirm) MRBear.tirm <- fitTirm(data=MRBear.Table, max.pop=200) MRBear.tirm # Perform likelihood ratio test to determine the best supported model MRBear.lik.ratio <- lrtCapwire(ecm=MRBear.ecm, tirm=MRBear.tirm, bootstraps=100) MRBear.lik.ratio # Calculate confidence interval for tirm MRBear.CI.tirm <- bootstrapCapwire(x=MRBear.tirm, bootstraps=1000, CI=c(0.025, 0.975)) MRBear.CI.tirm # Perform simulations to evaluate the precision of abundance estimates for different sample sizes # Create 1000 simulated capture histories using the parameters estimated via the better supported # tirm and your desired sample size (ml.na, ml.nb, alpha, sample size) holder <- list() for(i in 1:1000){ sim <- simTirm(25,110,5.338537,250) holder[[i]] <- sim } # Fit tirm to the simulated capture histories model <- list() for(i in 1:1000){ simTirm <- fitTirm(data <- holder[[i]], max.pop=200) model[[i]] <- simTirm } # Calculate confidence intervals (CI) for each simulated capture history CI <- list() for(i in 1:1000){ simTirm.boot <- bootstrapCapwire(x=model[[i]], bootstraps=1000, CI=c(0.025, 0.975)) CI[[i]] <- simTirm.boot } # Extract and calculate means across abundance and CI estimates for each simulated capture histories n <- list() lowCI <- list() highCI <- list() for(i in 1:1000){ runs <- CI[[i]] pop.size <- runs$ml.pop.size ci <- runs$conf.int low <- ci[1] high <- ci[2] n[[i]] <- pop.size lowCI[[i]] <- low highCI[[i]] <- high } mean.n <- mean(as.numeric(n)) mean.lowCI <- mean(as.numeric(lowCI)) mean.highCI <- mean(as.numeric(highCI)) MRBear.sim <- cbind(mean.n,mean.lowCI,mean.highCI) MRBear.sim