chat adjustment on derived parameter N - Huggins

[monicaarso]

Hi,

I'm running RDHFHet models accounting for trap dependence given GOF tests with R2ucare. Additionally, I am using adjust.chat as I still get a chat~1.6 after recalculating withouth the test2ct component (some may argue that is low enough to ignore, but for the sake of the question in place, let's assume I want to adjust).

What I have not managed to do is to use get.real for the derived parameter N so that the 95%CI around the estimates are inflated given the chat.

Is there an equivalent to popan.derived I can use (to which I'll have to add covariate data given the models)?

Many thanks,
Monica

[jlaake]

Use the chat argument in mark function to run a specific model in MARK or as argument to mark.wrapper to run a set of models. This will then adjust all std errors and confidence intervals for everything including derived parameters. It will also use quick with chat adjustment.

[monicaarso]

Hi Jeff,
Thanks for the prompt reply.
I reran the best model but with chat=1.56 in the mark call, but the CI around the N estimates are still the same when I call for the derived parameter:

Code: Select all

#Best model but with chat in the call mov3.chat=mark(Tay.process,Tay.ddl,model="RDHFHet",threads=1,model.parameters=list(S=S.dot,p=p.time.session.mixture.td,GammaDoublePrime=GammaDoublePrime.markov.constant,GammaPrime= GammaPrime.markov.constant,pi=pi.dot),chat=1.561125)
#get N estimates from model mov3
mov3$results$derived$`N Population Size`
#get N estimates from model mov3.chat
mov3.chat$results$derived$`N Population Size`

Cheers,
Monica

[jlaake]

My error. The beta and real estimate std error and intervals should have changed but apparently I forgot to do the derived parameters. I was thinking that adding chat to call would pass it to MARK but that only happens with profile intervals. The chat adjustments are done in the MARK interface rather than in mark.exe. Thus it has to be done in RMark. This is straightforward but the intervals will depend on the type of parameter so this will take some work on my end before it does it for all derived parameters. In the meantime, many of the derived parameters are abundance estimates so I wrote the following code that you can use for that instance which uses the log-normal distribution to construct the interval. I show it with a POPAN example and chat=1.5 but it will work for any model where all of the derived values use a log-normal confidence interval.

Code: Select all

adjust.derived=function(x,chat)
{
  C=exp(1.96*sqrt(log(1+chat*(x$se/x$estimate)^2)))
  return(data.frame(estimate=x$estimate,se=x$se*sqrt(chat),lcl=x$estimate/C,ucl=x$estimate*C))
}

data(dipper)
model=mark(dipper,model="POPAN")

model$results$derived

lapply(model$results$derived,adjust.derived,chat=1.5)

The code above just reports the new values. If you want to store them back into the model you can do that with

Code: Select all

model$results$derived=lapply(model$results$derived,adjust.derived,chat=1.5)

Note that if you store them back in the model object, you can't then adjust them again with a different chat so you may want to store them separately or with a different name (e.g., model$results$derived.chat).

Also, if you want to adjust the var-cov matrix of the derived parameters use

Code: Select all

lapply(model$results$derived.vcv,function(x) x*1.5)


--jeff

[monicaarso]

Hi Jeff,
Many thanks for the provided code.
A quick follow up question: when using model.average function in a CJS model for real parameter model average (e.g. survival) I take it the model.average will not automatically inflate 95% CI in the presence of chat>1 right? So I would have to manually adjust the averaged output if I wanted to?
Monica

[jlaake]

Monica

Sorry for the slow reply. I wanted to take time to create an example for you. The short answer is that if you use chat in the call to mark.wrapper for the models then the model averaged estimates will incorporate the chat value. Now what is important to realize is that chat adjusts model weights. Run the following example with the dipper data and you'll see that the results can be counterintuitive due to the change in weights. I use model.average and also show how the estimate and std error are computed for this simple example. Let me know if you have any questions.

Code: Select all

library(RMark)
data(dipper)
dp=process.data(dipper)
ddl=make.design.data(dp)
fit_models=function(chat=1)
{
  Phi.dot=list(formula=~1)
  Phi.time=list(formula=~time)
  p.dot=list(formula=~1)
  cml=create.model.list("CJS")
  results=mark.wrapper(cml,data=dp,ddl=ddl,chat=chat)
  return(results)
}

# fit models and look at the results - note that the weights have changed
# models with chat=1
results.1=fit_models()
results.1
# models with chat=2
results.2=fit_models(chat=2)
results.2

# look at Phi values and std error for constant model - shows that chat is incorporated
# survival for time 1 with constant survival model and chat=1 and chat=2
summary(results.1[[1]],se=TRUE)$reals$Phi[1,]
summary(results.2[[1]],se=TRUE)$reals$Phi[1,]
# .025133*sqrt(2)=0.03554343
# survival for time 1 with time-varying survival model and chat=1 and chat=2
summary(results.1[[2]],se=TRUE)$reals$Phi[1,]
summary(results.2[[2]],se=TRUE)$reals$Phi[1,]
# 0.1116453*sqrt(2)=0.1578903

# now look at model averaged results
model.average(results.1,parameter="Phi")[1,]
model.average(results.2,parameter="Phi")[1,]

# let's compute this for results.1 with chat=1
# first for model 1
m1.estimate=summary(results.1[[1]],se=TRUE)$reals$Phi[1,"estimate"]
m1.var=summary(results.1[[1]],se=TRUE)$reals$Phi[1,"se"]^2
# next for model 2
m2.estimate=summary(results.1[[2]],se=TRUE)$reals$Phi[1,"estimate"]
m2.var=summary(results.1[[2]],se=TRUE)$reals$Phi[1,"se"]^2
# the model averaged esimate for time 1
mavg.estimate=model.average(results.1,parameter="Phi")[1,"estimate"]
# the model averaged std error using revised=TRUE
mavg.se=sqrt(results.1$model.table$weight[1]*(m1.var+(m1.estimate-mavg.estimate)^2)+
               results.1$model.table$weight[2]*(m2.var+(m2.estimate-mavg.estimate)^2))
mavg.se
model.average(results.1,parameter="Phi")[1,"se"]

# now let's do the same thing for chat=2 (results.2)
# first for model 1
m1.estimate=summary(results.2[[1]],se=TRUE)$reals$Phi[1,"estimate"]
m1.var=summary(results.2[[1]],se=TRUE)$reals$Phi[1,"se"]^2
# next for model 2
m2.estimate=summary(results.2[[2]],se=TRUE)$reals$Phi[1,"estimate"]
m2.var=summary(results.2[[2]],se=TRUE)$reals$Phi[1,"se"]^2
# the model averaged esimate for time 1
mavg.estimate=model.average(results.2,parameter="Phi")[1,"estimate"]
# the model averaged std error using revised=TRUE
mavg.se=sqrt(results.2$model.table$weight[1]*(m1.var+(m1.estimate-mavg.estimate)^2)+
               results.2$model.table$weight[2]*(m2.var+(m2.estimate-mavg.estimate)^2))
mavg.se
model.average(results.2,parameter="Phi")[1,"se"]

# These results may seem counterintuitive because with chat=2,
# the std error is smaller than when chat=1. The reason is due to the re-weighting
# of the models using QAIC vs AIC.  With chat=2 the model with time dependent survival has less weight (0.03)
# than when chat=1 (0.17).  The survival estimate for time 1 with time dependent survival is
# larger (0.625 vs 0.56) and less certain (se 0.111 and 0.157 for chat 1 and 2) than the
# estimate for the constant survival model (0.025 and 0.035 for chat=1 and 2).


[monicaarso]

Jeff, once again, thanks very much for the post and the provided code. I had not thought about the effect of the weights on the model.average so this is a great example.
Many thanks,
Monica