[R] popbio and stochastic lambda calculation

[Shawn Morrison]

Hello R users,

I am trying to calculate the stochastic lambda for a published matrix 
population model using the popbio package.

Unfortunately, I have been unable to match the published results. Can 
anyone tell me whether this is due to slightly different methods being 
used, or have I gone wrong somewhere in my code?

Could the answer be as simple as comparing deterministic lambdas to 
stochastic lambdas? My stoch. lambda is lower (as expected) however the 
CI is substantially larger than published. The authors of the published 
example used a different method to calculate the CIs around lambda, but 
I would have expected the results of the vitalsim to be closer. (The 
authors used the 'delta' method from Caswell (2001)).

My main questions are:
1. Have I correctly calculated stochastic lambda and its 95% CIs?
2. Why am I unable to use any distributions other than beta? What can I 
do about that?

Thank you,
Shawn Morrison

#####################################################

#My example:

rm(list = ls())
objects()

library(popbio)


# Vital rate means and variances, and 'types' for the vrtypes argument 
in vitalsim
# 'names' is not used, but indicates what the vital
# rates represent: Sad = adult survival, Scub = cub survival
# Syrl = yearling survival, Ssub - subadult survival
# mx = number of female offspring per year

names = c("Sad", "Scub", "Syrl", "Ssub", "mx") # vital rate names, not used
mean = c(0.835, 0.640, 0.670, 0.765, 0.467) # vital rate means
var = c(0.213, 0.252, 0.241, 0.133, 0.0405) #variances of means
var.se = c(0.106, 0.107, 0.142, 0.149, 0.09) #standard errors of means
types = c(1,1,1,1,1) #for vrtypes argument


ex.vrs = data.frame(cbind(mean, var, var.se, types))
attach(ex.vrs)

## Define the matrix
ex.mxdef= function (vrs)
{
    matrix(c(0,            0,          0,              0,           
vrs[5]*vrs[1],
            vrs[2],          0,          0,              
0,                       0,
             0,         vrs[3],          0,              
0,                       0,
             0,            0,          vrs[4],          
0,                        0,
             0,            0,           0,             
vrs[4],             vrs[1]),
             nrow = 5, byrow = TRUE)
}



ex.mxdef(ex.vrs$mean)

# Matrix in published article is:
# 0,       0,   0,      0,  0.390,
# 0.640,   0,   0,      0,    0,
# 0,      0.67, 0,      0,    0, # 0,        0,  0.765,  0,    0,
# 0,        0,  0,      0.765, 0.835
# "My' result matches this.

###  TRIAL 1
###  Using variances estimated from published paper returns an error 
(ex.var$var), perhaps because variances are too large?

## run no.corr model
no.corr<-vitalsim(ex.vrs$mean, ex.vrs$var, diag(5),
diag(5), ex.mxdef, vrtypes=ex.vrs$types,
n0=c(200,130,90,80,490), yrspan=20 , runs=200)

###  TRIAL 2
###  Using standard error instead of variances (ex.var$var.se)
### This also produces an error when 'types' is set to anything but 1 
(betas).

## run no.corr model
no.corr<-vitalsim(ex.vrs$mean, ex.vrs$var.se, diag(5),
diag(5), ex.mxdef, vrtypes=ex.vrs$types,
n0=c(200,130,90,80,490), yrspan=20 , runs=200)

## deterministic lambda
no.corr[1:2]  # The published deterministic lambda is 0.953 so this 
appear to be working correctly


?popbio

# stochastic lambda
no.corr[2]
# 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
# sensitivities were used to weight the variances.

## log stochastic lambda
log(no.corr$stochLambda)
sd(no.corr$logLambdas)

se=function(x) sqrt(var(x)/length(x))
se(no.corr$logLambdas)  # The published article reports deterministic 
lambda ± SE to be 0.953 ± 0.079