[R] popbio and stochastic lambda calculation

[Chris Stubben]

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




-- 
View this message in context: http://n4.nabble.com/popbio-and-stochastic-lambda-calculation-tp1557647p1560745.html
Sent from the R help mailing list archive at Nabble.com.