[R] popbio and stochastic lambda calculation
Shawn Morrison
shawn.morrison at dryasresearch.com
Fri Feb 19 17:32:07 CET 2010
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
>
>
>
>
>
More information about the R-help
mailing list