library(readxl)
resp <- read_excel("../data/resp_gee.xlsx")
resp$trtp<-factor(resp$trtp, levels=c('P', 'A'))
resp$avisitn<-factor(resp$avisitn)
resp$outcome<-factor(resp$outcome)
resp$usubjid<-factor(resp$usubjid)
resp$respord<-factor(resp$respord)Generalized Linear Mixed Models (GLMM)
INTRODUCTION
Generalized Linear Mixed Models (GLMM) method combines the characteristics of the Generalized Linear Model (GLM), with mixed models (such a repeated measures over time). It extends the GLM framework using link functions that relate the predictors to transformed outcome variable.
\[ E(Y)=\mu \]
\[ g(\mu) = X\beta + Zb, \qquad b \sim N(0, G) \]
Where:
n: number of observations, p: number of fixed effects, q: number of random effects (subjects).
Y: vector of observed response variable (n x 1)
g: Link function that transforms Y to the linear scale (eg: logit)
X: matrix for fixed effects (n x p), Z: matrix of random effects, G: covariance matrix of the random effects.
B: vector of fixed effects coefficients (p x 1)., b: vector of random effects.
Link Function:
Dichotomous response variable: probit (in case of rare events complementary log-log may be preferable).
Outcomes with more than two categories:
Ordinal variable: cumulative
Nominal variable: generalized logit
Random Effects
GLMM are conditional models and estimate subject-average effects, and the intra-subject correlation is modelled via random effects. Unlike GEE models, GLMM models allow individual-level inference.
Estimation Methods
Maximum likelihood, based on approximations:
Gauss Hermite Quadrature (GHQ): Integral split in a given number of points.
Laplace: A specific case of GHQ, using 1 point.
Penalized Likelihood can be used too, but it is known that in binary data it underestimates variance components and biased results.
EXAMPLE DATA
A SAS data of clinical trial data comparing two treatments for a respiratory disorder available in “Gee Model for Binary Data” [1] in the SAS/STAT Sample Program Library [1] is used to create these examples.
To uniquely identify subjects, a new variable USUBJID was created by concatenating SITE and ID. Variables TREATMENT, BASELINE, and VISIT were renamed to TRTP, BASE, and AVISITN.
Additionally, two variables were created using randomly generated values to simulate variables with more than two categories. One was an ordinal variable with values 1, 2, and 3; the other was a nominal variable with categories ‘liver’, ‘lung’, and ‘bone’. Data were created in SAS (See SAS section), and imported in R.
The variables OUTCOME, TRTP, RESPORD, USUBJID and AVISITN were converted to factors. Since the modeling functions use the first (alphabetically) level as the reference category, TRT levels are ordered as ‘P’ (placebo) and ‘A’ (active), to ensure that placebo is used as the reference category in the models:
PACKAGES
library(lme4)Loading required package: Matrixlibrary(GLMMadaptive)
Attaching package: 'GLMMadaptive'The following object is masked from 'package:lme4':
negative.binomiallibrary(merDeriv)Loading required package: nonnest2This is nonnest2 0.5-8.
nonnest2 has not been tested with all combinations of supported model classes.Loading required package: sandwichLoading required package: lavaanThis is lavaan 0.6-21
lavaan is FREE software! Please report any bugs.library(parameters)
library(glmmTMB)
Attaching package: 'glmmTMB'The following objects are masked from 'package:sandwich':
meatHC, sandwichlibrary(clubSandwich)Registered S3 methods overwritten by 'clubSandwich':
method from
bread.lmerMod merDeriv
bread.mlm sandwichlibrary(MASS)
Attaching package: 'MASS'The following object is masked from 'package:GLMMadaptive':
negative.binomiallibrary(emmeans)Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'library(ordinal)GLMM WITH GHQ
GLMM with GHQ approximation can be fitted using lme4::glmer andGLMMadaptative::mixed_modelR functions.
In GLMMs, intra‑subject (within‑cluster) correlation is captured through the inclusion of random effects. Both functions estimate this correlation by assuming normally distributed random effects with infinite degrees of freedom and, by default, model the random‑effects covariance matrix GGG using a Variance Components (VC) structure.
lme::glmer
The syntax to fit a GLMM using GHQ with lme4::glme is displayed below, the random effects are specified as 1|usubjid , where the number 1 denotes a random intercept (different baseline for each individual). By specifying nAQG=<n> the GHQ method is applied with the specified number of points (in this example n=5).
model <- lme4::glmer(outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link ="logit"),
nAGQ=5)
summary(model)Generalized linear mixed model fit by maximum likelihood (Adaptive
Gauss-Hermite Quadrature, nAGQ = 5) [glmerMod]
Family: binomial ( logit )
Formula: outcome ~ trtp + avisitn + trtp * avisitn + (1 | usubjid)
Data: resp
AIC BIC logLik -2*log(L) df.resid
487.2 524.0 -234.6 469.2 435
Scaled residuals:
Min 1Q Median 3Q Max
-2.1952 -0.3744 0.2062 0.4089 2.1081
Random effects:
Groups Name Variance Std.Dev.
usubjid (Intercept) 6.963 2.639
Number of obs: 444, groups: usubjid, 111
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.06177 0.52741 -0.117 0.9068
trtpA 1.57505 0.78352 2.010 0.0444 *
avisitn2 -0.87192 0.54763 -1.592 0.1113
avisitn3 -0.28729 0.53694 -0.535 0.5926
avisitn4 -0.43176 0.53860 -0.802 0.4228
trtpA:avisitn2 1.04386 0.80358 1.299 0.1939
trtpA:avisitn3 0.63613 0.79997 0.795 0.4265
trtpA:avisitn4 -0.22003 0.78660 -0.280 0.7797
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtpA avstn2 avstn3 avstn4 trtA:2 trtA:3
trtpA -0.676
avisitn2 -0.497 0.317
avisitn3 -0.505 0.334 0.496
avisitn4 -0.503 0.330 0.499 0.500
trtpA:vstn2 0.338 -0.482 -0.683 -0.339 -0.341
trtpA:vstn3 0.338 -0.488 -0.336 -0.672 -0.337 0.495
trtpA:vstn4 0.346 -0.518 -0.334 -0.340 -0.681 0.495 0.496GLMMadapive::mixed_model
The syntax using GLMMadaptative::mixed_model is similar, but the random effects are specified by using the random argument.
model <- GLMMadaptive::mixed_model(fixed = outcome ~ trtp + avisitn + trtp*avisitn,
random = ~1|usubjid,
data = resp,
family = binomial(link = "logit"),
nAGQ=5)
summary(model)
Call:
GLMMadaptive::mixed_model(fixed = outcome ~ trtp + avisitn +
trtp * avisitn, random = ~1 | usubjid, data = resp, family = binomial(link = "logit"),
nAGQ = 5)
Data Descriptives:
Number of Observations: 444
Number of Groups: 111
Model:
family: binomial
link: logit
Fit statistics:
log.Lik AIC BIC
-234.5954 487.1909 511.5767
Random effects covariance matrix:
StdDev
(Intercept) 2.680307
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -0.0760 0.5806 -0.1308 0.89590
trtpA 1.6388 0.8728 1.8776 0.06044
avisitn2 -0.8746 0.5475 -1.5974 0.11017
avisitn3 -0.2890 0.5374 -0.5377 0.59080
avisitn4 -0.4335 0.5390 -0.8043 0.42123
trtpA:avisitn2 1.0461 0.8031 1.3026 0.19272
trtpA:avisitn3 0.6375 0.8000 0.7969 0.42553
trtpA:avisitn4 -0.2188 0.7875 -0.2778 0.78118
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 5
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE Results produced by the mixed_model function show slight deviations when compared to glmer, suggesting differences in the model implementation or estimation.
Even with Gauss–Hermite quadrature in both, lme4::glmer and GLMMadaptive::mixed_model use different likelihood parameterizations, adaptive quadrature implementations, and optimization routines, so they don’t evaluate the exact same objective and can yield different estimates (especially for non-Gaussian outcomes or sparse data). In practice, mixed_model() tends to provide more accurate marginal likelihood integration (subject-specific adaptive scaling), while glmer() is faster but more approximate, which can bias variance components and slightly shift fixed effects [2].
GLMM WITH LAPLACE
Laplace is a particular GHQ where only one point is used. In R, lme4::glmer can also be used to compute Laplace approximation, but not GLMMadaptative::mixed_model. The glmmTMB::glmmTMBfunction is also a viable alternative.
lme::glmer
The syntax using lme4::glmer is the same as in GHQ, but without specifying the number of points:
model <- lme4::glmer(formula =outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link ="logit"))Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00926468 (tol = 0.002, component 1)
See ?lme4::convergence and ?lme4::troubleshooting.summary(model)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: outcome ~ trtp + avisitn + trtp * avisitn + (1 | usubjid)
Data: resp
AIC BIC logLik -2*log(L) df.resid
493.4 530.2 -237.7 475.4 435
Scaled residuals:
Min 1Q Median 3Q Max
-2.1946 -0.3899 0.2097 0.4216 2.0956
Random effects:
Groups Name Variance Std.Dev.
usubjid (Intercept) 6.433 2.536
Number of obs: 444, groups: usubjid, 111
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.08269 0.52044 -0.159 0.8738
trtpA 1.64646 0.78036 2.110 0.0349 *
avisitn2 -0.85805 0.54287 -1.581 0.1140
avisitn3 -0.28067 0.53169 -0.528 0.5976
avisitn4 -0.42263 0.53331 -0.792 0.4281
trtpA:avisitn2 1.02843 0.79405 1.295 0.1953
trtpA:avisitn3 0.62119 0.79067 0.786 0.4321
trtpA:avisitn4 -0.20289 0.77543 -0.262 0.7936
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtpA avstn2 avstn3 avstn4 trtA:2 trtA:3
trtpA -0.675
avisitn2 -0.497 0.310
avisitn3 -0.508 0.331 0.495
avisitn4 -0.506 0.327 0.497 0.500
trtpA:vstn2 0.339 -0.474 -0.685 -0.339 -0.340
trtpA:vstn3 0.340 -0.482 -0.336 -0.674 -0.338 0.493
trtpA:vstn4 0.350 -0.510 -0.336 -0.342 -0.685 0.496 0.496
optimizer (Nelder_Mead) convergence code: 0 (OK)
Model failed to converge with max|grad| = 0.00926468 (tol = 0.002, component 1)
See ?lme4::convergence and ?lme4::troubleshooting.This initial attempt using lme4::glmerresulted in a convergence warning, which can be addressed by switching the optimizer to “bobyqa” and extending the maximum number of iterations, as shown in the code below. These adjustments not only suppress the warning but provides closer results to SAS [3].
model <- lme4::glmer(formula =outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link ="logit"),
nAGQ=1,
control=lme4::glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=100000)))
summary(model)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: outcome ~ trtp + avisitn + trtp * avisitn + (1 | usubjid)
Data: resp
Control:
lme4::glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+05))
AIC BIC logLik -2*log(L) df.resid
493.4 530.2 -237.7 475.4 435
Scaled residuals:
Min 1Q Median 3Q Max
-2.1964 -0.3901 0.2097 0.4213 2.0951
Random effects:
Groups Name Variance Std.Dev.
usubjid (Intercept) 6.44 2.538
Number of obs: 444, groups: usubjid, 111
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.08475 0.52067 -0.163 0.8707
trtpA 1.65139 0.78093 2.115 0.0345 *
avisitn2 -0.85752 0.54293 -1.579 0.1142
avisitn3 -0.28186 0.53178 -0.530 0.5961
avisitn4 -0.42356 0.53339 -0.794 0.4271
trtpA:avisitn2 1.02468 0.79414 1.290 0.1969
trtpA:avisitn3 0.62207 0.79088 0.787 0.4315
trtpA:avisitn4 -0.20495 0.77559 -0.264 0.7916
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtpA avstn2 avstn3 avstn4 trtA:2 trtA:3
trtpA -0.675
avisitn2 -0.497 0.310
avisitn3 -0.508 0.331 0.495
avisitn4 -0.506 0.326 0.497 0.500
trtpA:vstn2 0.339 -0.473 -0.685 -0.339 -0.340
trtpA:vstn3 0.340 -0.482 -0.336 -0.673 -0.338 0.493
trtpA:vstn4 0.350 -0.510 -0.336 -0.342 -0.685 0.496 0.496model <- glmer(formula =outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link ="logit"),
nAGQ=1,
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=100000)))glmmTMB::glmmTMB
The glmmTMB::glmmTMB function is also a viable alternative for Laplace approximation, with similar syntax to lme4::glmer.
model<- glmmTMB::glmmTMB(outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link = "logit"))
summary(model) Family: binomial ( logit )
Formula: outcome ~ trtp + avisitn + trtp * avisitn + (1 | usubjid)
Data: resp
AIC BIC logLik -2*log(L) df.resid
493.4 530.2 -237.7 475.4 435
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
usubjid (Intercept) 6.442 2.538
Number of obs: 444, groups: usubjid, 111
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.08476 0.52073 -0.163 0.8707
trtpA 1.65178 0.78148 2.114 0.0345 *
avisitn2 -0.85758 0.54300 -1.579 0.1143
avisitn3 -0.28190 0.53181 -0.530 0.5961
avisitn4 -0.42359 0.53344 -0.794 0.4271
trtpA:avisitn2 1.02474 0.79432 1.290 0.1970
trtpA:avisitn3 0.62216 0.79103 0.786 0.4316
trtpA:avisitn4 -0.20503 0.77578 -0.264 0.7916
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Although both glmmTMB and lme4::glmer implement the Laplace approximation for fitting GLMMs, they rely on different numerical frameworks and optimization strategies. glmmTMB uses TMB with automatic differentiation and its own Laplace implementation, whereas glmer uses lme4’s custom penalized deviance formulation and derivative‑free optimizers. As a result, the approximated likelihoods, gradients, and optimization paths differ, leading to non‑identical parameter estimates even under the same model specification. [2]
PENALIZED QUASI-LIKELIHOOD (PQL)
The PQL approach uses linear approximations instead of likelihood, making it less accurate for binary outcomes compared to the GHQ or Laplace methods described above. PQL computation can be obtained using the glmmPQL function form the MASS package, using the random argument to specify the random factors.
model <- MASS::glmmPQL(outcome ~ trtp + avisitn + trtp*avisitn,
random=list(~1|usubjid),
data = resp,
family = binomial(link = "logit"))iteration 1iteration 2iteration 3iteration 4iteration 5iteration 6summary(model)Linear mixed-effects model fit by maximum likelihood
Data: resp
AIC BIC logLik
NA NA NA
Random effects:
Formula: ~1 | usubjid
(Intercept) Residual
StdDev: 2.240666 0.7414888
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: outcome ~ trtp + avisitn + trtp * avisitn
Value Std.Error DF t-value p-value
(Intercept) 0.0000080 0.4125606 327 0.0000194 1.0000
trtpA 1.2202691 0.6029058 109 2.0239797 0.0454
avisitn2 -0.8149934 0.3942084 327 -2.0674177 0.0395
avisitn3 -0.2682000 0.3879416 327 -0.6913412 0.4898
avisitn4 -0.4033782 0.3890430 327 -1.0368474 0.3006
trtpA:avisitn2 0.9759879 0.5795035 327 1.6841796 0.0931
trtpA:avisitn3 0.5950261 0.5780995 327 1.0292797 0.3041
trtpA:avisitn4 -0.1998455 0.5678652 327 -0.3519241 0.7251
Correlation:
(Intr) trtpA avstn2 avstn3 avstn4 trtA:2 trtA:3
trtpA -0.684
avisitn2 -0.464 0.317
avisitn3 -0.467 0.319 0.495
avisitn4 -0.466 0.319 0.496 0.498
trtpA:avisitn2 0.315 -0.469 -0.680 -0.337 -0.338
trtpA:avisitn3 0.313 -0.467 -0.332 -0.671 -0.334 0.492
trtpA:avisitn4 0.319 -0.485 -0.340 -0.341 -0.685 0.500 0.497
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.8440535 -0.4341097 0.2723992 0.4824627 2.8108150
Number of Observations: 444
Number of Groups: 111 Note: glmmPQL is widely recognized as less reliable for binary outcomes, more robust approaches such as Laplace or GHQ discussed in previous sections are generally preferred.
SANDWICH SE AND DEGREES OF FREEDOM (DDFF)
Previous results are computed using default outputs, so naive SE and infinite ddff are used. However, Li. P. and Redden, D.T. (2015) [4], suggested using the Between-Within denominator degrees of freedom approximation method when using GLMMs in randomized trials with binary outcomes and small sample size.
Additionally, FDA [5] advises “sponsors to consider using of robust SE method such as the Huber-White”sandwich” SE, particularly when the model does not include treatment by covariate interaction.
In R, the functions described above do not include these options. Function parameters::dof_betwithin can be used to use between-withing ddff, does it would not return exactly the same results as shown in Li and Redden 2015 [4], but similar results are obtained [6].
The Sandwich S.E. can be obtained from different functions too: merDeriv::sandwich function can be used for glmer objects, and clubSandwich::vovHC for glmmTMB objects (recent versions of the package R (>= 3.6.0) support the computation of Sandwich S.E. using the using clubSandwich::vcovHC function [7].
A function ‘new_mode’ is created below, where results obtained from a glmer model can be modified to use the sandwich SE and the between-within approximated ddff:
lme::glmer
new_model <- function(model, est_fix , se , df){
#Re-calculate (2-sided) p-values using estimated parameters and its SE
t <- est_fix /se
pvalue <-2*pt(q=abs(t), df=df, lower.tail=FALSE)
#Combine results in a data frame and add row names
new_model <- round(cbind(est_fix, se, df, t, pvalue), digits=4)
colnames(new_model) <- c("Estimate", "Std. Error", "df", "tvalue", "P-value")
rownames(new_model) <- rownames(summary(model)$coefficients)
new_model
}#Run the model
model <- lme4::glmer(outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link ="logit"),
nAGQ=5)
#Get parameter estimation
est<-lme4::fixef(model)
#Get Sandwich covariance matrix
library(merDeriv)
vcov<-sandwich(model)
#Get S.E. (the diagonal from the covariance matrix), remove last value as it corresponds to random effects
se_sw0<- sqrt(diag(vcov))
se_sw <-head(se_sw0, -1)#Re-calculate p-values using Sandwich SE and Inf ddff
new_model(model, est_fix=est, se=se_sw , df=Inf ) Estimate Std. Error df tvalue P-value
(Intercept) -0.0618 0.5378 Inf -0.1149 0.9086
trtpA 1.5751 0.7914 Inf 1.9902 0.0466
avisitn2 -0.8719 0.5391 Inf -1.6172 0.1058
avisitn3 -0.2873 0.6102 Inf -0.4708 0.6378
avisitn4 -0.4318 0.5164 Inf -0.8361 0.4031
trtpA:avisitn2 1.0439 0.8220 Inf 1.2699 0.2041
trtpA:avisitn3 0.6361 0.8589 Inf 0.7407 0.4589
trtpA:avisitn4 -0.2200 0.7939 Inf -0.2772 0.7817#Re-calculate p-values using Sandwich SE and between-within ddff
new_model(model, est_fix=est, se=se_sw , df=parameters::dof_betwithin(model)) Estimate Std. Error df tvalue P-value
(Intercept) -0.0618 0.5378 324 -0.1149 0.9086
trtpA 1.5751 0.7914 324 1.9902 0.0474
avisitn2 -0.8719 0.5391 324 -1.6172 0.1068
avisitn3 -0.2873 0.6102 324 -0.4708 0.6381
avisitn4 -0.4318 0.5164 324 -0.8361 0.4037
trtpA:avisitn2 1.0439 0.8220 324 1.2699 0.2050
trtpA:avisitn3 0.6361 0.8589 324 0.7407 0.4594
trtpA:avisitn4 -0.2200 0.7939 324 -0.2772 0.7818glmmTMB::glmmTMB
new_model <- function(model, est_fix , se , df){
#Re-calculate (2-sided) p-values using estimated parameters and its SE
t <- est_fix /se
pvalue <-2*pt(q=abs(t), df=df, lower.tail=FALSE)
#Combine results in a data frame and add row names
new_model <- round(cbind(est_fix, se, df, t, pvalue), digits=4)
colnames(new_model) <- c("Estimate", "Std. Error", "df", "tvalue", "P-value")
rownames(new_model) <- rownames(summary(model)$coefficients$cond)
new_model
}#Run the model
model<- glmmTMB::glmmTMB(outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link = "logit"))
#Get estimator
est<-glmmTMB::fixef(model)$cond
#Get Sandwich S.E.
se_sw<- sqrt(diag(vcovHC(model)))#Re-calculate p-values with Sandwich SE and Ininite df
new_model(model, est_fix=est, se=se_sw , df=Inf ) Estimate Std. Error df tvalue P-value
(Intercept) -0.0848 0.5452 Inf -0.1555 0.8765
trtpA 1.6518 0.8182 Inf 2.0188 0.0435
avisitn2 -0.8576 0.5302 Inf -1.6175 0.1058
avisitn3 -0.2819 0.5986 Inf -0.4709 0.6377
avisitn4 -0.4236 0.5069 Inf -0.8357 0.4033
trtpA:avisitn2 1.0247 0.8026 Inf 1.2767 0.2017
trtpA:avisitn3 0.6222 0.8397 Inf 0.7410 0.4587
trtpA:avisitn4 -0.2050 0.7727 Inf -0.2653 0.7907#Re-calculate p-values with Sandwich SE between-withing aprox. df
new_model(model, est_fix=est, se=se_sw , df=parameters::dof_betwithin(model)) Estimate Std. Error df tvalue P-value
(Intercept) -0.0848 0.5452 324 -0.1555 0.8766
trtpA 1.6518 0.8182 324 2.0188 0.0443
avisitn2 -0.8576 0.5302 324 -1.6175 0.1067
avisitn3 -0.2819 0.5986 324 -0.4709 0.6380
avisitn4 -0.4236 0.5069 324 -0.8357 0.4039
trtpA:avisitn2 1.0247 0.8026 324 1.2767 0.2026
trtpA:avisitn3 0.6222 0.8397 324 0.7410 0.4592
trtpA:avisitn4 -0.2050 0.7727 324 -0.2653 0.7909PREDICTED PROBABILITIES AND ODDS RATIO (OR)
Estimated probabilities and ORs, obtained by using the emmeans and contrast functions. In the example below, the probabilities and OR are estimated for GLMM using GHQ approximation.
To obtain robust (sandwich) standard errors, the emmeans function allows specifying a sandwich covariance estimator using vcov.method = sandwich(model), as illustrated in the example below. To instead use the model‑based (naive) standard errors, simply supply the model’s own covariance matrix by specifying vcov.method = vcov(model).
#Run the model
model <- lme4::glmer(outcome ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp,
family = binomial(link ="logit"),
nAGQ=5)library(emmeans)
#Get predicted probabilities for each treatment, using the Sandwich matrix covariance
prob <- emmeans(model, ~ trtp*avisitn, data=resp, vcov.method=sandwich(model), type='response')
prob trtp avisitn prob SE df asymp.LCL asymp.UCL
P 1 0.485 0.1320 Inf 0.251 0.726
A 1 0.820 0.0854 Inf 0.594 0.934
P 2 0.282 0.1090 Inf 0.120 0.531
A 2 0.844 0.0772 Inf 0.631 0.944
P 3 0.414 0.1280 Inf 0.200 0.666
A 3 0.866 0.0691 Inf 0.668 0.954
P 4 0.379 0.1250 Inf 0.177 0.634
A 4 0.703 0.1150 Inf 0.445 0.875
Confidence level used: 0.95
Intervals are back-transformed from the logit scale #Get differences between treatments by visit (option "revpairwise" is used to compare A vs P)
diffs_mean<-emmeans(model, ~ trtp*avisitn, data=resp, vcov.method=sandwich(model))
diffs <- contrast(diffs_mean,"revpairwise", simple="trtp")
diffs <- as.data.frame(diffs)
#Calculate CI (alpha=0.05)
alpha<-0.05
z_crit <- qnorm(1 - alpha / 2)
diffs$low<-diffs$estimate - (z_crit*diffs$SE)
diffs$upp<-diffs$estimate + (z_crit*diffs$SE)
#Get OR applying exponential transformation;
or<-exp(diffs$estimate)
or_low<-exp(diffs$low)
or_upp<-exp(diffs$up)
#Get two-sided p-value
z <- diffs$estimate/diffs$SE
pvalue <- 2 * (1 - pnorm(z))
#Create a dataset with all the results
OR<-as.data.frame(cbind(diffs$avisitn, or,or_low, or_upp, z, round(pvalue, digits=4)))
colnames(OR)<-c('avisit', 'OR', 'lower.CL', 'upper.CL', 'Z', 'p-value')
OR avisit OR lower.CL upper.CL Z p-value
1 1 4.830986 1.0401682 22.43716 2.010228 0.0444
2 2 13.720725 2.8152984 66.86975 3.240837 0.0012
3 3 9.126521 1.8973419 43.90004 2.759125 0.0058
4 4 3.876834 0.8558308 17.56170 1.757983 0.0788OUTCOME WITH MORE THAN TWO CATEGORIES
Although less common than binary outcomes, endpoints with more than two categories may be the outcome of interest, which can be either ordinal or nominal. In those case, the multinomial distribution is used selecting the appropriated link function depending on the type of response variable.
Ordinal variable
The function ordinal::clmm can be used to fit a GLMM model when the outcome is an ordinal variable.
model<- ordinal::clmm(respord ~ trtp + avisitn + trtp*avisitn + (1 | usubjid),
data = resp)
summary(model)Cumulative Link Mixed Model fitted with the Laplace approximation
formula: respord ~ trtp + avisitn + trtp * avisitn + (1 | usubjid)
data: resp
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 444 -482.20 984.40 666(1357) 6.84e-04 1.7e+02
Random effects:
Groups Name Variance Std.Dev.
usubjid (Intercept) 0.1511 0.3887
Number of groups: usubjid 111
Coefficients:
Estimate Std. Error z value Pr(>|z|)
trtpA -0.43273 0.36062 -1.200 0.230
avisitn2 0.01155 0.35236 0.033 0.974
avisitn3 -0.10700 0.35079 -0.305 0.760
avisitn4 0.17114 0.34838 0.491 0.623
trtpA:avisitn2 0.26152 0.50072 0.522 0.601
trtpA:avisitn3 0.57356 0.50391 1.138 0.255
trtpA:avisitn4 0.31098 0.49794 0.625 0.532
Threshold coefficients:
Estimate Std. Error z value
1|2 -0.5630 0.2573 -2.188
2|3 0.6573 0.2589 2.539Nominal variable
No R functions have been identified for handling multinomial distributions for a nominal variable in a frequentist framework.
REFERENCES
[1] SAS Institute Inc.. SAS Help Center. The GEE procedure.
[2] Conceptual explanations were assisted using Microsoft Copilot (M365 Copilot, GPT‑5‑based model).
[3] Stack Overflow [Internet]. 2008 [Last visited: 2025 Sep 30]
[6] Documentation of package parameters. dof_betwithin
[7] Brooks, M. E., et al. (2025). glmmTMB: Generalized Linear Mixed Models using Template Model Builder (Version 1.1.12) [R package manual]. The Comprehensive R Archive Network (CRAN). https://cran.r-project.org/web/packages/glmmTMB/glmmTMB.pdf