RMARK:Model Averaged Confidence Intervals for Robust Design

[cb]

I am using RMark to run Robust Design (RDHuggins) models. I used this great forum to figure out how to model average the N estimates (derived parameters), but I can't figure how to get it to calculate the confidence intervals.

Can someone help me?

Thanks!

Here is the code we are using (basically taken from forum):

#Average models - derived parameters - Nhat
derNavg = function(fil)
{
num.models=nrow(fil$model.table) #stores num models for loop
num.years=length(fil[[1]]$results$derived$estimate)
estimate=matrix(0,ncol=num.years, nrow=num.models) #stores ests
se=matrix(0,ncol=num.years, nrow=num.models) #stores se
vcv = vector("list",length = num.models) #stores vcv
weight=fil$model.table$weight
#create list containing estimate and se for derived param for each model
for(i in 1:num.models)
{
# The actual model number is the row number for the model.table
model.numbers= as.numeric(row.names(fil$model.table))
# For each model extract the derived parameter valuess and their se
x=fil[[model.numbers[i]]]$results
estimate[i,]=x$derived$estimate
vcv[[i]]= x$derived.vcv
se[i,]=x$derived$se
}

#Call model.average using the list which includes est, wt, vcv list
# using vcv didn't work for 37213...don't know why, using se works
# vcv worked for all other meadows
# return(model.average(list(estimate=estimate,weight=weight,vcv=vcv)))
#didn't work, vcv wants list..not same as marklist
#return(model.average(list(estimate=estimate,weight=weight,se=se, vcv=TRUE)))

return(model.average(list(estimate=estimate,weight=weight,se=se)))
}

#call
m15Navg = derNavg(m15malesrobhSess.results_v2)

[cooch]

I am using RMark to run Robust Design (RDHuggins) models. I used this great forum to figure out how to model average the N estimates (derived parameters), but I can't figure how to get it to calculate the confidence intervals.

I'll leave the specifics to Jeff or one of the other 'usual suspects', but will remind you that calculating the SE or CI for modeled averaged N requires some extra work (relative to doing same for other parameters). This is described in detail in chapter 14 (chapter on abundance estimation). See especially section 14.9.1.

[jlaake]

It is relatively simple to compute the CI with 3 equations (see 14.9) from the results you have. If the estimates are in vector est and the std errors in vector se, then the following should work.

C=exp(1.96*sqrt(log(1+(se/est) ^2)))
lcl=est/C
ucl=est*C

Now your est is of N so if you want to base it on f0 as shown in 14.9, you'll need a vector of Mt+1 (number caught). I'll call that Mt1. Then you would do

f0=est-Mt1
C=exp(1.96*sqrt(log(1+(se/f0) ^2)))
lcl=Mt1+f0/C
ucl=Mt1+f0*C

Unless I made some typos the above should work with the appropriate values assigned to the vectors. I've not used RD Huggins so not sure what is being derived here and how you compute Mt1 in that case. This may seem strange but it is quite easy to program these models into RMark without really understanding what they do because MARK does all the work. So there are many models in RMark that I'm not familiar with because I've never used them other than with a dummy example.

--jeff

[cb]

Ok, This makes sense, and the values I get seem right. The code works fine.
To make sure I am thinking about this correctly:

For the Robust Design, since I am interested in the yearly N's, I use the closed capture procedures because this corresponds to the secondary closed periods. I calculate the CI's for each year using that year's est of N and se (as you describe). The Mt+1 would be for each year, and these are listed in the output file (and I know them).

One clarification on your instructions (and described in 14.9). I used the model averaged N for the est vector, and the model averaged se for the se vector. Is this correct? CI's seem reasonable and weighted.

Finally, I assume the fo routine is the version I want to use and not the profile confidence interval method which according to 14.9 - is not the best for closed capture.

Am I correct in this?

Thanks for your help!!

cb

[jlaake]

As I understand it, everything you said is correct. --jeff

[cb]

Very Good.
Thanks once again for your help!