CODEWITHSUNDEEP
rlme4
tomka
tomka

Reputation: 2638

Estimating risk ratio instead of odds ratio in mixed effect logistic regression in `R`

glmer is used to estimate effects on the logit scale of y when the data are clustered. In the following model

fit1 = glmer(y ~ treat + x + ( 1 | cluster), family = binomial(link = "logit"))

the exp of the coefficient of treat is the odds ratio of a binary 0-1 treatment variable, x is a covariate, and cluster is a clustering indicator across which we model a random effect (intercept). A standard approach in glm's to estimate risk ratios is to use a log link instead, i.e. family=binomial(link = "log"). However using this in glmer I get error

Error in (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L, maxit = 100L,  : 
  (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate

after calling

fit1 = glmer(y ~ treat + x + ( 1 | cluster), family = binomial(link = "log"))

A web search revealed other people had similar issues with the Gamma family.

This seems to be a general problem as the reproducible example below demonstrates. My question thus is: how can I estimate risk ratios using a mixed effect model like glmer?

Reproducible Example The following code simulates data that replicates the problem.

n = 1000                            # sample size
m = 50                              # number of clusters
J = sample(1:m, n, replace = T)     # simulate cluster membership
x = rnorm(n)                        # simulate covariate
treat = rbinom(n, 1, 0.5)           # simulate random treatment
u  = rnorm(m)                       # simulate random intercepts
lt = x + treat + u[J]               # compute linear term of logistic mixed effect model
p  = 1/(1+exp(-lt))                 # use logit link to transform to probabilities
y  = rbinom(n,1,p)                  # draw binomial outcomes
d  = data.frame(y, x, treat)

# First fit logistic model with glmer
fit1  = glmer( y ~ treat + x + (1 | as.factor(J)), 
               family = binomial(link = "logit"), data  = d) 
summary(fit1)

# Now try to log link    
fit2  = glmer( y ~ treat + x + (1 | as.factor(J)), 
               family = binomial(link = "log"), data  = d)

Upvotes: 2

Views: 798

Answers (1)

iacob
iacob

Reputation: 24201

This error is returned due to your model producing values > 1:

  • PIRLS step-halvings failed to reduce deviance in pwrssUpdate ...
    • When using lme4 to fit GLMMs with link functions that do not automatically constrain the response to the allowable range of the distributional family (e.g. binomial models with a log link, where the estimated probability can be >1, or inverse-Gamma models, where the estimated mean can be negative), it is not unusual to get this error. This occurs because lme4 doesn’t do anything to constrain the predicted values, so NaN values pop up, which aren’t handled gracefully. If possible, switch to a link function to one that constrains the response (e.g. logit link for binomial or log link for Gamma).

Unfortunately, the suggested workaround is to use a different link function.

The following paper surveys a number of alternative model choices for calculation for [adjusted] relative risk:

Upvotes: 2

Related Questions