I'm using the pop.est function to produce abundance estimates from a CJS model. These estimates are used to compare abundance estimates from other models and are not the main focus of the analysis.
A few questions:
1) I'm struggling with the specification of the input values. Specifically, my data are grouped using 2 variables, and I want to get occasion-specific abundance for each of the 4 groups, as well as summed up across groups within[/list] occasion. The help file says
If individual estimates are needed, use an nxn identity matrix for design where n is the length of ns. To get a total of all the estimates use a nx1 column matrix of 1s. Any other design matrix can be specified to subset, aggregate and/or average the estimates
I managed to get summation across occasions (using the n x 1 matrix), but can't figure out how to sum across groups within occasion.
2) 95% CIs - in the paper cited for pop.est (https://www.bearbiology.com/wp-content/uploads/2017/10/Taylor_13.pdf), the authors produced non-symmetric 95% CIs. How do I get from the vcv matrix provided by pop.est to the CIs?
3) Currently, I'm using an exceedingly clunky way to get the group-specific estimates. Is there a more succinct way to get those? This goes together with Q1 above, since either the clunky estimates need to be summed up, or they need to be less clunky to start with.
4) Is there a way to use this following model averaging? The examples I found were all for a single model.
Thanks so much!
example code: 1) load libraries, data, create grouping vars, and process data, run a CJS model
- Code: Select all
library(plyr)
library(dplyr)
library(RMark)
data(dipper)
# make up another grouping variable
dipper$Var1 <- sample(c("A", "B"), nrow(dipper), replace = TRUE)
data1.proc <- process.data(dipper, model="CJS", begin.time = 2000, groups = c("sex", "Var1"))
data1.ddl <- make.design.data(data1.proc)
data1.analysis=function(){
Phi.1 = list(formula= ~ 1)
p.t = list(formula=~ time)
cml=create.model.list("CJS")
mark.wrapper(cml,data=data1.proc, ddl=data1.ddl, output=FALSE, delete = FALSE)
}
data1.results <- data1.analysis()
2) clunky group-specific pop.est
- Code: Select all
md <- data1.results[[1]]
k <- 1
pop <- list()
for(sex in unique(dipper$sex)){
for(var1 in unique(dipper$Var1)){
sub <- dipper[dipper$sex == sex & dipper$Var1 == var1,]
xmat <- matrix(as.numeric(unlist(strsplit(sub$ch,""))), ncol=nchar(sub$ch[1]))
ns <- colSums(xmat)[-1]
lab <- paste("sex", sex, sep = "")
lab1 <- paste("Var1", var1, sep = "")
lab <- paste(lab, lab1, sep = ".")
group <- which(md$group.labels == lab)
p.indices <- extract.indices(md, "p", df = data.frame(group = rep(group,6), row=1:6,col=1:6))
p.list <- covariate.predictions(md,data=data.frame(index=p.indices))
# call pop.est using diagonal design matrix to get
# separate estimate for each occasion
pop[[k]] <- as.data.frame(pop.est(ns,p.list$estimates$estimate, design=diag(1, ncol=6, nrow=6),p.list$vcv)$Nhat)
pop[[k]]$sex <- sex
pop[[k]]$Var1 <- var1
k <- k + 1
}
}
pop <- ldply(pop, data.frame)
apply(pop[, 1:6], 2, sum)
3) total estimates - as per the help files. The values above aren't the same as the totals here
- Code: Select all
xmat=matrix(as.numeric(unlist(strsplit(dipper$ch,""))), ncol=nchar(dipper$ch[1]))
ns=colSums(xmat)[-1]
p.indices=extract.indices(md,"p",df=data.frame(group=rep(1,6), row=1:6,col=1:6))
p.list=covariate.predictions(md,data=data.frame(index=p.indices))
pop.est(ns,p.list$estimates$estimate, design=diag(1,ncol=6,nrow=6),p.list$vcv)