Delta Method et Var/Cov Matrix

[cooch]

At some point tonite, I'll go through and make sure all pages for that example have the same parameter estimates in place....

Done (and just uploaded). I also added some text earlier in the Appendix alerting the reader that the optimization method, and starting values, can influence the values of the estimates and VC matrix MARK returns.

[B.K. Sandercock]

Appendix B of the Mark manual illustrates use of the delta method for calculating variance of parameters that are a function of multiple parameters or are on different scales (logit vs. probability).
The following R script will complete the calculations for all of the different examples in Appendix B.
With thanks to Evan Cooch and Ben Bolker for helping me to develop the script.
Note that the page numbers may change if any further revisions are made to Appendix B in future editions of the manual.

Code: Select all

# delta method for Appendix B of Mark manual
# Brett K. Sandercock and Ben M. Bolker
# need to install package emdbook
library(emdbook)

# Page B-11 Example 1 variance of product of survival probabilities
# period = 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 = TRUE)
(period=prod(par.mean))
(period.var = deltavar(phi1*phi2*phi3, meanval=par.mean, Sigma=par.varcovar, verbose = FALSE))
# check answers on pages B-11 and B-14, period=0.138871, var(period) = 0.0025565

# Page B-14 Example 2 variance of estimate of reporting rate
# No parameter estimates reported

# Page B-15 Example 3 variance of back-transformed estimates - simple
# phihat = exp(beta)/(1+exp(beta))
par.mean = c(beta=0.2648275)
par.var = (0.1446688)^2
(phihat=plogis(as.numeric(par.mean)))
# (phihat=1/(1+exp(-par.mean)))  # alternative
(phihat.var = deltavar(exp(beta)/(1+exp(beta)), meanval=par.mean, Sigma=par.var, verbose=FALSE))
(phihat.SE = sqrt(phihat.var))
# check answers on pages B-15 and B-16, phihat=0.5658226, var(phihat) = 0.001263, SE(phihat) = 0.0355404

# Page B-18 Example 4 variance of back-transformed estimates - somewhat harder
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 = TRUE)
(phihat=plogis(as.numeric(sum(par.mean))))
(phihat.var = deltavar(exp(int+beta)/(1+exp(int+beta)), meanval=par.mean, vars=par.name, Sigma=par.varcovar, verbose=FALSE))
(phihat.SE = sqrt(phihat.var))
# check answers on pages B-18 and B-20, phihat=0.48005, var(phihat) = 0.0040742678, SE(phihat) = 0.0638300

# Page B-21 Example 5 variance of back-transformed estimates - harder still
#mass = 110
mass = 120
mass.mean = 109.9680; mass.SD = 24.7926
mass2 = mass^2; mass2.mean = 12707.4640; mass2.SD = 5532.0322
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.2567341, beta1=1.1750463, beta2=-1.0554957)
par.varcovar = matrix(c(0.0009006967, -0.0004110129, 0.0003662771,
-0.0004110129, 0.0373928740, -0.0364201221,
0.0003662771, -0.0364291221, 0.0362817338), nrow = 3, byrow = TRUE)
(phihat.logit = par.mean %*% c(1, mass.z, mass2.z))
(phihat = plogis(phihat.logit))
(phihat.var = deltavar(exp(int+beta1*mass.z+beta2*mass2.z)/(1+exp(int+beta1*mass.z+beta2*mass2.z)), meanval=par.mean, Sigma=par.varcovar, verbose=FALSE))
(phihat.SE = sqrt(phihat.var))
# check answers on pages B-22 and B-23, mass=110, phihat.logit = 0.3742, phihat=0.5925, var(phihat) = 0.000073867, SE(phihat) = 0.0085946
# check answers on page B-24, mass=120, phihat=(not given), var(phihat) = 0.000074246, SE(phihat) = 0.0086166

