CJS abundance - pop.est from RMark

posts related to the RMark library, which may not be of general interest to users of 'classic' MARK

CJS abundance - pop.est from RMark

Postby SoConfused » Thu Mar 22, 2018 1:52 pm

Hello,

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)

SoConfused
 
Posts: 53
Joined: Wed Nov 05, 2014 8:25 am

Re: CJS abundance - pop.est from RMark

Postby jlaake » Mon Mar 26, 2018 2:36 pm

pop.est was just an extra that I threw in. If you look at the code it is not a big function so you are always welcome to do it however you want. Notice that in the example I provided there are no groups. Your clunky method tries to use the code for no groups and expand to groups and it is not correct. All of the p's have covariances so you can't work with each group independently and then expect to sum up the results and get the correct vcv matrix for the totals. For the dipper example you provided you'll need to get 24 (4 groups x 6 occasions) n's and p's and the 24x24 vcv matrix of the p's. Then to get all 24 estimates you would just use a 24x24 identity matrix. To get sums across 6 groups at each time, construct the matrix that sums groups at each time - you should have 6 columns. And to get the total across time for all groups use a 24x1 matrix, although that makes no sense to me unless you are going to divide by 6 to get a mean.

I highly recommend that you read the material on the delta method in Cooch and White so you understand what is going on. Although in this case the application is simply sums and as such the variance of a sum is the sum of vcv matrix of the Ns that you are summing. You may have heard that the variance of a sum is the sum of the variance + 2x the sum of the covariances but since the covariances are in the upper and lower diagonal of the matrix, if you some the entire matrix it is the same thing. So what you can do is to get the 24 individual estimates and the 24 x 24 vcv matrix of the N's with a single call of pop.est and do the summations yourself because that is all the design matrix is doing when you give a 24x1 or 24x6 matrix.

Also, I believe the paper(s) should have said that a log normal confidence interval was used to construct the confidence intervals. Plenty of references for that (ie distance sampling, AFS Monograph #5) and I even think you may find it on wiki.

--jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA

Re: CJS abundance - pop.est from RMark

Postby SoConfused » Mon Apr 02, 2018 9:36 am

Ok, I think I got it - here's a reworked dipper example (same setup as in the question). The estimates and variances sum up correctly now; I looked up the distance 95% CIs and I think this is how they would work.

Is this correct? Also - can this be used on averaged models? Or only on a single model at a time?

# data and 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()

md <- data1.results[[1]]


# pop.est estimates for each sampling event (across groups)
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))
original.total <- pop.est(ns, p.list$estimates$estimate, design=diag(1,ncol=6,nrow=6),p.list$vcv)


# now by group
Code: Select all
n.occ <- nchar(dipper$ch[1]) - 1
n.group <- length(unique(dipper$sex)) * length(unique(dipper$Var1))

ns1 <- colSums(xmat[dipper$sex == "Female" & dipper$Var1 == "A",])[-1]
ns2 <- colSums(xmat[dipper$sex == "Male" & dipper$Var1 == "A",])[-1]
ns3 <- colSums(xmat[dipper$sex == "Female" & dipper$Var1 == "B",])[-1]
ns4 <- colSums(xmat[dipper$sex == "Male" & dipper$Var1 == "B",])[-1]
ns <- c(ns1, ns2, ns3, ns4)

p.indices <- extract.indices(md,"p",df=data.frame(group=rep(1:n.group, each = n.occ), row=rep(1:n.occ, n.group), col=rep(1:n.occ, n.group)))
p.list <- covariate.predictions(md,data=data.frame(index=p.indices))
by.group <- pop.est(ns, p.list$estimates$estimate, design=diag(1,ncol=n.occ * n.group,nrow=n.occ * n.group),p.list$vcv)


# to sum up by occasion
Code: Select all
out <- list()
vars <- list()

for(i in 1:n.occ){
   out[[i]] <- sum(by.group$Nhat[c(i, i + (1:(n.group-1)) * n.occ)])
   vars[[i]] <- sum(by.group$vcv[c(i, i + (1:(n.group-1)) * n.occ), c(i, i + (1:(n.group-1)) * n.occ)])
        }   

vars <- unlist(vars)
out <- unlist(out)
vars
out
original.total ### the total original estimates is the same as the summed up values

# 95% CIs, from http://distancesampling.org/downloads/dist_encyc_env.pdf
VarLn <- log(1 + vars / out^2)
C <- exp(1.96 * sqrt(VarLn))

# 95% CI is D/C, D*C
out/C
out * C
SoConfused
 
Posts: 53
Joined: Wed Nov 05, 2014 8:25 am

Re: CJS abundance - pop.est from RMark

Postby jlaake » Mon Apr 02, 2018 9:46 am

I didn't run your code but it all looks correct to me. You might want to compute one by hand as a check on your code. If you want to do over multiple models, I would use covariate.predictions with a model.list to model average the p's and then use those values and their vcv matrix where you used the values for a single model. Obviously the n's don't change over models.

--jeff
jlaake
 
Posts: 1417
Joined: Fri May 12, 2006 12:50 pm
Location: Escondido, CA


Return to RMark

Who is online

Users browsing this forum: No registered users and 14 guests