analyzing binary mortality data collected at different remeasurement intervals. Bonus for R implementation

Question

I am replicating an analysis that models tree mortality data. Data are structured such that forest sites are revisted at some random interval, which is recorded. It is then determined if a tree lived or died over that random interval, generating 0 1 mortality data (if a tree dies, it gets a 1 in the dependent variable). The interval between initial and final observation varies continuously, from 5-15 years. This is relevant, as the more time that passes, the more likely a tree will die.

Here are some pseudo data for R:

mort <- c(0,1,0,0,0,0,0,1,1,1,1,0,1,0,0,0,0,0,0,0,1,0)
interval <- runif(length(mort), 5, 15)
pollution <- rnorm(length(mort), 25,5)
data<- data.frame(mort, interval, pollution)

I am trying to replicate an analysis which uses a logistic regression model for binary mortality data using the the logit transformation. Authors then model how pollution affects tree mortality rates. In the manuscript the authors write, "because recensus is not annual, we relate annual mortality probability, pi, of tree i to the observed binomial data on whether that tree lived or died Mi via a Bernoulli likelihood,

where ti is the time interval between successive censuses."

My question: How would I implement this using the glm function, or something analagous, in R? Note: I understand modeling this as a hazard function would also be appropriate, but it is not what I am interested in.

Answer 1

Ok, if each tree $i$ has an age-independent hazard $\lambda_i$, then its lifetime $T_i$ follows an exponential distribution such that the probability that the tree is dead after a time interval of length $t_i$ becomes \begin{equation} p_i = P(T_i \le t_i) = 1-\exp(-\lambda_i t_i). \end{equation} If we in turn assume that pollution has an additive effect on the log of the hazard, that is, $\log\lambda_i=\beta_0 + \beta_1 x_i$, we get the overall model \begin{equation} p_i = 1-\exp(-\exp(\beta_0 + \beta_1 x_i + \log t_i)), \end{equation} or \begin{equation} \log(-\log(1-p_i) = \beta_0 + \beta_1 x_i + \log t_i \end{equation} which is a glm with cloglog link-function and $\log t_i$ as offset.

This is fitted in R by doing

glm(mort ~ pollution, offset=log(interval), binomial(link="cloglog"))

This model is similar to the one you're trying to reproduce in that survival probabilities $e^{-\lambda_i t_i}=(e^{-\lambda_i})^{t_i}$ decrease with powers of annual survival probabilities $e^{-\lambda_i}$ but the exact effect of pollution is slightly different from what what follows from the logit link assumption.