# Page B-25 Delta method and model averaging
# period = phi1*phi2*phi3
par.mean = c(phi1=0.5673, phi2=0.5332, phi3=0.5332)
par.varcovar = matrix(c(0.0019410083, 0.0001259569, 0.0001259569,
0.000125569, 0.0033727452, 0.0033727423,
0.0001259569, 0.0033727423, 0.0033727452), nrow = 3, byrow = TRUE)
(period=prod(par.mean))
(period.var = deltavar(phi1*phi2*phi3, meanval=par.mean, Sigma=par.varcovar, verbose = FALSE))
# check answer on pages B-25 and B-26, period=0.1613, var(period) = 0.001435



[cooch]

The following R script will complete the calculations for all of the different examples in Appendix B.

Thanks Brett. Very useful and generous contribution.

[B.K. Sandercock]

Appendix C for the Mark manual describes a second function that can be used to apply the delta method.

Options for delta method calculations in Program R include:
deltavar() function in package emdbook (used above)
deltamethod() function in package msm (see page C-101)

[GHORTON]

I am attempting to estimate the SE (and 95% CI) for an estimate of cumulative apparent survival. The point estimate for Phi_cum comes from the product of 6 Phi-hats from a 7 occasion CJS model (SIN link). To estimate the variance of Phi_cum, I used the approach in Example 1 starting on B-19 of the MARK manual. I first made these computations using the Dipper data (Phi(.) p(.) logit) using my own code as well as with the Sandercock/Bolker R code and came up with the same values for Phi_cum and var-hat in both cases (which also means it matched what is in the MARK manual).

As pointed out in the sidebar at the bottom of B-24 (which I realize is for a single point estimate and not the product of multiple estimates), it seems the next step would be to take the square root of var-hat for Phi_cum to get SE-hat and use it to construct the 95% CI. However, also based on the sidebar at the bottom of B24 (i.e., the 95% CI should not be symmetrical for a single estimate of Phi-hat if different than 0.5 and it should fall between 0 and 1) it seems that this might be incorrect(?). What makes me uncertain is given that there is 1 variance calculated using this approach, the CI will be symmetrical and therefore could fall outside the 0-1 range. Perhaps I should be using the beta estimates but if so how? Any guidance would be most appreciated.

[cooch]

This is a slight hijack of the thread, but subject/intent is close enough I'll keep it in place.

Basically, you want a 95% CI. Which leads to the question -- why? Give what a frequentist CI is, that is not always an easy questyion to answer. It might be simpler to simply use MCMC, and calculate a highest posterior density for the product. This is easy, and lets you say 'MCMC' in your methods, which gives you street cred with the cool kids who think that everything needs to be done using a Bayesian approach to things. More importantly, its easy, and robust. This approach is dscussed last few pages of Appendix B.

[cooch]

