2
$\begingroup$

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,

enter image description here

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.

$\endgroup$
  • 2
    $\begingroup$ If I understand your problem correctly I would have thought survival analysis was closer to what you want. The CRAN task view on SurvivalAnalysis provides a long list of software. $\endgroup$ – mdewey May 3 '16 at 14:42
  • $\begingroup$ Do you just want a hazard that depends on covariates but that is constant over time for each tree (no age dependence)? If so you can do this with a glm with cloglog link-function and log(interval) as an offset variable. $\endgroup$ – Jarle Tufto Jun 6 '16 at 20:51
  • $\begingroup$ @Jarle thats precisely what I want! Can you respond to the question with an example of how to do this in R? Will this still account for the fact that my response variable in a binomial outcome (live or die, 0 or 1)? $\endgroup$ – colin Jun 7 '16 at 14:31
5
+50
$\begingroup$

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.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ this is perfect. Thank you for the detailed explanation and the implementation. $\endgroup$ – colin Jun 7 '16 at 16:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.