[R] popbio and stochastic lambda calculation

[Shawn Morrison]

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
>
>
>
>
>