Mixed regression models and custom link functions in R?

Question

It seems like the current revision of lmer does not allow for custom link functions.

If one needs to fit a logistic linear mixed effect model with a custom link function what options are available in R?
If none - what options are available in other statistics/programming packages?
Are there conceptual reasons lmer does not have custom link functions, or are the constraints purely pragmatic/programmatic?

It is unfortunate that nobody could chime in with a package where this is easy to do. However, it turned out that adding the appropriate functions to lme4 wasn't as difficult as I had feared. Now the only challenge that remains is to find the "correct" link function (see: http://stats.stackexchange.com/questions/1583/appropriate-link-function-for-2afc-data). So I'll accept ars' answer, but someone may want to ask this again in the future since next time the answer may be different. — russellpierce, Aug 12 '10 at 07:09
An example from Ben Bolker: https://stackoverflow.com/a/15935574/2727349 — swihart, Jul 08 '17 at 00:28
@swihart that answer corresponds to custom link functions for glm, lme4::glmer() is a different enough situation the same answer does not apply I'm afraid (unless big changes have happened to the glmer internals... Which maybe is true) — russellpierce, Jul 08 '17 at 00:32
... and is true, see my answer below along with a link to a relevant answer by Ben. — russellpierce, Jul 25 '17 at 12:07

score 3 · Answer 1 · answered Aug 08 '10 at 19:26

Douglas Bates addressed this on the sig-ME list a while back:

using glmer with user-defined link function

I'm not aware of significant changes since, but his recommendation (using a quasi family with specified link and variance) might be of use. Hopefully this addresses your first and third questions. I'm not aware of other packages - sorry.

score 2 · Answer 2 · answered Jul 12 '17 at 01:05

repeated::gnlmix()
PROC NLMIXED -- If the link or its inverse can be represented with the functions available in SAS.
I need to dig more into the code.

Some repeated::gnlmix() examples:

library(repeated)
set.seed(99)
dat <- data.frame(id = c(1,1,1, 2,2,2, 3,3,3, 4,4,4),
                   y = c(1,0,1, 0,1,1, 0,0,0, 1,1,1),
                   x = rnorm(12)
                    )
attach(dat)
y_cbind <- cbind(y, 1-y)

# probit regression
gnlmix(y=y_cbind,
       distribution = "binomial",
       nest = id,
       mixture = "normal",
       random = "b_id",
       mu = ~ pnorm(beta0 + beta1*x + b_id),
       pmu = list(beta0=0, beta1=0),
       pmix = log(4)
)
# Location parameters:
#        estimate      se
# beta0    0.2825  0.7878
# beta1   -0.7628  1.1498
# 
# Mixing dispersion parameter:
#    estimate     se
#      0.4249  2.072

# cauchy link
gnlmix(y=y_cbind,
       distribution = "binomial",
       nest = id,
       mixture = "normal",
       random = "b_id",
       mu = ~ pcauchy(beta0 + beta1*x + b_id),
       pmu = list(beta0=0, beta1=0),
       pmix = log(4)
)
# Location parameters:
#        estimate     se
# beta0    0.6973  1.791
# beta1   -1.2328  2.812
# 
# Mixing dispersion parameter:
#    estimate     se
#       1.805  3.038

# logistic regression
gnlmix(y=y_cbind,
       distribution = "binomial",
       nest = id,
       mixture = "normal",
       random = "b_id",
       mu = ~ plogis(beta0 + beta1*x + b_id),
       pmu = list(beta0=0, beta1=0),
       pmix = log(4)
)
# Location parameters:
#        estimate     se
# beta0    0.4981  1.356
# beta1   -1.2740  2.036
# 
# Mixing dispersion parameter:
#    estimate     se
#       1.496  2.177

# custom link / inverse (logistic regression)
custom_inv <- function(eta) exp(eta) / (exp(eta) + 1)
gnlmix(y=y_cbind,
       distribution = "binomial",
       nest = id,
       mixture = "normal",
       random = "b_id",
       mu = ~ custom_inv(beta0 + beta1*x + b_id),
       pmu = list(beta0=0, beta1=0),
       pmix = log(4)
)
# Location parameters:
#        estimate     se
# beta0    0.4981  1.356
# beta1   -1.2740  2.036
# 
# Mixing dispersion parameter:
#    estimate     se
#       1.496  2.177

# custom link / inverse 
custom_inv <- function(eta) log(exp(eta)+1) / (log(exp(eta)+1) + pi/2)
gnlmix(y=y_cbind,
       distribution = "binomial",
       nest = id,
       mixture = "normal",
       random = "b_id",
       mu = ~ custom_inv(beta0 + beta1*x + b_id),
       pmu = list(beta0=0, beta1=0),
       pmix = log(4)
)
# Location parameters:
#        estimate     se
# beta0     4.333  5.477
# beta1    -2.291  5.571
# 
# Mixing dispersion parameter:
#    estimate     se
#       3.879  2.386

For GLMs, see gnlm::gnlr() as shown here.

FWIW `lme4::glmer` hasn't had this limitation in a long time ... see my answer [here](https://stackoverflow.com/questions/19012128/user-defined-link-function-for-glmer-for-known-fate-survival-modeling) — Ben Bolker, Jul 12 '17 at 02:10

score 2 · Accepted Answer · answered Jul 12 '17 at 10:42

As pointed out by @Ben Bolker, this isn't a problem for versions of lme4 >= 1.0-4 (circa 09-2013). Link functions can now be specified in R code - which seems to answer 1) and 3) from above with the answer being that the constraints for link functions were due to how lme4 pre 1.0 was programmed as opposed to conceptual.

Mixed regression models and custom link functions in R?

3 Answers3