That is a big help Chris - thank you very much for your assistance. I
also ran your code (below) for the example in Caswell (2001, ~p300) and
got matching results.
I'd like to add my vote for including this method in the next version of
popbio - I think it would be a very useful addition.
Thanks again,
Shawn
Chris Stubben wrote:
Shawn Morrison-2 wrote:
# The paper reports a 95% CI of 0.79 - 1.10
# "My" reproduced result for the CIs is much larger, especially on the
upper end. Why would this be?
# The authors report using the 'delta' method (Caswell, 2001) to
calculate the CI in which the
Shawn,
I probably can't help much with the vitalsim example, but I would check Box
8.10 in Morris and Doak (2002).
I do have a few ideas about the delta method below.
# List the vital rates
vr<-list(cub= 0.64, yly = 0.67, sub=0.765, adt=0.835, mx=0.467)
# and the matrix using an expression in R
el<- expression(
0, 0, 0, 0, adt*mx,
cub,0, 0, 0, 0,
0, yly,0, 0, 0,
0, 0, sub,0, 0,
0, 0, 0, sub,adt)
# this should get the projection matrix
A<-matrix( sapply( el, eval, vr), nrow=5, byrow=TRUE)
lambda(A)
[1] 0.9534346
# use the vitalsens function to get the vital rate sensitivites and
save the second column
vitalsens(el, vr)
estimate sensitivity elasticity
cub 0.640 0.1236186 0.08298088
yly 0.670 0.1180835 0.08298088
sub 0.765 0.2068390 0.16596176
adt 0.835 0.7628261 0.66807647
mx 0.467 0.1694131 0.08298088
sens<-vitalsens(el, vr)[,2]
# I'm not sure about the covariance matrix next. In Step 7 in Slakski et al
2007 ("Calculating the variance of the finite rate of population change from
a matrix model in Mathematica") they just use the square of the standard
errors, so I'll do the same...
se<-list(cub= 0.107, yly = 0.142, sub=0.149, adult=0.106, mx=0.09)
cov<-diag(unlist(se)^2)
## and then the variance of lambda from step 8
var<-t(sens) %*% ( cov%*%sens)
[,1]
[1,] 0.008176676
# and the confidence intervals
lambda(A) - 1.96*sqrt(var)
lambda(A) + 1.96*sqrt(var)
CI of 0.78 and 1.13 is close to paper
Hope that helps,
Chris
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
