www.phidot.org

Hello,
I am currently trying to apply delta method following Mark book procedure presented in appendix B, using the last version of ESURGE. To test the validity of my method, I worked on classical European Dipper dataset and run the model (StP.). Survival estimates with ESURGE are quite similar to estimates found with Mark and presented in Mark book. Nevertheless, var-cov matrix is quite different. I applied delta method on three first survival estimates (Y=Phi1*Phi2*Phi3) using var-cov matrix presented in the result in ESURGE and found
ESURGE (generalized logit link function): Y=0.13598±0.14453 (SE)
ESURGE (Identity link function): Y=0.13598±0.03488 (SE)
MARK (from data presented in Mark book): Y=0.13887127±0.05106 (SE)
Note that dataset used with Mark program is likely to differ slightly from the one used with ESURGE.
Nevertheless, it seems that valid results (close to the one obtained with Mark) are obtained with ESURGE using identity link function only. Why?
Thanks for your help
Cheers
Jocelyn

Did you take into account the fact that the var-cov matrix is
for the beta and not for the survival in the case of the logit?
Rémi

Fine! thanks for that.
But in case of logit, applying delta method is very hard!
- For instance, taking the example presented before: Y=Phi1*Phi2*Phi3. If I take var-cov matrix with Beta values. So I transform in the case of dipper data Phi=exp(B)/(1+exp(B)). So, I need to apply delta method to estimate var(Y) from derivatives of Y=[exp(B1)/(1+exp(B1)]*[exp(B2)/(1+exp(B2)]*[exp(B3)/(1+exp(B3)], which is quite complicated?
Sorry for what is perhaps simple question ...
Many thanks
Jocelyn

Jocelyn wrote:Fine! thanks for that.
But in case of logit, applying delta method is very hard!
- For instance, taking the example presented before: Y=Phi1*Phi2*Phi3. If I take var-cov matrix with Beta values. So I transform in the case of dipper data Phi=exp(B)/(1+exp(B)). So, I need to apply delta method to estimate var(Y) from derivatives of Y=[exp(B1)/(1+exp(B1)]*[exp(B2)/(1+exp(B2)]*[exp(B3)/(1+exp(B3)], which is quite complicated?
Sorry for what is perhaps simple question ...
Many thanks
Jocelyn

MARK lets you output the V-C matrix for either the betas or the reals, which avoids the problem you describe (which actually isn't too big a problem -- the derivative for the logit is trivial, and can even be done in R. If you want a more 'symbolic' approach to handling the math, l'd have a look at Maxima).

See also the delta method using R:

http://rss.acs.unt.edu/Rdoc/library/msm ... ethod.html

Sincerely,

Rémi

The emdbook package in Program R has a function deltavar() that will use the delta method to calculate the variance of any function (Bolker 2008). Arguments needed include a function, a list of the parameters, and estimates of the mean value and SD for each input parameter. The argument Sigma can be a vector of variances if the parameters are independent, or a var-covar matrix if parameters are correlated.

# need to install package emdbook
library(emdbook)
# calculate variance of lambda
# lambda = (fecundity * juvenile survival) + adult survival
par.name = c("F", "Sj", "Sa")
par.mean = c(F=3.4, Sj=0.25, Sa=0.60)
par.SD = c(0.5, 0.02, 0.05)
par.var = par.SD^2
deltavar(F*Sj+Sa, meanval=par.mean, vars=par.name, Sigma=par.var)
# check answer, var(lambda) = 0.022749

B.K. Sandercock wrote:The emdbook package in Program R has a function deltavar() that will use the delta method to calculate the variance of any function (Bolker 2008). Arguments needed include a function, a list of the parameters, and estimates of the mean value and SD for each input parameter. The argument Sigma can be a vector of variances if the parameters are independent, or a var-covar matrix if parameters are correlated.

# need to install package emdbook
library(emdbook)
# calculate variance of lambda
# lambda = (fecundity * juvenile survival) + adult survival
par.name = c("F", "Sj", "Sa")
par.mean = c(F=3.4, Sj=0.25, Sa=0.60)
par.SD = c(0.5, 0.02, 0.05)
par.var = par.SD^2
deltavar(F*Sj+Sa, meanval=par.mean, vars=par.name, Sigma=par.var)
# check answer, var(lambda) = 0.022749

Bolkers method works pretty well, but I'd generally recommend getting comfortable with a CAS (computer algebra system) for more complicated transforms. R is very powerful, but a 'symbolic math environment' it isn't (it can do rudimentary calculus, but not much more). Maxima is a very good (and free - which of course is also a major selling point for R) CAS (Maple, and Mathematica are *very, very* good, but also *very, very* expensive). And don't forget, the Delta approximation is just that - it works pretty well *in some cases*, but not others (see Appendix B of the MARK book for details).

As an aside, lambda = Sa + F*Sj is a valid population model, but you need to be careful. Frequently, this model is attributed to Pulliam (1988) - his Am Nat paper on 'source/sink' dynamics. Specifically, his eqn. (5a) and (5b), which relate to his figure 1. Here is an expanded explanation of where the equation for lambda comes from. Consider a pre-breeding census for a population with 2 possible age classes (adults and babies). Prior to breeding, there is only a single age class (adults). Thus is the number of adults at time t. Consider the following expression projecting to :



In the first term on the RHS, we see that adults this year (time t) survive to be adults next year (time t+1), with probability . The second term is the fertility contribution: number of adults this year tims the per adult (per capita) fertility rate (i.e., total number of babies produced) , and what proportion of those babies survive the year (). If they survive, then this years babies are adults next year. So, is the number of adults next year contributed by fertility contributions of adults from this year.

Now, some simple algebra yields the equation for :



OK, but how do you get there from here, using a matrix approach? The answer is - for a simple annual life cycle diagram, you don't. You need to set things up as a seasonal life-cycle. Doing so also makes Pulliam's assumptions clear. Consider the following seasonal life-cycle:



Pre-breeding, year t - single-age class (adults). Post-breeding, 2 age classes (adults and babies). A key assumption here is that there is no mortality of adults during breeding (thus, the arc connecting the two adult nodes between the pre- and post-breeding censuses is parameterized with 1). Adult fertility contribution between pre-breeding (year t) to post-breeding (year t) is given by the parameter. Both age-classes post-breeding (in year t) contribute to the adult e class in year (t+1). Adults post-breeding survive with probability , while babies survive with probability .

Setting this up as a seasonal matrix model, in the following PRE indicates the matrix mapping pre-breeding to post-breeding (in year t), while POST maps the transitions from post-breeding (year t) to pre-breeding (year t+1).



So, the projection from pre-breeding (t) to pre-breeding (t+1) is



which of course is the equation we've been considering all along. (In other words, the projection matrix is , with the single element being equal to . Since is the dominant eigenvalue of the projection matrix, and since the dominant eigenvalue of a matrix is simply the value of the single element of the matrix, then .

OK -- a somewhat lengthy aside. But it is important to know where it comes from, and what assumptions you need to make to 'get there from here'.

And while we're at at (i.e., as long as I'm on my soapbox.

You frequently see the following:

Nope. More correctly, it should be . The difference (1+r) is an approximation to . This approximation comes from the Taylor series for :



It is only really good if r is small (<0.07) - after that, the value of 1+r is biased low wrt to true - you'd need to include the second order term for better fit. But, if you have r, then simply use . The only possible reason to use the approximation is if you don't have even a simple calculator at your disposal, and you want to express projected growth is terms of , and not r.

Evan's comments and cautions are well taken. Still, it might be useful for folks to know that the function deltavar() works for all the empirical examples given in Appendix 2 of the Mark Manual on Delta Method. The following R scripts complete the calculations and give exactly the same answers reported for the three examples on pages B-11 through B-19.

I get different answers for example 3 on pages B-20 to B-22, but suspect there are errors in the Appendix. On page B-20, values of the covariate mass are reported as mean = 109.97 with a SD = 24.79. The values for mass-squared are given as mean = 12707.46 and SD = 5532.03, but I think the correct values should be mean = 12093.4 and SD = 614.54. If these corrections are made, the estimates of the parameter and variance come out a little differently.

Appendix 2 was a great set of materials for testing the delta method function, thanks again for putting this resource together. Feel free to add the R scripts to the Appendix if you think they'd be helpful. If you decide to add some parameter estimates to example 2 on reporting rates on page B-14, I can prepare a script for that one too.

Cheers, Brett.

# delta method
# bookkeeping
rm(list = ls())

# need to install package emdbook
library(emdbook)

# Page B-11 Example 1 variance of product of survival probabilities
# period = phi1*phi2*phi3
par.name = c("phi1", "phi2", "phi3")
par.mean = c(phi1=0.6109350, phi2=0.458263, phi3=0.4960239)
par.varcovar = matrix(c(0.0224330125, -0.0003945405, 0.0000654469,
-0.0003945405, 0.0099722201, -0.0002361998,
0.0000654469, -0.0002361998, 0.0072418858), nrow = 3, byrow = T)
(period=prod(par.mean))
(period.var = deltavar(phi1*phi2*phi3, meanval=par.mean, vars=par.name, Sigma=par.varcovar, verbose=T))
# check answers, period=0.138871, var(period) = 0.0025565

# Page B-14 Example 2
# No parameter estimates reported

# Page B-15 Example 3a variance of back-transformed estimates - simple
# phihat = exp(beta)/(1+exp(beta))
par.name = c("beta")
par.mean = c(beta=0.2648275)
par.var = (0.1446688)^2
(phihat=plogis(as.numeric(par.mean)))
(phihat.var = deltavar(exp(beta)/(1+exp(beta)), meanval=par.mean, vars=par.name, Sigma=par.var, verbose=T))
# check answers, phihat=0.5658226, var(phihat) = 0.001263

# Page B-18 Example 3b variance of back-transformed estimates - somewhat harder
# phihat = exp(-1*(int+beta))/(1+exp(-1*(int+beta)))
par.name = c("int", "beta")
par.mean = c(int= 0.4267863, beta=-0.5066372)
par.varcovar = matrix(c(0.0321405326, -0.0321581167,
-0.0321581167, 0.0975720877), nrow = 2, byrow = T)
(phihat=plogis(as.numeric(sum(par.mean))))
(phihat.var = deltavar(exp(-1*(int+beta))/(1+exp(-1*(int+beta))), meanval=par.mean, vars=par.name, Sigma=par.varcovar, verbose=T))
# check answers, phihat=0.48005, var(phihat) = 0.0040742678

# Page B-20 Example 4 variance of back-transformed estimates - a bit harder still
# phihat = exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2)))
(mass = 110)
(mass.mean = 109.97)
(mass.SD = 24.79)
(mass2 = mass^2)
(mass2.mean = mass.mean^2)
(mass2.SD = mass.SD^2)
(mass.z = (mass - mass.mean)/mass.SD)
(mass2.z = (mass2 - mass2.mean)/mass2.SD)
par.name = c("int", "beta1", "beta2")
par.mean = c(int= 0.256732, beta1=1.1750358*mass.z, beta2=-1.0554864*mass2.z)
par.varcovar = matrix(c(0.0009006921, -0.0004109710, 0.0003662359,
-0.0004109710, 0.0373887267, -0.0364250288,
0.0003662359, -0.0364250288, 0.0362776933), nrow = 3, byrow = T)
(phihat=plogis(as.numeric(sum(par.mean))))
(phihat.var = deltavar(exp(-1*(int+beta1+beta2))/(1+exp(-1*(int+beta1+beta2))), meanval=par.mean, vars=par.name, Sigma=par.varcovar, verbose=T))
# values below may be incorrect
# check answers, mass=110, phihat=0.5925, var(phihat) = 0.00007387
# check answers, mass=120, phihat=(not given), var(phihat) = 0.000074214

Thanks - very helpful. I'll be in touch offline if/when I decide to add R-scripts to the appendix.

B.K. Sandercock wrote:I get different answers for example 3 on pages B-20 to B-22, but suspect there are errors in the Appendix. On page B-20, values of the covariate mass are reported as mean = 109.97 with a SD = 24.79. The values for mass-squared are given as mean = 12707.46 and SD = 5532.03, but I think the correct values should be mean = 12093.4 and SD = 614.54. If these corrections are made, the estimates of the parameter and variance come out a little differently.

Nope, values are correct - se p. 13 of chapter 11: means of mass is 109.96803 (SD 24.79), mean of mass^2=12707.464 (SD=5532.032). Back-transformed estimate

So, differences are likely due to something else. I'll check further.