Solved – How to interpret the variance of random effect in a generalized linear mixed model

lme4-nlmer

In a logistic Generalized Linear Mixed Model (family = binomial), I don't know how to interpret the random effects variance:

Random effects:
 Groups   Name        Variance Std.Dev.
 HOSPITAL (Intercept) 0.4295   0.6554  
Number of obs: 2275, groups: HOSPITAL, 14

How do I interpret this numerical result?

I have a sample of renal trasplanted patients in a multicenter study. I was testing if the probability of a patient being treated with a specific antihypertensive treatment is the same among centers. The proportion of patients treated varies greatly between centers, but may be due to differences in basal characteristics of the patients. So I estimated a generalized linear mixed model (logistic), adjusting for the principal features of the patiens.
This are the results:

Generalized linear mixed model fit by maximum likelihood ['glmerMod']
 Family: binomial ( logit )
Formula: HTATTO ~ AGE + SEX + BMI + INMUNOTTO + log(SCR) + log(PROTEINUR) + (1 | CENTER) 
   Data: DATOS 

     AIC      BIC   logLik deviance 
1815.888 1867.456 -898.944 1797.888 

Random effects:
 Groups   Name        Variance Std.Dev.
 CENTER (Intercept) 0.4295   0.6554  
Number of obs: 2275, groups: HOSPITAL, 14

Fixed effects:
                           Estimate Std. Error z value Pr(>|z|)    
(Intercept)               -1.804469   0.216661  -8.329  < 2e-16 ***
AGE                       -0.007282   0.004773  -1.526  0.12712    
SEXFemale                 -0.127849   0.134732  -0.949  0.34267    
BMI                        0.015358   0.014521   1.058  0.29021    
INMUNOTTOB                 0.031134   0.142988   0.218  0.82763    
INMUNOTTOC                -0.152468   0.317454  -0.480  0.63102    
log(SCR)                   0.001744   0.195482   0.009  0.99288    
log(PROTEINUR)             0.253084   0.088111   2.872  0.00407 **

The quantitative variables are centered.
I know that the among-hospital standard deviation of the intercept is 0.6554, in log-odds scale.
Because the intercept is -1.804469, in log-odds scale, then probability of being treated with the antihypertensive of a man, of average age, with average value in all variables and inmuno treatment A, for an "average" center, is 14.1 %.
And now begins the interpretation: under the assumption that the random effects follow a normal distribution, we would expect approximately 95% of centers to have a value within 2 standard deviations of the mean of zero, so the probability of being treated for the average man will vary between centers with coverage interval of:

exp(-1.804469-2*0.6554)/(1+exp(-1.804469-2*0.6554))

exp(-1.804469+2*0.6554)/(1+exp(-1.804469+2*0.6554))

Is this correct?

Also, how can I test in glmer if the variability between centers is statistically significant?
I used to work with MIXNO, an excellent software of Donald Hedeker, and there I have an standard error of the estimate variance, that I don't have in glmer.
How can I have the probability of being treated for the "average" man in each center, with a confidene interval?

Thanks

Best Answer

It's probably most helpful if you show us more information about your model, but: the baseline value of the log-odds of whatever your response is (e.g. mortality) varies across hospitals. The baseline value (the per-hospital intercept term) is the log-odds of mortality (or whatever) in the baseline category (e.g. "untreated"), at a zero value of any continuous predictors. That variation is assumed to be Normally distributed, on the log-odds scale. The among-hospital standard deviation of the intercept is 0.6554; the variance (just the standard deviation squared -- not a measure of the uncertainty of the standard deviation) is $0.6554^2=0.4295$.

(If you clarify your question/add more detail about your model I can try to say more.)

update: your interpretation of the variation seems correct. More precisely,

cc <- fixef(fitted_model)[1] ## intercept
ss <- sqrt(unlist(VarCorr(fitted_model))) ## random effects SD
plogis(qnorm(c(0.025,0.975),mean=cc,sd=ss))

should give you the 95% interval (not really quite confidence intervals, but very similar) for the probabilities of a baseline (male/average age/etc.) individual getting treated across hospitals.

For testing the significance of the random effect, you have a variety of choices (see http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html for more information). (Note that the standard error of a RE variance is usually not a reliable way to test significance, since the sampling distribution is often skewed/non-Normal.) The simplest approach is to do a likelihood ratio test, e.g.

pchisq(2*(logLik(fitted_model)-logLik(fitted_model_without_RE)),
       df=1,lower.tail=FALSE)/2

The final division by 2 corrects for the fact that the likelihood ratio test is conservative when the null value (i.e. RE variance=0) is on the boundary of the feasible space (i.e. the RE variance cannot be <0).

Related Solutions

Solved – Where is the correlation parameter in the linear mixed-effect model equation

Maybe this will help you: in lme4 you can specify correlated intercept and slope:

y ~ x + (x | g) that translates to: y ~ 1 + x + (1 + x | g)

where y is a response variable and x is a predictor, while g is some grouping variable for random effects. lme4 by default assumes that the terms could be correlated, but you can also define them as uncorrelated:

y ~ x + (x || g) that translates to: y ~ 1 + x + (1 | g) + (0 + x | g)

In the second case intercept and slope are treated as independent, i.e. their correlation is constrained to be zero. So yes, correlated parameter is a part of variance-covariance matrix. Check the article by Bates et al. (in press) for more information on this.

Model Fit – Advice on Improving the Model Fit of a Mixed Model on Repeated Measurements

