1
$\begingroup$

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.

I am not sure if this is due to my data or a general problem. My question thus is: how can I estimate risk ratios using a mixed effect model like glmer, preferably in R?

I saw that stata has a function meglm which seems to allow this link and family. Not sure if this works there.

$\endgroup$
6
  • $\begingroup$ I'm not sure it makes sense to use a log link, since that will not constrain the response properly. Have you tried other packages such as brms and GLMMAdaptive` ? $\endgroup$ Mar 5, 2021 at 19:00
  • $\begingroup$ @RobertLong Thanks, I am also not absolutely sure. My idea was that I have to use the log link as I thought this is the general way it's done in glm. Maybe I am wrong? I cannot use Bayesian statistics for this analysis. GLMMAdaptive may be an option, I will check. So my general question is indeed how to do it rather than what's the cause of this specific error. $\endgroup$
    – tomka
    Mar 5, 2021 at 19:05
  • $\begingroup$ I would recommend that you stick with the odds ratio that is a natural output from logistic regression. $\endgroup$ Mar 5, 2021 at 20:00
  • $\begingroup$ Here is an example, but without random effects stats.stackexchange.com/questions/289008/… $\endgroup$ Sep 17, 2023 at 0:53
  • 1
    $\begingroup$ Changing from a probably well-fitting link function to a surely poorly-fitting one is dubious. Use the logistic model to get covariate-specific risk ratios as shown here. Unlike odds ratios, risk ratios are incapable of being constant over a wide range of baseline risk. Consider a variable that has a risk ratio of 2. This ratio cannot apply to someone with a base risk $> \frac{1}{2}$. $\endgroup$ Sep 17, 2023 at 11:52

2 Answers 2

3
$\begingroup$

This is an older question, but might benefit from some newer approaches using marginaleffects. The difference here is that we use post estimation techniques to compute the relative risk rather than have the relative risk be an estimand of the model.

Let's simulate some similar data and fit the model

library(tidyverse)
library(marginaleffects)
library(lme4)

set.seed(0)
n <- 100
cluster <- sample(1:3, size = n, replace=T)
txt <- rbinom(n, 1, 0.5)
x <- rnorm(n)

X <- model.matrix(~x+txt)
Z <- model.matrix(~factor(cluster)-1)
gamma <- rnorm(ncol(Z), 0, 0.25)
eta <- X%*% c(-1, 0.1, 1) + Z %*% gamma
p <- plogis(eta)
y <- rbinom(n, 1, p)


data <- tibble(x, txt, cluster, y)
fit <- glmer(y ~ txt + x + ( 1 | cluster), family = binomial(link = "logit")) 

Let's assume now we want to get the relative risk for a new cluster not in our data. We can do that with


avg_comparisons(
  fit, 
  variables = 'txt',
  by='x',
  newdata = datagrid(x=seq(-2, 2, 0.1), txt=0:1, cluster = -1),
  allow.new.levels=T,
  re.form = NULL,
  comparison = 'lnratioavg',
  transform = exp
  
)

This returns the relative risk of the treatment within strata defined by $x$. We can plot this easily using the dataframe returned, or the plot_comparisons1 function

plot_comparisons(
  fit, 
  variables = 'txt',
  by='x',
  newdata = datagrid(x=seq(-2, 2, 0.1), txt=0:1, cluster = -1),
  allow.new.levels=T,
  re.form = NULL,
  comparison = 'lnratioavg',
  transform = exp
  
)

enter image description here

To get the cluster relative risks, change the newdata and by argument as follows


plot_comparisons(
  fit, 
  variables = 'txt',
  by=c('x', 'cluster'),
  newdata = datagrid(x=seq(-2, 2, 0.1), txt=0:1, by='cluster'),
  comparison = 'lnratioavg',
  transform = exp
  
)

More info can be found here on the {marginaleffects} website

$\endgroup$
1
  • $\begingroup$ Nice, thanks for picking this up. $\endgroup$
    – tomka
    Apr 30 at 18:40
0
$\begingroup$

GLMMAdaptive solved this annoying problem for me, it seems to converge better for binomial with log link, but may need adjusting your control parameters, I set update_GH_every=50)

GLMMadaptive::mixed_model(fixed=dependent_var~independent_var+independent_var2
                      ,control = list(
                                     verbose=T,
                                    
                                     update_GH_every=40
                                     ),

                      random=~1|yourrandom_effect
                      na.action=na.fail,
                      data=youdata
                      family=binomial(link='log'),
         

)

$\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.