So far we studied linear mixed models (LMM), where the responses are continuous and normally distributed. It is natural to extend LMM to handle nonnormally distributed responses.
Recall that, in GLM, the distribution of \(Y_i\) is from the exponential family of distributions of form \[
f(y_i \mid \theta_i, \phi) = \exp \left[ \frac{y \theta_i - b(\theta_i)}{a(\phi)} + c(y_i, \phi) \right].
\] If we use the canonical link, then \[
\theta_i = g(\mu_i) = \eta_i,
\] where \(\mu_i = \mathbb{E} Y_i\). Now let \[
\eta_i = \mathbf{x}_i^T \boldsymbol{\beta} + \mathbf{z}_i^T \boldsymbol{\gamma}_i,
\] where \(\boldsymbol{\beta}\) is the fixed effects and \(\boldsymbol{\gamma}\) is the random effects with density \(h(\gamma \mid \boldsymbol{\Sigma})\). (Typically we assume multivariate normal with mean 0 and covariance \(\boldsymbol{\Sigma}\).) Then the likelihood is \[
L(\boldsymbol{\beta}, \phi, \boldsymbol{\Sigma}) = \prod_{i=1}^n \int f(y_i \mid \boldsymbol{\beta}, \phi, \boldsymbol{\gamma}) h(\boldsymbol{\gamma} \mid \boldsymbol{\Sigma}) \, d \boldsymbol{\gamma}
\] and the log-likelihood is \[
\ell(\boldsymbol{\beta}, \phi, \boldsymbol{\Sigma}) = \sum_{i=1}^n \log \int f(y_i \mid \boldsymbol{\beta}, \phi, \boldsymbol{\gamma}) h(\boldsymbol{\gamma} \mid \boldsymbol{\Sigma}) \, d \boldsymbol{\gamma}.
\] This model, combining GLM and mixed effects model, is called the generalized linear mixed effects model (GLMM).
2 Overview of estimation and inference methods
Estimation and inference of GLMM have been active research. We give an overview of a few commonly used methods.
Maximum likelihood estimate (MLE) is obtained by maximizing the log-likelihood function. However, each iteration of the optimization algorithm needs to evaluate numerical integration, which quickly becomes infeasible when the dimension of random effects \(\boldsymbol{\gamma}\) is large.
Gauss-Hermite quadrature can be used for numerical integration. It approximates the integral by a weighted sum of integrands at different knots. The more knots, the more accurate the approximation, but at higher computational expense.
Laplace approximation is a special case of Gauss-Hermite quadrature with only one knot. It is less accurate than Guass-Hermite quadrature but computationally cheaper.
Bayesian method is another approach for estimation and inference of GLMM. It’s not covered in this course due to time constraint. Interested students can study ELMR Chapter 12.
Penalized quasi-likelihood (PQL). In the IRWLS (iteratively reweighted least squares) procedure for estimating the GLM, each iteration uses the pseudo-response or working response \[
\mathbf{z}^{(t)} = \boldsymbol{\eta}^{(t)} + (\mathbf{W}^{(t)})^{-1} (\mathbf{y} - \widehat{\boldsymbol{\mu}}^{(t)}).
\] It can be adapted to the mixed effects model. PQL has computational advantage but generally considered less accurate than other approaches.
Generalized estimation equations (GEE). To be discussed later.
3 Binary response example
ctsib data in EMLR studies the balance of 40 individuals.
Response: stability (1 = stable, 0 = instable).
Predictors: sex (binary), age (continuous), height (continuous), surface (binary, norm or foam), vision (factor, open or closed or dome).
Each subject is tested twice on a combination of surface and vision. 12 measurements per subject.
# A tibble: 480 × 9
Subject Sex Age Height Weight Surface Vision CTSIB stable
<int> <fct> <int> <dbl> <dbl> <fct> <fct> <int> <dbl>
1 1 male 22 176 68.2 norm open 1 1
2 1 male 22 176 68.2 norm open 1 1
3 1 male 22 176 68.2 norm closed 2 0
4 1 male 22 176 68.2 norm closed 2 0
5 1 male 22 176 68.2 norm dome 1 1
6 1 male 22 176 68.2 norm dome 2 0
7 1 male 22 176 68.2 foam open 2 0
8 1 male 22 176 68.2 foam open 2 0
9 1 male 22 176 68.2 foam closed 2 0
10 1 male 22 176 68.2 foam closed 2 0
11 1 male 22 176 68.2 foam dome 2 0
12 1 male 22 176 68.2 foam dome 2 0
13 2 male 22 181 67.6 norm open 1 1
14 2 male 22 181 67.6 norm open 1 1
15 2 male 22 181 67.6 norm closed 2 0
16 2 male 22 181 67.6 norm closed 2 0
17 2 male 22 181 67.6 norm dome 2 0
18 2 male 22 181 67.6 norm dome 2 0
19 2 male 22 181 67.6 foam open 2 0
20 2 male 22 181 67.6 foam open 2 0
21 2 male 22 181 67.6 foam closed 3 0
22 2 male 22 181 67.6 foam closed 3 0
23 2 male 22 181 67.6 foam dome 3 0
24 2 male 22 181 67.6 foam dome 3 0
# ℹ 456 more rows
summary(ctsib)
Subject Sex Age Height Weight
Min. : 1.00 female:240 Min. :18.00 Min. :142.0 Min. : 44.20
1st Qu.:10.75 male :240 1st Qu.:21.75 1st Qu.:166.8 1st Qu.: 60.65
Median :20.50 Median :25.50 Median :173.0 Median : 68.00
Mean :20.50 Mean :26.80 Mean :172.1 Mean : 71.14
3rd Qu.:30.25 3rd Qu.:33.00 3rd Qu.:180.2 3rd Qu.: 83.50
Max. :40.00 Max. :38.00 Max. :190.0 Max. :102.40
Surface Vision CTSIB stable
foam:240 closed:160 Min. :1.000 Min. :0.0000
norm:240 dome :160 1st Qu.:2.000 1st Qu.:0.0000
open :160 Median :2.000 Median :0.0000
Mean :1.919 Mean :0.2375
3rd Qu.:2.000 3rd Qu.:0.0000
Max. :4.000 Max. :1.0000
Graphical summary.
subsum <- ctsib %>%group_by(Subject) %>%summarise(Height = Height[1], Weight = Weight[1], stable =mean(stable), Age = Age[1], Sex = Sex[1])subsum %>%ggplot() +geom_point(mapping =aes(x = Height, y = stable))
subsum %>%ggplot() +geom_point(mapping =aes(x = Weight, y = stable))
subsum %>%ggplot() +geom_point(mapping =aes(x = Age, y = stable))
subsum %>%ggplot() +geom_boxplot(mapping =aes(x = Sex, y = stable))
seizures id treat expind timeadj
Min. : 0.00 Min. : 1 Min. :0.0000 Min. :0.0 Min. :2.0
1st Qu.: 3.00 1st Qu.:15 1st Qu.:0.0000 1st Qu.:1.0 1st Qu.:2.0
Median : 6.00 Median :30 Median :1.0000 Median :1.0 Median :2.0
Mean : 12.86 Mean :30 Mean :0.5254 Mean :0.8 Mean :3.2
3rd Qu.: 14.50 3rd Qu.:45 3rd Qu.:1.0000 3rd Qu.:1.0 3rd Qu.:2.0
Max. :151.00 Max. :59 Max. :1.0000 Max. :1.0 Max. :8.0
age period drug phase
Min. :18.00 Min. :0 placebo :140 baseline : 59
1st Qu.:22.00 1st Qu.:1 treatment:155 experiment:236
Median :28.00 Median :2
Mean :28.85 Mean :2
3rd Qu.:35.00 3rd Qu.:3
Max. :57.00 Max. :4
Numerical summary. Is there any difference in the number of seizures per week in the placebo and treatment groups?
Linear mixed-effects model fit by maximum likelihood
Data: epilo
AIC BIC logLik
NA NA NA
Random effects:
Formula: ~1 | id
(Intercept) Residual
StdDev: 0.6820012 1.605385
Variance function:
Structure: fixed weights
Formula: ~invwt
Fixed effects: seizures ~ offset(log(timeadj)) + expind + treat + I(expind * treat)
Value Std.Error DF t-value p-value
(Intercept) 1.0761832 0.09990514 230 10.772050 0.0000
expind 0.1125119 0.07412152 230 1.517939 0.1304
I(expind * treat) -0.3037615 0.10819095 230 -2.807642 0.0054
Correlation:
(Intr) expind
expind -0.198
I(expind * treat) -0.014 -0.656
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2934834 -0.5649468 -0.1492931 0.3224895 6.3123337
Number of Observations: 290
Number of Groups: 58
5 GEE (generalized estimation equation)
The quasi-likelihood (quasi-binomial, quasi-Poisson, etc) approach can be generalized to handle correlated, nonnormal responses.
Let \(\mathbf{Y}_i\) be the response vector of \(i\)-th individual. \[\begin{eqnarray*}
\mathbb{E} \mathbf{Y}_i &=& \boldsymbol{\mu}_i \\
g(\boldsymbol{\mu}_i) &=& \mathbf{X}_i \boldsymbol{\beta} \\
\mathbf{V}_i &=& \text{Var}(\mathbf{Y}_i) = \phi \mathbf{A}_i^{1/2} \mathbf{R}_i(\boldsymbol{\alpha}) \mathbf{A}_i^{1/2},
\end{eqnarray*}\] where \(\mathbf{A}_i = \text{diag}(a(\boldsymbol{\mu}))\) captures the individual variances and \(\mathbf{R}_i(\boldsymbol{\alpha})\) is a working correlation matrix.
Given estimates of \(\phi\) and \(\boldsymbol{\alpha}\), we solve the estimation equation\[
\sum_i [D_{\boldsymbol{\beta}} \boldsymbol{\mu}_i(\boldsymbol{\beta})]^T \mathbf{V}_i^{-1}(\mathbf{Y}_i - \boldsymbol{\mu}_i) = \mathbf{0}.
\]
Let’s revisit the ctsib (stability) data set. We use the exchangeable correlation, or equivalently, compound symmetry. Standard errors and Wald test are constructed using the sandwich estimators. \(\boldsymbol{\beta}\) estimates in GEE are about half the size of those from GLMM. The scale.fix = TRUE instructs the function geeglm not estimate the dispersion parameter \(\phi\), since for binary response the dispersion is always 1.
