3
$\begingroup$

I am reading a paper in which the authors models tree survival (mortality). They go and remeasure tagged trees for decades to establish "survival functions" for the given tree species and note whether they are still alive or dead. Initially when the trees are young, they measure every year but as the trees get older the measurement interval gets larger. In order to account for the varying measurement intervals (labeled as $L$) they added it as an exponent to the logistic link function (which based on this paper is called a generalized logistic model: A Generalized Mixed Logistic Model for Predicting Individual Tree Survival Probability with Unequal Measurement Intervals ) by Yuqing Yang and Shongming Huang:

$$P_{ijk} = \left( \frac{e^c}{1 + e^c} \right)^{L_{jk}}$$

where $j$ is the probability of the $j^{th}$ tree at the $k^{th}$ measurement in the $i{th}$ plot, and $c$ is $\beta_{0}+\beta_{1}*DIAMETER + \beta_{1}*DIAMETER^2 + \beta_{3}*TEMPERATURE$.

My question is how can this be modeled in R? Without $L_{jk}$ as the exponent I would do:

glm(survival ~ DIAMETER + I(DIAMETER )^2 + TEMPERATURE, 
    family = binomial(link = "logit"))

Could it be as simple as adding an offset:

glm(survival ~ DIAMETER + I(DIAMETER )^2 + TEMPERATURE, 
    offset = log(L), family = binomial(link = "logit"))

?

There is a very similar question here but it doesn't help me understanding whether setting $L$ as an exponent in the model description, is the same as specifying the offset when defining and running the model in R.

EDIT

Thinking about this a bit more, it could be modeled as:

model <-
  glm(
    survival ~ DIAMETER + I(DIAMETER) ^ 2 + TEMPERATURE,
    family = binomial(link = "logit"),
    data = d
  )

# Calculate p0 using logistic link function and L as an exponent
d$pred <- predict(model, type = "link")
d$p0 <- (exp(d$pred) / (1 + exp(d$pred))) ^ d$L

Or is there yet a more approriate way?

EDIT 2

I did some more research and this publication by @BenBolker provides an answer I was looking for: https://rpubs.com/bbolker/logregexp

I leave it here if people come across the same question.

$\endgroup$
6
  • 2
    $\begingroup$ Have you come across this thread? Custom link function needed for generalised linear model in R $\endgroup$
    – dipetkov
    Commented Jan 25 at 19:35
  • $\begingroup$ @dipetkov thanks for the link - definitely helpful! $\endgroup$
    – Stefan
    Commented Jan 25 at 19:56
  • $\begingroup$ If you found an answer, it would be best to post it as an actual answer so that the question actually has an "official" answer. $\endgroup$ Commented Jan 26 at 16:51
  • $\begingroup$ @COOLSerdash I was thinking about it but paraphrasing everything from Ben's post is a bit time consuming and simply posting a short note with a link isn't a good answer, no? What would you do? I can paste his post as an answer and making clear that it is a complete copy of his original post... $\endgroup$
    – Stefan
    Commented Jan 26 at 16:55
  • 1
    $\begingroup$ Personally, I'd post the link to Ben's post as well as a few words of how I adapted the solution to my specific dataset and problem. $\endgroup$ Commented Jan 26 at 17:23

1 Answer 1

2
$\begingroup$

Just wanted to share what I did in the end as suggested by @COOLSerdash :

First, I took @BenBolker 's hack for a custom power-logistic function which can be found here: https://rpubs.com/bbolker/logregexp

This is how it looks like (straight copy from the link):

library(MASS)
logexp <- function(exposure = 1) {
  ## hack to help with visualization, post-prediction etc etc
  get_exposure <- function() {
    if (exists("..exposure", env=.GlobalEnv))
      return(get("..exposure", envir=.GlobalEnv))
    exposure
  }
  linkfun <- function(mu) qlogis(mu^(1/get_exposure()))
  ## FIXME: is there some trick we can play here to allow
  ##   evaluation in the context of the 'data' argument?
  linkinv <- function(eta) plogis(eta)^get_exposure()
  logit_mu_eta <- function(eta) {
    ifelse(abs(eta)>30,.Machine$double.eps,
           exp(eta)/(1+exp(eta))^2)
  }
  mu.eta <- function(eta) {       
    get_exposure() * plogis(eta)^(get_exposure()-1) *
      logit_mu_eta(eta)
  }
  valideta <- function(eta) TRUE
  link <- paste("logexp(", deparse(substitute(exposure)), ")",
                sep="")
  structure(list(linkfun = linkfun, linkinv = linkinv,
                 mu.eta = mu.eta, valideta = valideta, 
                 name = link),
            class = "link-glm")
}

Note: as you can see there is a hack to allow for help with visualization, post-prediction etc. since glm() and glmer() do not evaluate the exposure variable in the context of the data argument.

So this is how I adapted this solution to my problem using glm(), glmer() and ggeffects::ggpredict().

First for glm():

library(lme4)
library(ggeffects)

m <- glm(
  survival ~ scale(dbh) + scale(dbh ^ 2) + scale(temp)  ,
  family = binomial(logexp(d$measurement_interval)),
  data = d
)
# pull in the mean for the exposure variable. In my case `measurement_interval`
..exposure <- mean(d$measurement_interval)

# run model predictions (that's were `..exposure` is needed)
pred_dbh <- ggpredict(m, terms = c("dbh [all]"))
pred_temp <- ggpredict(m, terms = c("temp [all]"))

# remove `..exposure` to avoid issues when using different exposure values
rm(..exposure)

# plot predictions
plot(pred_dbh)
plot(pred_temp)

The same also works for glmer():

m2 <- glmer(
  survival ~ scale(dbh) + scale(dbh ^ 2) + scale(temp) + (1 | plot) ,
  family = binomial(logexp(d$measurement_interval)),
  data = d
)
# pull in the mean for the exposure variable. In my case `measurement_interval`
..exposure <- mean(d$measurement_interval)

# run model predictions (that's were `..exposure` is needed)
pred_dbh2 <- ggpredict(m2, terms = c("dbh [all]"))
pred_temp2 <- ggpredict(m2, terms = c("temp [all]"))

# remove `..exposure` to avoid issues when using different exposure values
rm(..exposure)

# plot predictions
plot(pred_dbh2)
plot(pred_temp2)

It is generally much easier to fit these types of models using a Bayesian non-linear approach for example with the brms package. The reason it didn't work for me though is because of the large size of the dataset (+1M observations) and it takes days to fit those models.

Anyway, here the Bayesian approach that I tried:

library(brms)

m3 <-
  brm(
    bf(
      # modeling survival using the "power-logistic function"
      survival ~ (exp(eta) / (1 + exp(eta))) ^ measurement_interval,
      eta ~ scale(dbh) + scale(I(dbh ^ 2)) + scale(temp) + (1 | plot),
      nl = TRUE,
      # set link to "identity" since we defined the "power-logistic function" above already
      family = bernoulli(link = "identity")
    ),
    data = d)
  )
# plot predictions (via `brms::conditional_effects()`)
conditional_effects(m3)

# get predictions (via `ggeffects::ggpredict()`)
pred_dbh3 <- ggpredict(m3, terms = "dbh [all]"))
pred_temp3 <- ggpredict(m3, terms = "temp [all]"))

# plot predictions
plot(pred_dbh3)
plot(pred_temp3)

I hope that's useful!

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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