RMark model averaged f(0) < 1.0

[howeer]

Hi,
Whenever I try to get model averaged f(0) estimates (in order to get model averaged N by adding f(0) to M(t+1)), the f(0) estimates are < 1. Jeff has heard about this before, and identified the problem as a problem with the link function for N, but even running newest versions of R and RMark I always get f(0) < 1.0. Upper confidence limits for f(0) were > 1.0, however.

Are other folks able to get reasonable estimates of model averaged f(0)from RMark? Is there still a bug, or am I missing something?

Thanks

[jlaake]

Eric-

Below is an example using a modified version of the function contained in the help for the edwards.eberhardt example data set. As you can see the model averaged result is 19.078 for f(0) which when added to M(t+1) gives an average N of 95 which is the result for the third model which has almost all of the weight. Have you checked to see if the estimated N is greater than M(t+1) for your models. It is quote common for f(0) <1 if the model thinks all animals have been caught over the course of the study (even if that is not the case). There is a recent paper by Jason Baker on monk seals that examines the problems with heterogeneity with closed capture models. It was quite common that his estimate of N was M(t+1) even though he had observations of missed seals that were seen outside the sampling occasions.

As far as I know the code is working correctly. I hadn't heard from you in awhile and assumed you had worked out any issues. I used the current version for the example which is v1.7.7 posted 6 March 2008.

regards --jeff

> run.edwards.eberhardt
function()
{
#
# Define parameter models
#
pdotshared=list(formula=~1,share=TRUE)
ptimeshared=list(formula=~time,share=TRUE)
#
# Run assortment of models
#
#
# Capture Closed models
#
# constant p=c
ee.closed.m0 =mark(edwards.eberhardt,model="Closed",model.parameters=list(p=pdotshared))
# constant p and constant c but different
ee.closed.m0c =mark(edwards.eberhardt,model="Closed")
# time varying p=c
ee.closed.mt =mark(edwards.eberhardt,model="Closed",model.parameters=list(p=ptimeshared))
#
# Return model table and list of models
#
return(collect.models() )
}
> ee.results=run.edwards.eberhardt()

> ee.results
model npar AICc DeltaAICc weight Deviance
3 p(~time)c()N(~1) 19 354.5968 0.00000 9.999949e-01 305.2648
1 p(~1)c()N(~1) 2 379.6029 25.00616 3.715170e-06 364.8259
2 p(~1)c(~1)N(~1) 3 381.5498 26.95301 1.403543e-06 364.7640
> summary(ee.results[[1]])$real$N
[,1]
96.25877
> summary(ee.results[[2]])$real$N
[,1]
93.89185
> summary(ee.results[[3]])$real$N
[,1]
95.07817

> model.average(ee.results,"N")
par.index estimate se fixed group
N g1 36 19.07818 6.613443 1
> 19.078+76
[1] 95.078

[cooch]

> model.average(ee.results,"N")
par.index estimate se fixed group
N g1 36 19.07818 6.613443 1
> 19.078+76
[1] 95.078

Just out of curiosity - how does model.average handle the unconditional CI calculation? The generic approach used in MARK is incorrect for model averaged abundance estimates (as described in sec. 15.8.1). I know I could look at the RMark source, but I'm too lazy.

[jlaake]

Following on with my example from the last posting, the code and results below show how the ci is computed for f(0) with RMark and then how it is done with a log-normal as shown in 15.8.1 in the book.

In RMark, the unconditional se for the real parameters is converted to the link space (log in this case) using the delta method, a 95% ci is computed in the link space, and then the end-points are computed via back-transformation. For f(0) this is close to the log-normal but not the exactly the same.

> model.average(ee.results,"N",vcv=T)
$estimates
par.index estimate se lcl ucl fixed group
N g1 36 19.07818 6.613443 9.67081 37.63664 1

$vcv.real
36
36 43.73763

# confidence interval in RMark
# computes se of log(f0), normal ci on log(f0), back-transformation of endpoints
> exp(log(19.07818)-1.96*(6.613443/19.07818))
[1] 9.67081
> exp(log(19.07818)+1.96*(6.613443/19.07818))
[1] 37.63665

# log-normal confidence interval
> exp(1.96*sqrt(log(1+(6.61344/19.078)^2))) # computation of C
[1] 1.935291
# lower limit - divide estimate by C
> 19.078/1.935
[1] 9.859432
# upper limit - multiply estimate by C
> 19.078*1.935
[1] 36.91593

[howeer]

Hi and thanks,
I hadn't posted because in the analysis I had been doing, the top-ranked model had 98% of the total weight, so I didn't model average.

Yes, with the current analysis, N-hat from individual models was considerably > 1, and 1.7.7 is the version I'm using. Results from individual models are reasonable, but inconsistent with model averaged estimates. Maybe it's a problem in my sript; I sent examples to Jeff.
Thanks again,
Eric

[jlaake]

Eric's assessment was correct that there was a bug in RMark. I was using the Closed model in my example and it was setup ok. However, for the FullHet and HetClosed models, the default link for N was being incorrectly set to logit which makes no sense for N because it restricts f(0) <1 as he found. An error caused from copying code and I should know better. Now MARK correctly ignores any setting for the link with N and uses the log link which is why his MARK runs looked fine. However, when he used model.average it incorrectly used the logit link with the beta for N which restricted his f(0) <1. I checked all of the settings for N in the function setup.parameters and they are all now set to log.

On another point in regard to my previous posting, the reason for the difference in the CI for N between RMark and the log-normal is in the approximation of the variance with the delta method. The delta method uses cv^2 for the variance and the log-normal variance is log(1+cv^2). The delta method variance is slightly larger but as long as the cv<0.5 the difference doesn't exceed 5% for the se (see below).
cv sqrt(ln(1+cv^2)) ratio
0.05 0.05 1.00
0.10 0.10 1.00
0.15 0.15 1.01
0.20 0.20 1.01
0.25 0.25 1.02
0.30 0.29 1.02
0.35 0.34 1.03
0.40 0.39 1.04
0.45 0.43 1.05
0.50 0.47 1.06
0.55 0.51 1.07
0.60 0.55 1.08
0.65 0.59 1.09
0.70 0.63 1.11

[jlaake]

Forgot to say that the bug has been fixed in v1.7.8 which has been posted to phidot.

[cooch]

On another point in regard to my previous posting, the reason for the difference in the CI for N between RMark and the log-normal is in the approximation of the variance with the delta method. The delta method uses cv^2 for the variance and the log-normal variance is log(1+cv^2). The delta method variance is slightly larger but as long as the cv<0.5 the difference doesn't exceed 5% for the se

Which, in effect, is due to the non-linearity of the log transform. The standard application of the 'Delta method' involves a simple first order Taylor expansion, and thus is a bit 'twitchy' when the transform is non-linear, and if the variation in the data is 'too big'. This is discussed in a bit of detail in Appendix B (see section B.3.1). In such cases, a second-order Taylor expansion can help, but it doesn't always solve the problem.