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