Summing abundance estimates across groups

[jlaake]

I have been working with someone offlist that is using the RDHHet model and they have groups in the data. I noticed that MARK produces the abundance by group for each session in the derived parameters but does not provide a total abundance for each session (Gary if I have this wrong please let me know. May need to modify things in RMark). I wrote the code below that produces the estimates, std errors and confidence intervals for the total (sum) abundance for each session. Posting here as it may be useful for others and it also shows how log-normal confidence intervals are constructed for abundance for these types of models. Also, I show an example of how to model average as I have gotten questions like this before for derived parameters. Don't try to run the example for model averaging because it will not work. But you can use the functions and example with your models. Let me know if you have any problems.

Code: Select all

# compute coefficient of variation and log-normal confidence interval on f0
# and then add MT1 to end points of confidence interval because they are constants
# returns data.frame with lower and upper confidence interval values
log_normal_ci_f0=function(f0,se,MT1)
{
  cv=se/f0
  C=exp(1.96*sqrt(log(1+cv^2)))
  LCL=f0/C+MT1
  UCL=f0*C+MT1
  return(data.frame(LCL=LCL,UCL=UCL))
}

# sum abundance estimates over groups for each session and return data.frame with estimates,std errors, and log-normal
# confidence intervals
combine_group_abundances=function(model)
{
  # get number of groups and number of occasions
  ngrps=model$number.of.groups
  nocc=model$nocc
  # get abundance estimates from derived parameters results
  abundance_estimates=cbind(grp=rep(1:ngrps,each=nocc),session=rep(1:nocc,times=ngrps),model$results$derived[[1]])
  # get total across groups for each session 
  total_abundance_by_session=cbind(session=abundance_estimates[1:nocc,2],total=sapply(split(abundance_estimates$estimate,abundance_estimates$session),sum))
  # Compute variance-covariance matrix of total for each session
  derivatives=matrix(0,nrow=nocc,ncol=ngrps*nocc)
  for(i in 1:nrow(derivatives))
    derivatives[i,seq(i,ngrps*nocc,nocc)]=1
  vc=derivatives%*%model$results$derived.vcv[[1]]%*%t(derivatives)
  # sqrt of diagonal is the std error of the total
  stderrors=sqrt(diag(vc))
  #now getting the confidence interval is a little more involved because the interval is set using
  # a log-normal distribution on f0 where f0=N-Mt1 where N is the abundance estimate and Mt1 is the total
  # number seen. This makes sure that the lower bound on the confidence interval is at least the number you saw in the session
  # get Mt+1 for each group by primary occasion
  x=model$data$data
  sec=model$nocc.secondary
  # get number seen by group and session
  get_mt1=function()
  {
    Mt1=matrix(0,nrow=length(sec),ncol=length(levels(x$group)))
    start= cumsum(c(1,sec))[-(length(sec)+1)]
    stop=cumsum(c(0,sec))[-1]
    for(i in 1:length(sec))
      Mt1[i,]=sapply(split(as.numeric(colSums(sapply(strsplit(x$ch,""),function(x)as.numeric(x[start[i]:stop[i]])))>0)*x$freq,x$group),sum)
    return(Mt1)
  }
  # sum across groups
  MT1=rowSums(get_mt1())
  # compute number not seen which is abundance estimate minus number seen
  f0_by_session=total_abundance_by_session[,"total"]-MT1
  Total_abundance=data.frame(total_abundance_by_session,se=stderrors,log_normal_ci_f0(f0_by_session,stderrors,MT1))
  return(list(MT1=MT1,Total_abundance=Total_abundance,vcv=vc))
}


Here is an example to show how you would compute model averaged values

Code: Select all

# Now you will most likely want to model average across a set of models. You will want to run combine_group_abundances
# for each model and construct a list containing estimate matrix, AIC measure or weight vector, and se matrix
# and then call model.average for list format

# demo with first 3 models
# get model numbers of first 3 models (if you want to do for all models, just replace 3 with the number of models)
number_of_models=3
model_num=as.numeric(rownames(rd.noemigration.results$model.table))[1:number_of_models]
# get AICc values
AICc=rd.noemigration.results$model.table$AICc[1:number_of_models]
# construct empty matrix for estimates and std errors
estimate=matrix(NA,nrow=number_of_models,ncol=rd.process$nocc)
se=matrix(NA,nrow=number_of_models,ncol=rd.process$nocc)
for(i in 1:number_of_models)
{
  xlist=combine_group_abundances(rd.noemigration.results[[i]])
  if(i==1)MT1=xlist$MT1
  estimate[i,]=xlist$Total_abundance$total
  se[i,]=xlist$Total_abundance$se
}
# call model.average for a list
mavgs=model.average(x=list(estimate=estimate,AICc=AICc,se=se))

# now compute confidence intervals for model average estimates and combine all into a data.frame
mavg_df=data.frame(estimate=mavgs$estimate,se=mavgs$se,log_normal_ci_f0(mavgs$estimate-MT1,mavgs$se,MT1))

mavg_df


[jlaake]

Thanks to Gary for adding the code to MARK for these models. He will release soon. Note that the total abundance across groups for each session (All) will only appear in the output and not in the derived parameters, so it will not appear in the RMark derived object. So you can use this code or just snip results from the MARK output. If for some reason you want to use the v-c matrix of the All values, you will have to use the code posted here.

--jeff

[gwhite]

I have added code to the robust design and multi-state robust design models to compute the sum of the N estimates across groups, both Huggins and full likelihood data types. This version is on my web site, and I'm sure Evan will have it copied to phidot soon.

Gary

[cooch]

I have added code to the robust design and multi-state robust design models to compute the sum of the N estimates across groups, both Huggins and full likelihood data types. This version is on my web site, and I'm sure Evan will have it copied to phidot soon.

Gary

So done, as of 30 seconds ago.