A common situation in ecology is a survival model with a binary outcome (0=die, 1=survive) where individuals (for this example think of individual nesting attempts by birds) differ in the number of days they are exposed to risk of mortality. To account for this, we use a modified logistic regression which incorporates the # of exposure days into the link function.
As described by Shaffer (2004):
“The daily survival rate is modeled in terms of x through the choice of an appropriate predictor function, which in our case should yield values between zero and one. As is done in logistic regression, we use the S-shaped logistic function:
The systematic component of our generalized linear model is then [s(x)]t. Next, we consider the function:
The above function is monotonic and differentiable with respect to θ, and it can be shown that g(θ) = β0 + β1x, which satisfies the criteria for a link function in a generalized linear model. Those three components: the binomial response distribution, the predictor function given in Expression 1, and the link function given in Expression 2, completely specify our generalized linear model. The model (hereafter “the logistic-exposure model”) is similar to the logistic regression model but differs in the form of the link function. The logistic-exposure link function contains an exponent (1/t) in the numerator and denominator that is not present in the logistic-regression link function. The exponent is necessary to account for the fact that probability of surviving an interval depends on interval length.”
Code for this link function is available on the web, and is also one of the example link functions described in R if you type "help(family)":
logexp <- function(days = 1)
{
linkfun <- function(mu) qlogis(mu^(1/days))
linkinv <- function(eta) plogis(eta)^days
mu.eta <- function(eta) days * plogis(eta)^(days-1) *
.Call("logit_mu_eta", eta, PACKAGE = "stats")
valideta <- function(eta) TRUE
link <- paste("logexp(", days, ")", sep="")
structure(list(linkfun = linkfun, linkinv = linkinv,
mu.eta = mu.eta, valideta = valideta, name = link),
class = "link-glm")
}
It works just fine in a model such as this:
glm(survive ~ date, family=binomial(link=logexp(days=dat$Days)),data=dat)
The problem I run into is when trying to use this custom link function in a GLMER model with the addition of a random effect (I have found one example online of this method being implemented here: http://rstudio-pubs-static.s3.amazonaws.com/4082_51aa699bd9f041c7b3f7cf7b9252f60c.html).
In our case we would like to include site as a random effect. Models are formulated the same way as the GLMs from before:
glmer(survive ~ date + (1|site), family=binomial(link=logexp(days=dat$Days)),data=dat)
However, now I get an error message:
Error in famType(glmFit$family) : unknown link: ‘logexp(3)’unknown link: ‘logexp(4)’unknown link: ‘logexp(3)’unknown link: ‘logexp(2)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(4)’unknown link: ‘logexp(3)’unknown link: ‘logexp(2)’unknown link: ‘logexp(1)’unknown link: ‘logexp(4)’unknown link: ‘logexp(5)’unknown link: ‘logexp(4)’unknown link: ‘logexp(3)’unknown link: ‘logexp(4)’unknown link: ‘logexp(5)’unknown link: ‘logexp(3)’unknown link: ‘logexp(4)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(3)’unknown link: ‘logexp(2)’unknown link: ‘logexp(1)’unknown link: ‘logexp(3)’unknown link: ‘logexp(1)’unknown link: ‘logexp(1)’unknown link: ‘logexp(1)’unknown link: ‘logexp(1)’unkn In addition: Warning message: In if (!(lTyp <- match(family$link, linkNms, nomatch = 0))) stop(gettextf("unknown link: %s", : the condition has length > 1 and only the first element will be used
The error message lists an unknown link for every row of data with the numbers corresponding to the exposure days for that nest visit (or row of data).
ex: The first 'logexp(3)' corresponds to the first row of data which has 3 exposure days.
Has anyone else been able to use this custom link function in a GLMER model? Or if not does anyone have an idea about what’s causing the error?
######UPDATE######
Many thanks to Ben Bolker for solving this problem. I updated to 3.0.2 and the latest version of lme4 and used the link function Ben's R related post (https://www.rpubs.com/bbolker/logregexp), which is this one:
library(MASS)
logexp <- function(exposure = 1)
{
linkfun <- function(mu) qlogis(mu^(1/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)^exposure
mu.eta <- function(eta) exposure * plogis(eta)^(exposure-1) *
.Call(stats:::C_logit_mu_eta, eta, PACKAGE = "stats")
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")
}