If you really want to use the Delta method, here is the 'addendum example' from Appendix B (product of phi1, phi2, and phi3 from male dipper done as an R script. Fairly self-explanatory, with a bit of study. The beta estimates come from model phi(t)p(.), based on an identity design matrix:

Code: Select all

# Given beta values on the logit scale
beta_1 <- 0.4512440
beta_2 <- -0.1673375
beta_3 <- -0.0159047

# Given variance-covariance matrix
var_cov_matrix <- matrix(c(
  0.397057267, -0.0066860755, 0.001101428,
  -0.0066860755, 0.161802478, -0.0038059824,
  0.001101428, -0.0038059824, 0.1158848706
), nrow = 3, byrow = TRUE)

# Step 1: Convert beta parameters to the probability scale
p1 <- exp(beta_1) / (1 + exp(beta_1))
p2 <- exp(beta_2) / (1 + exp(beta_2))
p3 <- exp(beta_3) / (1 + exp(beta_3))

# The product of probabilities
p <- p1 * p2 * p3

# Step 2: Compute the gradients of g(beta_1, beta_2, beta_3)
grad_g_beta_1 <- (exp(beta_1) / (1 + exp(beta_1))^2) * p2 * p3
grad_g_beta_2 <- (exp(beta_2) / (1 + exp(beta_2))^2) * p1 * p3
grad_g_beta_3 <- (exp(beta_3) / (1 + exp(beta_3))^2) * p1 * p2

# Gradient vector
gradients <- c(grad_g_beta_1, grad_g_beta_2, grad_g_beta_3)

# Step 3: Use the delta method to calculate the variance of the product
var_product <- t(gradients) %*% var_cov_matrix %*% gradients

# Step 4: Calculate the 95% confidence interval
ci_lower <- max(0, p - 1.96 * sqrt(var_product))  # constrained to [0, 1]
ci_upper <- min(1, p + 1.96 * sqrt(var_product))  # constrained to [0, 1]

# Output the calculated results
list(product = p, CI_lower = ci_lower, CI_upper = ci_upper)



Thus, the 95% CI for the product on the real probability scale is approximately [0.040,0.238].

[cooch]

Which, I'll note, is very similar to the 95% HPD from the MCMC approach: [0.05256, 0.24017]. Not surprising, perhaps. So you can decide which is 'cooler' (MCMC or Delta method), since most of the time results will largely be the same.

[GHORTON]

Thanks for your reply. I will likely stick with the 95% CI but may try the MCMC approach as well.

Because some of the phi-hats from my fish data set are near 1, I would like to use the sin link instead of the logit. Using the dipper data, model phi(t)p(.), I was not able to get the variance from the sin model (my adaptation of your R code, see below) to match the variance from the logit model (your original R code with the logit). I believe the problem is with the gradient calculations. Would you be able send me the gradient formula for the sin link?

Again, thanks again for your help.

Code: Select all

#This is from phi(t)p(.) for the dipper data using the sin link and identity matrix
# Given beta values on the sin scale
beta_1 <- 0.2237318
beta_2 <- -0.0835713
beta_3 <- -0.0079523

# Given variance-covariance matrix
var_cov_matrix <- matrix(c(
  0.0943779291, -0.0016241721, 0.0002684908,
  -0.0016241721, 0.0401687802, -0.0009481403,
  0.0002684908, -0.0009481403, 0.0289693842
), nrow = 3, byrow = TRUE)

# Step 1: Convert beta parameters to the probability scale
p1 <- (sin(beta_1)+1)/2
p2 <- (sin(beta_2)+1)/2
p3 <- (sin(beta_3)+1)/2

# The product of probabilities
p <- p1 * p2 * p3

# Step 2: Compute the gradients of g(beta_1, beta_2, beta_3)
grad_g_beta_1 <- ((sin(beta_1)+1)^2 / 2) * p2 * p3 #???
grad_g_beta_2 <- ((sin(beta_2)+1)^2 / 2) * p1 * p3 #???
grad_g_beta_3 <- ((sin(beta_3)+1)^2 / 2) * p1 * p2 #???

# Gradient vector
gradients <- c(grad_g_beta_1, grad_g_beta_2, grad_g_beta_3)

# Step 3: Use the delta method to calculate the variance of the product
var_product <- t(gradients) %*% var_cov_matrix %*% gradients

# Step 4: Calculate the 95% confidence interval
ci_lower <- max(0, p - 1.96 * sqrt(var_product))  # constrained to [0, 1]
ci_upper <- min(1, p + 1.96 * sqrt(var_product))  # constrained to [0, 1]

# Output the calculated results
list(product = p, CI_lower = ci_lower, CI_upper = ci_upper)


[cooch]

In other words, you don't know how to take a derivative of a function?

For the logit

If

then



If you look carefully, that is what is in the R code I posted, for each for the beta parameters in turn. You're using this derivative to form the gradient vector -- i.e., you're using the derivative of the logit transform.

For the sin, consult the table of transforms and back-transforms in the -sidebar- at the top of p. 22 in Chapter 6 of the MARK book. The back-transform on the sin link (again, assuming beta estimates from identity matrix) is



So, you simply differentiate g with respect to (as was the case for the logistic, above).



The rest is up to you.

Incidentally, if you have boundary estimates (which seems to be motivating your interest in using the sin link), then the MCMC approach might actually be more robust here.