If you're willing to use a Bayesian approach instead, you could try brms. The model seems to converge then:

Input

library(dplyr)
library(brms)
mydata <- mydata %>% mutate(age = c(scale(age)))
m <- brm(bf(outcome_log ~ age + sex*age + smoking*age + obesity*age + diab*age + hypt*age + hyperchol*age + ckd*age + (1+age | pat_id), 
   hu ~ sex + smoking + obesity + diab + hyperchol + ckd + (1|pat_id)), 
   data=mydata,
   family=hurdle_lognormal())
summary(m)

Output

 Family: hurdle_lognormal 
  Links: mu = identity; sigma = identity; hu = logit 
Formula: outcome_log ~ age + sex * age + smoking * age + obesity * age + diab * age + hypt * age + hyperchol * age + ckd * age + (1 + age | pat_id) 
         hu ~ sex + smoking + obesity + diab + hyperchol + ckd + (1 | pat_id)
   Data: mydata (Number of observations: 304) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Group-Level Effects: 
~pat_id (Number of levels: 146) 
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)          0.29      0.03     0.23     0.35 1.00     1156     1968
sd(age)                0.07      0.05     0.00     0.20 1.01      403      520
sd(hu_Intercept)       7.03      1.97     4.09    11.48 1.00      874     1446
cor(Intercept,age)    -0.12      0.47    -0.90     0.86 1.00     1848     1946

Population-Level Effects: 
                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept                  1.47      0.12     1.24     1.70 1.00     1204     2014
hu_Intercept               4.68      2.54     0.56    10.67 1.00      985     1301
age                        0.02      0.14    -0.25     0.29 1.00     1192     1922
sexMale                    0.12      0.08    -0.03     0.28 1.00      984     1678
smokingNeversmoker        -0.33      0.08    -0.49    -0.17 1.01     1130     1588
obesityYes                -0.10      0.10    -0.30     0.09 1.00     1079     1792
diabYes                    0.30      0.13     0.05     0.57 1.00     1161     1839
hyptYes                    0.11      0.07    -0.04     0.25 1.00     1055     1835
hypercholYes               0.07      0.10    -0.14     0.27 1.00     1187     1871
ckdBelow90                -0.02      0.08    -0.17     0.14 1.01      864     1225
age:sexMale                0.04      0.08    -0.12     0.21 1.00     1406     2314
age:smokingNeversmoker     0.23      0.10     0.04     0.43 1.00     1386     2032
age:obesityYes            -0.08      0.14    -0.35     0.19 1.00     1780     2447
age:diabYes                0.17      0.15    -0.13     0.47 1.00     1943     2514
age:hyptYes               -0.16      0.09    -0.33     0.00 1.00     1554     2319
age:hypercholYes           0.16      0.13    -0.09     0.41 1.00     1372     2314
age:ckdBelow90            -0.03      0.09    -0.21     0.14 1.00     1390     1857
hu_sexMale                -2.93      1.76    -6.87     0.07 1.00     1174     1367
hu_smokingNeversmoker      2.83      1.77    -0.24     6.99 1.00     1119     1477
hu_obesityYes             -0.22      2.42    -5.20     4.35 1.00     1258     1722
hu_diabYes               -26.76     17.68   -71.71    -5.47 1.00     1516      912
hu_hypercholYes           -6.18      2.58   -12.12    -2.22 1.00     1075     1809
hu_ckdBelow90             -1.64      1.74    -5.20     1.47 1.00     1051     1638

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.15      0.01     0.13     0.17 1.00     1871     2724

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

This appears to identify some of the problem. In particular, it looks like there is a (quasi-) separation problem with diab in the hurdle equation. It's coefficient of -26 on the logic scale is really far away from zero. My guess is that this is one of the reasons that the frequentist GLMM is not converging with diab in the model. You can see the same, though to a lesser degree, with hyperchol.

Edit - diagnosing the problem

Essentially, the hurdle model is estimating a login on the 0 vs. not zero on the outcome. Then estimating a different model (in this case, a log-normal GLM) on the non-zero observations. One potential problem in logistic regression modes (or any model for binary dependent variables) is separation. The simplest example of this problem is when there is no variation in y for a particular category in a categorical x. Consider the cross-tabulation between zero/not-zero on the outcome and the diab variable.

table(
  factor(I(mydata$outcome == 0), 
     levels=c(FALSE, TRUE), 
     labels=c("Not Zero", "Zero")), 
  mydata$diab)

#           No Yes
# Not Zero 178  19
# Zero     107   0

If you think about what's happening here. The model is trying to recover the log of the odds ratio for Yes vs No observations:

If we substitute the probabilities in from the cross-tab, we would get:

The odds ratio would be (0/1)/(.375/.625) = 0. When we take the log of that value, it's negative infinity. So, the model is trying to send the coefficient out toward negative infinity. In Frequentist models, the telltale sign of this is coefficients very far away from zero (for the scale of the dependent variable) with very large standard errors. This is essentially what is happening here. In the Bayesian case, the prior may provide enough information for the model to converge.

In the frequentist world, these problems are often solved with regularization (e.g., Firth login, a penalized-likelihood solution). In the Bayesian setting, these results can be approximated. Normal priors with smaller variances would do something like an L-2 penalty while using Laplace priors would generate something akin to the L-1 penalty.

Best Answer

Related Solutions

Solved – Where is the correlation parameter in the linear mixed-effect model equation

Model Fit – Advice on Improving the Model Fit of a Mixed Model on Repeated Measurements

Edit - diagnosing the problem

Related Question