HSCI 410 · Lesson 5

Modelling Dependent Data

Exploratory Data Analysis For Epidemiology

Learning objectives for this lesson:

  • Recognise clustered and repeated-measures data, identify the cluster variable, and explain why ignoring the clustering makes standard errors too small, especially for cluster-level predictors.
  • Estimate the intracluster correlation coefficient (ICC) from a mixed model with no predictors, and use the design effect and the effective sample size to describe how much information the clustering removes.
  • Fit and interpret a linear mixed model with a random intercept in R using lmer(), reading its fixed effects, its clinic and residual variances, and the ICC after adjustment.
  • Check a linear mixed model with residual plots, a Q-Q plot of the cluster effects, the number of clusters and a test for a singular fit.
  • Fit a logistic mixed model (GLMM) with glmer() and a GEE model with geeglm() for a binary outcome, and explain the difference between cluster-specific and population-averaged odds ratios.
  • Arrange repeated measures in long format, fit a mixed model and a GEE model with a group-by-time interaction, and explain how missed visits affect each analysis.

This course was developed by Dr. Kiffer G. Card, Faculty of Health Sciences, Simon Fraser University based on Dohoo, I. R., Martin, S. W., & Stryhn, H. (2012). Methods in Epidemiologic Research. VER Inc.

Reference

Glossary: Key Terms, People & Concepts

📚 Reference page, available throughout the lesson

This glossary collects the key concepts, methods and people in this lesson. It can be used as a reference while working through the material or as a review before assessments. Typing in the search box filters the entries.

Key Concepts & Ideas
Independence Independence is the assumption that one observation tells us nothing about another. It fails when observations share a cluster, such as patients in the same clinic or visits by the same person.
Dependent Data Dependent data are data in which some observations are more alike than others because they share a cluster. Clustered data and repeated measures are the two common forms.
Cluster A cluster is the unit that groups observations together, such as a clinic, a school, a neighbourhood or a person who is measured several times. The cluster variable identifies which cluster each observation belongs to.
Clustered Data Clustered data arise when people are grouped inside larger units, such as patients in clinics or students in schools. Such data are also called multilevel or hierarchical data.
Repeated Measures Repeated measures are measurements of the same outcome on the same person at several times. Each person forms a cluster of visits, and studies of this kind are also called longitudinal studies.
Cluster-Level Predictor A cluster-level predictor takes the same value for everyone in a cluster, such as whether a clinic is urban. Its standard error is the most affected when clustering is ignored.
Person-Level Predictor A person-level predictor can differ between people in the same cluster, such as age or smoking. Its standard error is usually less affected by clustering.
Long and Wide Format In long format a repeated-measures dataset has one row per visit, so each person appears on several rows. In wide format it has one row per person with a column for each visit. Mixed models and GEE in R need long format.
Missing at Random (MAR) Data are missing at random when the chance that a value is missing depends only on information the model uses, such as a person's study arm or earlier measurements. Mixed models that use every attended visit are valid under this assumption, while GEE needs the stronger condition that missingness depends only on the predictors in the model.
Complete-Case Analysis A complete-case analysis keeps only people with no missing values. In repeated-measures data it can discard many people and can bias the results when the people who miss visits differ from those who do not.
Methods & Statistical Concepts
Intracluster Correlation Coefficient (ICC) The intracluster correlation coefficient is the share of the total variation in an outcome that lies between clusters. It equals the between-cluster variance divided by the sum of the between-cluster and within-cluster variances, and it ranges from 0 (no clustering) to 1. HSCI 410 Lesson 7 uses the intraclass correlation coefficient, the same kind of variance ratio, for a different purpose: the reliability of a score.
Design Effect The design effect is the factor by which clustering inflates the variance of an estimate compared with a simple random sample of the same size. For clusters of average size m it equals 1 + (m − 1) × ICC.
Effective Sample Size The effective sample size is the number of independent observations that would carry the same information as the clustered sample. It equals the sample size divided by the design effect.
Mixed Model A mixed model is a regression model that contains both fixed effects for the predictors and random effects for the clusters. It is fitted in R with lmer() for a continuous outcome and glmer() for a binary outcome.
Fixed Effect A fixed effect is a coefficient for a predictor in a mixed model, such as the effect of smoking. It is read in the same way as a coefficient of an ordinary regression model.
Random Effect A random effect is a shift for each cluster, up or down from the overall average, that is assumed to come from a normal distribution. The model estimates the variance of these shifts, and ranef() prints the estimated shift for each cluster.
Random Intercept A random intercept gives each cluster its own baseline level of the outcome. In lmer() and glmer() it is written (1 | cluster).
Random Slope A random slope lets the effect of a predictor, such as time, differ from cluster to cluster. It is written (1 + time | cluster), which is the same as (time | cluster), and needs more data than a random intercept, so it often fails to converge in small studies.
Variance Components The variance components of a mixed model are the between-cluster variance and the residual (within-cluster) variance. They are printed under Random effects in the model summary and by VarCorr().
Linear Mixed Model (LMM) A linear mixed model is a linear regression with random effects for the clusters. It is fitted with lmer() from the lme4 package, and loading lmerTest adds p-values for the fixed effects.
Singular Fit A singular fit occurs when a mixed model estimates a variance of zero (or a correlation of exactly ±1) for a random effect. R reports "boundary (singular) fit", and isSingular() returns TRUE. It usually means the random-effects part is more complex than the data can support.
Generalized Linear Mixed Model (GLMM) A generalized linear mixed model is a generalized linear model, such as logistic regression, with random effects for the clusters. It is fitted with glmer(), and its exponentiated fixed effects are cluster-specific odds ratios.
Generalized Estimating Equations (GEE) GEE fits a regression model, such as logistic regression, without cluster effects and calculates standard errors that allow for the correlation within clusters. A working correlation sets how the data are weighted. It is fitted with geeglm() from the geepack package, and its odds ratios are population-averaged.
Working Correlation The working correlation is the pattern of correlation within clusters that a GEE model assumes. The exchangeable structure assumes the same correlation between any two members of a cluster, and an autoregressive structure, AR(1), assumes that measurements closer in time are more alike.
Sandwich Standard Error A sandwich standard error is the standard error that GEE reports. It is calculated by adding up the residuals within each cluster and measuring how much these cluster totals vary from one cluster to another, so it needs a reasonable number of clusters. It remains valid when the working correlation is not exactly right, provided there are enough clusters.
Cluster-Specific (Conditional) Effect A cluster-specific effect compares two people in the same cluster, so the cluster's shift cancels out. For a cluster-level predictor, such as an urban clinic, it compares two clusters with the same baseline. The fixed effects of a mixed model are cluster-specific.
Population-Averaged (Marginal) Effect A population-averaged effect compares the average outcome of two groups across all clusters. GEE estimates population-averaged effects, and for a binary outcome they are usually closer to 1 than the cluster-specific odds ratios of a GLMM.
Group-by-Time Interaction A group-by-time interaction is the term arm:visit, which R includes together with arm and visit when the formula contains arm * visit. It tests whether the outcome changes at a different rate in two groups. In a trial with repeated measures it usually answers the main question.
Key People
Ronald A. Fisher (1890–1962) Fisher was a British statistician and geneticist who described the intraclass correlation, the statistic on which the intracluster correlation coefficient is based, in his 1925 book Statistical Methods for Research Workers.
Leslie Kish (1910–2000) Kish was a Hungarian-born American survey statistician who introduced the design effect in his 1965 book Survey Sampling.
Nan Laird and James Ware Laird and Ware were biostatisticians at Harvard University whose 1982 paper set out the random-effects (mixed) model for longitudinal data in the form that is still used today.
Kung-Yee Liang and Scott Zeger Liang and Zeger were biostatisticians at Johns Hopkins University when they introduced generalized estimating equations for longitudinal data in 1986.
No matching entries. Try a different search term.
Section 1 of 4

Recognising Dependent Data

⏱ Estimated time: 40 minutes
Lesson 5 · Section 1

Recognising Dependent Data

When observations share a clinic, a household or a person, they are no longer independent of one another.

Two forms

Clustered data and repeated measures

Clustered data

People are grouped inside larger units, such as patients in clinics, students in schools or residents in neighbourhoods.

Repeated measures

The same person is measured several times, so each person's measurements form a cluster.

In both forms, observations from the same cluster tend to be more alike than observations from different clusters.

Why it matters

Ignoring dependence makes a model overconfident

  • The model counts every observation as a separate piece of information.
  • Its standard errors are too small, its confidence intervals too narrow and its p-values too small.
  • The problem is largest for predictors that describe the cluster, such as an urban or rural clinic.
Running example

966 patients in 30 clinics

Dot plot of systolic blood pressure for 966 patients, one column per clinic, sorted by clinic average. Clinic averages range from about 102 to 125 mmHg. Urban clinics are teal and rural clinics are orange.
Each column is one clinic and each dot is one patient; the diamonds mark the clinic averages.
Measuring clustering

The intracluster correlation coefficient (ICC)

Intracluster correlation coefficient
\[ \text{ICC} = \frac{\color{#C2410C}{\sigma^2_{\text{between}}}}{\color{#C2410C}{\sigma^2_{\text{between}}} + \color{#0B7B6B}{\sigma^2_{\text{within}}}} \]
σ²between variance between cluster averages σ²within variance among people in the same cluster

The ICC is also the correlation between the outcomes of two people in the same cluster.

What the ICC looks like

Three levels of clustering

Three simulated panels of ten clusters each, with ICC 0, 0.2 and 0.8. With ICC 0 the cluster averages are all near zero; with ICC 0.8 the clusters are tight and far apart.
The panels show simulated clusters with an ICC of 0, 0.2 and 0.8, and the diamonds mark the cluster averages.
ICC in the clinic data

About 17.5% of the variation lies between clinics

22.1
variance between clinics
104.3
variance within clinics
0.175
ICC = 22.1 ÷ (22.1 + 104.3)

ICCs for clinics, schools and neighbourhoods are often between 0.01 and 0.1, so 0.175 indicates strong clustering.

Design effect

From 966 patients to about 150 independent people

Design effect and effective sample size
\[ \text{DEFF} = 1 + (\color{#0B7B6B}{\bar m} - 1) \times \color{#C2410C}{\text{ICC}} \qquad n_{\text{eff}} = \frac{n}{\text{DEFF}} \]
m̄ average cluster size (32.2) ICC intracluster correlation coefficient (0.175)
6.45
design effect ≈ 1 + 31.2 × 0.175
150
effective sample size = 966 ÷ 6.45
What goes wrong

The same difference, with and without allowing for clinics

Two estimates of the urban minus rural difference in systolic blood pressure. Ignoring clinics: 1.49 with 95% CI 0.04 to 2.95. Allowing for clinics: 1.67 with 95% CI minus 2.17 to 5.51.
Ignoring clinics gives p = 0.044; allowing for clinics gives p = 0.40.
Which predictors are affected?

Cluster-level predictors lose the most precision

Cluster-level predictors

These take the same value for everyone in a cluster, such as an urban clinic or a school assigned to one trial arm. Their standard errors grow the most.

Person-level predictors

These vary within each cluster, such as age or smoking. Their standard errors change little and can even shrink.

Two families of methods

Mixed models and GEE

MethodWhat it doesR function
Linear mixed modelGives each cluster its own baseline (a random intercept) for a continuous outcomelme4::lmer()
Generalized linear mixed model (GLMM)Adds a random intercept to a logistic or Poisson modellme4::glmer()
Generalized estimating equations (GEE)Fits a regression for the whole population, with standard errors that allow for clusteringgeepack::geeglm()
Carry forward

What to take into the next section

  • Observations are dependent when they share a cluster, such as a clinic or a person measured several times.
  • The ICC measures how alike observations in a cluster are, and the design effect shows the information lost.
  • Ignoring dependence makes a model overconfident, and mixed models and GEE are the standard remedies.

Introduction and Overview

Every regression model in Lessons 3 and 4 assumed that the observations are independent, so that one person's outcome carries no information about another person's outcome. This lesson deals with data in which that assumption fails. Observations are dependent when they share something, such as patients attending the same clinic or the same person measured at several visits. This first section explains how to recognise dependent data, how to measure the strength of the dependence with the intracluster correlation coefficient (ICC) and the design effect, and what goes wrong when the dependence is ignored. The later sections fit models that allow for it.

Learning Objectives

  • Recognise clustered data and repeated measures, and identify the variable that defines the clusters.
  • Explain why ignoring dependence makes standard errors, confidence intervals and p-values too small.
  • Calculate and interpret the intracluster correlation coefficient (ICC), the design effect and the effective sample size.
  • Distinguish cluster-level from person-level predictors and explain why cluster-level predictors are affected most.
  • Compare an ordinary regression with a mixed model in R and describe how the standard error changes.

What Are Dependent Data?

Data are dependent when observations can be grouped into clusters whose members are more alike than members of different clusters. Two forms are common in public health. In clustered data, people are grouped inside larger units: patients in clinics, students in schools, residents in neighbourhoods, or workers in workplaces. People in the same cluster share a setting, staff, policies and local conditions, so their outcomes tend to be similar. In repeated measures (also called longitudinal data), the same person is measured more than once, and each person's measurements form a cluster. Clusters can also be nested inside one another, such as patients within doctors within clinics, which this lesson mentions only briefly.

Table 5.1 gives four examples of studies that produce dependent data, with the cluster, the observation and the reason why observations in the same cluster are alike.

Table 5.1. Examples of clustered data and repeated measures in public health studies.

StudyClusterObservationWhy observations in a cluster are alike
A survey of patients at 30 clinicsClinicPatientPatients share doctors, practices and a local population.
A school-based physical activity programSchoolStudentStudents share teachers, facilities and a neighbourhood.
A household survey of food securityHouseholdHousehold memberMembers share income, meals and housing.
A trial with visits at 0, 6, 12 and 18 monthsPersonVisitEach person brings their own genes, habits and baseline health to every visit.

The three flip cards below define the cluster, independence and the effective sample size, terms that the rest of the lesson uses repeatedly.

ClusterClick to explore
IndependenceClick to explore
Effective Sample SizeClick to explore

Dependent data therefore arise whenever a study samples people in groups or measures the same people repeatedly, and every method in this lesson needs the variable that identifies the clusters. The next part explains why this grouping matters for a regression model.

Why Dependence Matters

A regression model that assumes independence treats every observation as a separate piece of information. When observations in a cluster are alike, part of what each one says has already been said by the others, so the model overstates how much information it has. The point estimates are often similar, but the standard errors are too small, the confidence intervals too narrow and the p-values too small. In practice, this means that an analysis that ignores clustering can report a statistically significant association that the data do not support (Dohoo, Martin & Stryhn, 2012, Chapter 20).

The clinic data introduced next show how large this effect can be.

The Clinic Data

The first three sections of this lesson use phaa_clinics.csv, a simulated dataset of 966 patients attending 30 primary care clinics, with between 18 and 45 patients per clinic (32.2 on average). It records each patient's age, gender (female = 1 for women), smoking status, BMI, systolic blood pressure (sbp) and whether they were referred to a specialist (referred = 1). Two variables describe the clinic itself: clinic_urban (20 urban and 10 rural clinics) and clinic_size.

Figure 5.1 plots the systolic blood pressure of every patient by clinic, so that the differences between clinic averages can be seen before any model is fitted.

Dot plot of systolic blood pressure for 966 patients in 30 clinics, sorted by clinic average. Each column is a clinic; diamonds show clinic averages from about 102 to 125 mmHg. Urban clinics are teal and rural clinics orange.
Figure 5.1. Each column shows the patients of one clinic, sorted from the lowest clinic average (left) to the highest (right), and the diamonds mark the clinic averages. The averages range from about 102 to 125 mmHg, so a patient's blood pressure depends partly on the clinic they attend.

Figure 5.1 shows that a patient’s blood pressure depends partly on the clinic the patient attends. The intracluster correlation coefficient, introduced next, measures how much of the variation in blood pressure lies between clinics.

Measuring Clustering: The Intracluster Correlation Coefficient

The intracluster correlation coefficient (ICC) measures how alike observations in the same cluster are. It divides the total variance of the outcome into a part between clusters (how much the cluster averages differ) and a part within clusters (how much people in the same cluster differ from one another), and it is the share of the total that lies between clusters. Equation 5.1 expresses this share as a ratio of the two variances.

Intracluster correlation coefficient
\[ \text{ICC} = \frac{\color{#C2410C}{\sigma^2_{\text{between}}}}{\color{#C2410C}{\sigma^2_{\text{between}}} + \color{#0B7B6B}{\sigma^2_{\text{within}}}} \]Eq 5.1
The ICC is the variance between clusters divided by the total of the variance between clusters and the variance within clusters. It runs from 0 (no clustering) to 1 (everyone in a cluster has the same value).

The ICC can also be read as the correlation between the outcomes of two people chosen from the same cluster. Lesson 7 uses a statistic of the same form, the intraclass correlation coefficient, for a different application: measuring the reliability of a score across occasions or raters. Figure 5.2 shows simulated clusters with three values of the ICC.

Three simulated panels, each with ten clusters of fifteen observations. With ICC 0 the cluster averages are all close to zero. With ICC 0.2 they spread a little. With ICC 0.8 each cluster is tight and the clusters are far apart.
Figure 5.2. The panels show simulated data with ICCs of 0, 0.2 and 0.8, and the diamonds mark the cluster averages. As the ICC rises, the clusters become tighter and their averages move further apart.

Worked Example 5.1 applies Equation 5.1 to the clinic data, using the two variances that Activity 5.1 estimates at the end of this section.

Worked Example 5.1: The ICC for Blood Pressure

A mixed model with no predictors (Activity 5.1 below) estimates a between-clinic variance of 22.10 and a within-clinic variance of 104.34. The ICC is 22.10 ÷ (22.10 + 104.34) = 22.10 ÷ 126.44 = 0.175. About 17.5% of the variation in systolic blood pressure lies between clinics, and the blood pressures of two patients from the same clinic have a correlation of about 0.175.

ICCs for clinics, schools and neighbourhoods in health research are often between 0.01 and 0.1, so a value of 0.175 indicates strong clustering. Small ICCs still matter when clusters are large, as the design effect shows.

The ICC of 0.175 measures the strength of the clustering in the clinic data. The next part converts it into a loss of information with the design effect and the effective sample size.

The Design Effect and the Effective Sample Size

The design effect (DEFF) expresses the ICC as a loss of information. It depends on both the ICC and the average cluster size. Equation 5.2 gives the design effect and the effective sample size that follows from it.

Design effect and effective sample size
\[ \text{DEFF} = 1 + (\color{#0B7B6B}{\bar m} - 1) \times \color{#C2410C}{\text{ICC}} \qquad n_{\text{eff}} = \frac{n}{\text{DEFF}} \]Eq 5.2
The design effect is 1 plus the average cluster size minus 1, multiplied by the ICC. Dividing the sample size by the design effect gives the effective sample size.

For the clinic data, DEFF = 1 + (32.2 − 1) × 0.1748 = 6.45 (using the unrounded ICC), and the effective sample size is 966 ÷ 6.45 = 150. For questions about clinic characteristics, the 966 clustered patients therefore carry about as much information as 150 independent people. The formula also shows why cluster size matters: with 32 people per cluster, even an ICC of 0.05 gives a design effect of 1 + 31 × 0.05 = 2.55, which reduces the effective sample size by more than half. Design effects are also used when planning studies, to inflate the sample size needed for a cluster-sampled survey or a cluster-randomized trial. Students who took HSCI 230 or HSCI 341 met this formula in HSCI 230 Lesson 5, Section 6, and HSCI 341 Lesson 2, Section 5, where it inflates the sample size of a planned study; this section uses the same formula to describe the information lost in an analysis.

The effective sample size of about 150 summarises how much information the clustering removes from the clinic data for questions about clinic characteristics. The next part shows how this loss appears in a regression model that ignores the clinics.

What Goes Wrong When Clustering Is Ignored

The question of whether patients in urban clinics have higher blood pressure than patients in rural clinics shows the problem directly. An ordinary linear regression with lm() treats the 966 patients as independent. A linear mixed model with lmer() gives each clinic its own baseline, as Section 2 explains.

Figure 5.3 compares the two estimates of the urban minus rural difference with their 95% confidence intervals, and Table 5.2 gives the estimates, standard errors and p-values from the two models.

Two estimates of the urban minus rural difference in systolic blood pressure with 95% confidence intervals. Ignoring clinics: 1.49 (0.04 to 2.95). Allowing for clinics: 1.67 (minus 2.17 to 5.51).
Figure 5.3. The two lines show the estimated difference in blood pressure between urban and rural clinics with its 95% confidence interval. The estimates are similar, but the interval from the mixed model is about two and a half times as wide and includes 0.

Table 5.2. The urban minus rural difference in systolic blood pressure from an ordinary regression and a linear mixed model.

ModelUrban minus rural (mmHg)Standard errorp-value
Ordinary regression, lm()1.490.740.044
Mixed model, lmer()1.671.960.40

The ordinary regression suggests a statistically significant difference, and the mixed model shows that the data cannot distinguish urban from rural clinics. The standard error rises from 0.74 to 1.96 because clinic_urban is a cluster-level predictor: it takes one value per clinic, so the 966 patients supply only 30 independent values of it. A person-level predictor such as smoking varies among patients within each clinic, and its standard error changes little when clustering is taken into account (in the Section 2 model, which adds age, gender and smoking, it falls from 0.79 to 0.70, and Section 3 shows a similar pattern for referrals). When a predictor varies within clusters, allowing for clustering can make its estimate more precise, because the model can compare people within the same clinic.

Box 5.1 adds a caution about how this problem is detected in practice.

⚠ Box 5.1: The problem cannot be seen in the ordinary regression output

Nothing in the output of lm() or glm() warns that the observations are clustered. Dependence has to be recognised from the study design and from a variable that identifies the clusters, which is why the first step in any analysis is to ask how the data were collected.

A clustered dataset therefore needs a method that allows for the clustering, and the next part introduces the two families of methods that the lesson uses.

Methods for Dependent Data

Two families of methods are used throughout this lesson. Mixed models add a random effect for each cluster and estimate how much the clusters vary; GEE fits a regression with no cluster effects, so the results describe the population as a whole, and calculates standard errors that allow for the correlation within clusters. Both need a variable that identifies the clusters. Simpler approaches also exist, such as analysing cluster averages or adding a separate indicator variable for every cluster, but they lose information or cannot estimate the effects of cluster-level predictors, and this lesson does not use them.

Table 5.3 lists the three models that the lesson fits, with the outcomes they handle, the R function for each and the sections in which they are used.

Table 5.3. Methods for dependent data used in this lesson.

MethodOutcomeR functionUsed in
Linear mixed modelContinuouslme4::lmer() (with lmerTest for p-values)Sections 2 and 4
Generalized linear mixed model (GLMM)Binary or countlme4::glmer()Section 3
Generalized estimating equations (GEE)Continuous, binary or countgeepack::geeglm()Sections 3 and 4

This section has established that dependence arises from the design of a study, that the ICC and the design effect measure its strength, and that ignoring it understates the uncertainty for cluster-level predictors. The narrated R walkthrough below demonstrates all of the methods in Table 5.3 on the lesson’s two datasets. Activity 5.1 then reproduces the main results of this section, from the ICC to the comparison in Table 5.2, before the knowledge check for the section.

Learn to do this in R

Narrated R walkthrough: Mixed Models and GEE in R

This walkthrough uses the same two datasets as this lesson. It measures clustering with the intracluster correlation coefficient, fits and checks a linear mixed model, models a binary outcome with a generalized linear mixed model and with GEE, and analyses repeated measurements over time.

Open the Mixed Models and GEE in R walkthrough
R Activity 5.1: detecting clustering in the clinic data

This activity measures the clustering in the clinic data and compares an ordinary regression with a mixed model. Download phaa_clinics.csv and save it in the folder that holds your R script, then choose Session → Set Working Directory → To Source File Location in RStudio. The lesson uses the lme4, lmerTest and geepack packages, which can be installed once with the commented line at the top of the code.

# install.packages(c("lme4", "lmerTest", "geepack"))   # run once, if not yet installed
library(lmerTest)    # loads lme4 (for lmer()) and adds p-values
clinics <- read.csv("phaa_clinics.csv")   # file must be in the working directory
clinics$clinic_id <- factor(clinics$clinic_id)
clinics$smoker <- factor(clinics$smoker, levels = c("No", "Yes"))
clinics$clinic_urban <- factor(clinics$clinic_urban, levels = c("rural", "urban"))
nlevels(clinics$clinic_id)                        # how many clinics?
summary(as.vector(table(clinics$clinic_id)))      # patients per clinic
Console output
[1] 30 Min. 1st Qu. Median Mean 3rd Qu. Max. 18.00 26.25 32.50 32.20 39.75 45.00

The factor() line tells R to treat clinic_id as a label for each clinic. The file holds 30 clinics with between 18 and 45 patients each (mean 32.2).

m0 <- lmer(sbp ~ 1 + (1 | clinic_id), data = clinics)   # a model with clinics only
vc <- as.data.frame(VarCorr(m0))
vc[, c("grp", "vcov")]                    # variance between clinics and within clinics
icc <- vc$vcov[1] / sum(vc$vcov)          # intracluster correlation coefficient
icc
Console output
grp vcov 1 clinic_id 22.10395 2 Residual 104.34414 [1] 0.1748065

The term (1 | clinic_id) asks for a separate baseline for each clinic. The variance between clinics is 22.10 and the variance within clinics (Residual) is 104.34, so the ICC is 22.10 ÷ 126.44 = 0.175.

m_bar <- mean(table(clinics$clinic_id))   # average clinic size
deff  <- 1 + (m_bar - 1) * icc             # design effect
deff
nrow(clinics) / deff                       # effective sample size
Console output
[1] 6.453962 [1] 149.6755

The design effect is 6.45, and the effective sample size is about 150 independent people.

naive <- lm(sbp ~ clinic_urban, data = clinics)                     # ignores clinics
mixed <- lmer(sbp ~ clinic_urban + (1 | clinic_id), data = clinics)  # allows for clinics
summary(naive)$coefficients
summary(mixed)$coefficients
Console output
Estimate Std. Error t value Pr(>|t|) (Intercept) 110.671233 0.5853981 189.052936 0.00000000 clinic_urbanurban 1.493493 0.7421687 2.012336 0.04446231 Estimate Std. Error df t value Pr(>|t|) (Intercept) 110.65727 1.592412 26.91822 69.4903668 6.566440e-32 clinic_urbanurban 1.66819 1.959339 27.40708 0.8514043 4.019229e-01

Reading the output. The first table is from lm(): urban clinics average 1.49 mmHg higher, with a standard error of 0.74 and a p-value of 0.044. The second table is from lmer(): the difference is 1.67 mmHg, but the standard error is 1.96 and the p-value is 0.40. The df column (about 27) reflects the 30 clinics, which are the real units of information for a clinic-level predictor.

R Reflect on what you just ran

Use the questions below to interpret the output you produced. Look at your console output before answering.

1. Report the ICC from m0 and explain in one sentence what it means for two patients who attend the same clinic.

Model answerThe ICC is 0.175. About 17.5% of the variation in systolic blood pressure lies between clinics, so the blood pressures of two patients from the same clinic are correlated at about 0.175, and they tend to be more similar than the blood pressures of two patients from different clinics.

2. Report the design effect and the effective sample size, and explain what the effective sample size tells you about the 966 patients.

Model answerThe design effect is 6.45 and the effective sample size is 966 ÷ 6.45 = about 150. For questions about clinic characteristics, the 966 clustered patients carry about as much information as 150 independent people, because patients in the same clinic partly repeat the same information.

3. Compare the estimate, standard error and p-value for clinic_urbanurban from naive and mixed. Which result would you report, and why?

Model answerThe ordinary regression gives 1.49 mmHg (SE 0.74, p = 0.044) and the mixed model gives 1.67 mmHg (SE 1.96, p = 0.40). The mixed model should be reported, because clinic_urban is a clinic-level predictor with only 30 independent values, and the ordinary regression treats it as if it had 966. The ordinary regression understates the uncertainty, and the mixed model shows that the data cannot tell whether urban clinics differ from rural clinics.
Saved.
Knowledge check: this section

1. Which study produces repeated-measures data?

When the same person is measured several times, the measurements from that person are a cluster and are not independent of one another.

2. What usually happens to the standard error of a cluster-level predictor when clustering is ignored?

Ignoring clustering treats every observation as independent information, so the standard error of a cluster-level predictor is too small, and its confidence interval and p-value are too small as well.

3. A dataset has a between-cluster variance of 5 and a within-cluster variance of 45. What is the ICC?

ICC = 5 ÷ (5 + 45) = 5 ÷ 50 = 0.10.

4. Clusters contain 21 people on average and the ICC is 0.05. What is the design effect?

DEFF = 1 + (21 − 1) × 0.05 = 1 + 1.0 = 2.0, so the effective sample size is half the actual sample size.

5. In a study of patients in clinics, which variable is a cluster-level predictor?

Whether the clinic has evening hours takes the same value for every patient in a clinic, so it describes the cluster. Age and smoking vary among patients within a clinic.

✎ Reflection

This section explained that data are dependent when observations share a cluster, such as patients in a clinic or repeated visits by the same person. The intracluster correlation coefficient (ICC) is the share of the outcome's variation that lies between clusters, and the design effect, 1 + (average cluster size − 1) × ICC, shows how much information is lost; the effective sample size is the sample size divided by the design effect. Ignoring dependence makes standard errors too small, especially for cluster-level predictors that take one value per cluster. A health unit surveys 40 students in each of 25 schools (1,000 students in all) about their physical activity and wants to know whether schools with a daily physical education policy have more active students. The ICC for physical activity is 0.08. Calculate the design effect and the effective sample size, explain whether the school policy is a cluster-level or a person-level predictor, and describe what would go wrong if the analysis treated the 1,000 students as independent.

Model answerThe design effect is 1 + (40 − 1) × 0.08 = 1 + 3.12 = 4.12, so the effective sample size is 1,000 ÷ 4.12 = about 243 students. The daily physical education policy is a cluster-level predictor, because it takes the same value for every student in a school, so the data hold only 25 independent values of it, one per school. If the analysis treated the 1,000 students as independent, it would act as if it had about four times as much information as it does. The standard error for the policy would be too small, the confidence interval too narrow and the p-value too small, so the analysis could report that the policy is associated with more activity when the 25 schools cannot support that conclusion. A mixed model with a random intercept for school, or GEE with school as the cluster, would give an appropriate standard error.
✓ Reflection saved!
● Complete the quiz and reflection to continue.
Section 2 of 4

Linear Mixed Models

⏱ Estimated time: 40 minutes
Lesson 5 · Section 2

Linear Mixed Models

A model for a continuous outcome that gives each cluster its own baseline.

The idea

Each clinic gets its own baseline

Random-intercept linear mixed model
\[ \text{SBP}_{ij} = \color{#6D28D9}{\beta_0} + \color{#C2410C}{u_j} + \color{#1D4ED8}{\beta_1}\,\text{age}_{ij} + \dots + \color{#0B7B6B}{\varepsilon_{ij}}, \qquad \color{#C2410C}{u_j} \sim N(0, \sigma^2_u) \]
β0 overall intercept uj shift for clinic j (random effect) β1, … fixed effects, read as in Lesson 4 εij difference for patient i within the clinic

The model estimates one variance, σ²u, that describes how much the clinic baselines differ.

Two kinds of effects

Fixed effects and random effects

Fixed effects

These are the coefficients for the predictors, read as differences in average blood pressure, as in Lesson 4.

Random effects

These are the clinic shifts, summarised by their variance or standard deviation.

lmer(sbp ~ age + female + smoker + clinic_urban + (1 | clinic_id), data = clinics)

The term (1 | clinic_id) asks for a random intercept for each clinic.

Reading the output

Fixed effects from lmer()

PredictorEstimate (mmHg)95% CIp-value
Age (per year)0.410.37 to 0.45< 0.001
Female−2.80−3.88 to −1.73< 0.001
Smoker5.003.63 to 6.37< 0.001
Urban clinic1.85−1.85 to 5.540.34

The random effects show a standard deviation of 4.6 mmHg between clinic baselines and 8.3 mmHg among patients within a clinic.

The clinic effects

Clinics differ by several mmHg, after adjustment

Estimated clinic effects for 30 clinics, sorted from lowest to highest, with 95% intervals. Most lie within 5 mmHg of zero; three clinics lie about 8 to 11 mmHg above average.
The adjusted ICC is 21.6 ÷ (21.6 + 68.9) = 0.24, and each point is one clinic's estimated shift from the overall average.
Interpretation

What a fixed effect means

  • The smoking coefficient compares a smoker and a non-smoker of the same age and gender in the same clinic.
  • For a continuous outcome, this is also the average difference in the whole population.
  • The main gain from the mixed model is a standard error that reflects the clustering.
Model checking

Three plots check most of the assumptions

Three diagnostic plots for the linear mixed model: residuals versus fitted values with a flat smoother, a Q-Q plot of residuals close to the line, and a Q-Q plot of the 30 clinic effects with three high points above the line.
Panels A and B look acceptable; panel C shows three clinics with higher baselines than a normal distribution predicts.
Further checks

Enough clusters, and no singular fit

AssumptionHow to check itThis example
Linearity and equal varianceResiduals versus fitted valuesAcceptable
Normal residualsQ-Q plot of residualsAcceptable
Normal clinic effectsQ-Q plot of the random effectsThree high clinics
Enough clusters (about 20 to 30)nlevels(clinic_id)30 clinics
Clinic variance estimatedisSingular(model)FALSE
Going further

Random slopes

Random intercept

Clinic baselines differ, and each predictor has the same effect in every clinic.

(1 | clinic_id)

Random slope

The effect of a predictor, such as age, is also allowed to differ between clinics.

(age | clinic_id)

Random slopes need more clusters and often trigger warnings with small samples.

Carry forward

What to take into the next section

  • A random intercept gives each cluster its own baseline, summarised by one variance.
  • Fixed effects are read as in an ordinary regression, and random effects describe the clusters.
  • The checks are residual plots, a Q-Q plot of cluster effects, the number of clusters and the singular-fit warning.

Introduction and Overview

This section presents the linear mixed model, the standard model for a continuous outcome when observations are clustered. It continues with systolic blood pressure in the 30 clinics of phaa_clinics.csv, now with four predictors: age, gender, smoking and whether the clinic is urban. The section explains the random intercept, shows how to read the fixed and random parts of the R output, lists the checks that a mixed model needs, and describes random slopes briefly.

Learning Objectives

  • Describe a random intercept and explain how a linear mixed model differs from an ordinary linear regression.
  • Fit a random-intercept model in R with lmer() and interpret its fixed effects.
  • Interpret the random-effect variances and calculate the adjusted ICC.
  • Check a linear mixed model with residual plots, a Q-Q plot of the cluster effects, the number of clusters and the singular-fit check.
  • Explain what a random slope adds and when it is used.

The Random-Intercept Model

An ordinary linear regression has a single intercept. A random-intercept model gives each cluster its own intercept, made up of the overall intercept plus a cluster-specific shift. The shifts are treated as a random sample from a normal distribution with mean 0, and the model estimates the variance of that distribution, which describes how much the clusters differ (Laird & Ware, 1982). Equation 5.3 writes the model for systolic blood pressure in the clinic data.

Random-intercept model for systolic blood pressure
\[ \text{SBP}_{ij} = \color{#6D28D9}{\beta_0} + \color{#C2410C}{u_j} + \color{#1D4ED8}{\beta_1}\,\text{age}_{ij} + \color{#1D4ED8}{\beta_2}\,\text{female}_{ij} + \color{#1D4ED8}{\beta_3}\,\text{smoker}_{ij} + \color{#1D4ED8}{\beta_4}\,\text{urban}_{j} + \color{#0B7B6B}{\varepsilon_{ij}} \]Eq 5.3
The blood pressure of patient i in clinic j equals the overall intercept, plus a shift for clinic j, plus the fixed effects of the predictors, plus a difference for the patient within the clinic. The clinic shifts uj have variance σ²u, and the patient differences εij have variance σ²e. The subscript j on urban shows that it is a clinic-level predictor.

The name "mixed" refers to the two kinds of effects in the model. The fixed effects (β0 to β4) are the coefficients for the predictors, read as differences in average blood pressure, as in Lesson 4. The random effects (uj) are the clinic shifts, summarised by their variance. Summarising the 30 clinics with one variance lets the model estimate the effect of a clinic-level predictor such as clinic_urban.

The three flip cards below restate the random intercept, the fixed effects and the variance components, with the lmer() notation for the random intercept.

Random InterceptClick to explore
Fixed EffectsClick to explore
Variance ComponentsClick to explore

With the parts of the model defined, the next part reads them from the R output for the clinic data.

Reading the Output

Activity 5.2 below fits the model with lmer() after loading the lmerTest package, which adds p-values to the fixed effects. The output has two main blocks. Table 5.4 summarises the fixed-effects block, with a 95% confidence interval and an interpretation for each predictor.

Table 5.4. Fixed effects from the random-intercept model for systolic blood pressure.

Fixed effectEstimate (mmHg)95% CIp-valueInterpretation
Age (per year)0.410.37 to 0.45< 0.001Each extra year of age is associated with 0.41 mmHg higher blood pressure.
Female−2.80−3.88 to −1.73< 0.001Women average 2.8 mmHg lower than men.
Smoker5.003.63 to 6.37< 0.001Smokers average 5.0 mmHg higher than non-smokers.
Urban clinic1.85−1.85 to 5.540.34There is no clear evidence of a difference between urban and rural clinics.

Each estimate holds the other predictors constant. The df column of the output shows about 937 degrees of freedom for the patient-level predictors and about 28 for urban clinics, because the information about a clinic-level predictor comes from the 30 clinics.

The random effects block reports a clinic variance of 21.57 (standard deviation 4.64 mmHg) and a residual variance of 68.94 (standard deviation 8.30 mmHg). The adjusted ICC is 21.57 ÷ (21.57 + 68.94) = 0.24. It is larger than the ICC of 0.175 from the model with no predictors, because age, gender and smoking explain much of the variation among patients within clinics but little of the variation between clinics. Figure 5.4 shows the estimated shifts of the 30 clinic baselines, which the clinic variance summarises.

Estimated clinic effects for 30 clinics with 95% intervals, sorted from lowest (about minus 8 mmHg) to highest (about plus 11 mmHg). Most lie within 5 mmHg of zero.
Figure 5.4. Each point is the estimated shift of one clinic's baseline from the overall average, after adjusting for the predictors, with a 95% interval. Most clinics lie within 5 mmHg of the average, and three lie about 8 to 11 mmHg above it.

Worked Example 5.2 uses the intercept and the fixed effects of the model to predict blood pressure for two patients, and shows how a clinic shift such as those in Figure 5.4 changes the predictions.

Worked Example 5.2: Predicted Blood Pressure for Two Patients

For a 50-year-old male non-smoker in a rural clinic, the model predicts 90.01 + 0.408 × 50 = 110.4 mmHg in an average clinic. A 50-year-old male smoker in the same clinic is predicted to be 5.0 mmHg higher, at 115.4 mmHg. In a clinic whose estimated shift is +8 mmHg, both predictions rise by 8 mmHg, but the difference between them stays 5.0 mmHg, because the clinic shift is shared by everyone in the clinic.

These interpretations depend on the model being appropriate for the data, so the next part sets out the assumptions of a linear mixed model and how to check them.

Assumptions and How to Check Them

A linear mixed model makes the assumptions of linear regression, applied to the residuals within clusters, and adds assumptions about the clusters. Table 5.5 lists each assumption with what it means and how to check it, and Figure 5.5 shows the three diagnostic plots for the blood pressure model.

Table 5.5. Assumptions of a linear mixed model and how to check them.

AssumptionWhat it meansHow to check it
Linearity and equal varianceEach continuous predictor has a straight-line relationship with the outcome, and the residuals have the same spread at every fitted value.Plot of residuals against fitted values.
Normal residualsThe differences among patients within a clinic are roughly normal.Q-Q plot of the residuals.
Normal cluster effectsThe clinic shifts come from a roughly normal distribution.Q-Q plot of ranef(model). With few clusters this check is rough.
Enough clustersThe cluster variance is estimated from the clusters, so a reasonable number is needed.nlevels(clinic_id); about 20 to 30 clusters is a common minimum.
Independent clustersClinics are unrelated to one another, and patients are independent within a clinic once the clinic effect is allowed for.Judged from the study design.
Cluster variance estimatedThe model is not singular, so the cluster variance is not estimated as zero.isSingular(model) returns FALSE; a "boundary (singular) fit" message signals a problem.
Three diagnostic plots for the linear mixed model. A: residuals versus fitted values form a flat, even cloud. B: residuals follow the Q-Q line closely. C: the 30 clinic effects follow the Q-Q line except for three high points.
Figure 5.5. Panel A supports linearity and equal variance, and panel B supports normal residuals. Panel C shows that three clinics have higher baselines than a normal distribution would predict, which is worth reporting and investigating.

In this example the residual checks look acceptable, there are 30 clinics, and the model is not singular. With fewer than about 40 clusters, a small-sample correction to the tests is advised (Hemming & Taljaard, 2023), and the lmerTest output in this section already uses one (Satterthwaite's method). The three high clinics in panel C are a mild departure from normal cluster effects. Fixed-effect estimates are usually not very sensitive to this assumption, but the clinics themselves might be examined, for example to see whether they serve an older population or record blood pressure differently.

The checks support the random-intercept model for the clinic data. The final part of this section considers a more flexible model, in which the effect of a predictor may differ between clinics.

Random Slopes

The random-intercept model lets clinic baselines differ but assumes that each predictor has the same effect in every clinic. A random slope lets the effect of a predictor vary between clusters as well. In lmer(), a random slope for age is written (age | clinic_id), which gives each clinic its own intercept and its own age slope. Random slopes suit research questions about how an effect differs between settings, such as whether a program works better in some schools than others. They need more clusters and more observations per cluster, and with 30 clinics they often produce a singular fit or a convergence warning. This lesson fits random intercepts only.

This section has fitted a linear mixed model with a random intercept for each clinic, shown how to read its fixed effects, variance components and adjusted ICC, and checked it against the assumptions in Table 5.5. Activity 5.2 fits the same model in R and runs the checks, and the knowledge check that follows tests the main ideas of the section. Section 3 applies the same approach to a binary outcome.

R Activity 5.2: a random-intercept model for blood pressure

This activity fits the linear mixed model and runs its checks. It loads the data again, so it can be run on its own.

library(lmerTest)   # loads lme4 and adds p-values to lmer() output
clinics <- read.csv("phaa_clinics.csv")   # file must be in the working directory
clinics$clinic_id <- factor(clinics$clinic_id)
clinics$smoker <- factor(clinics$smoker, levels = c("No", "Yes"))
clinics$clinic_urban <- factor(clinics$clinic_urban, levels = c("rural", "urban"))
lmm <- lmer(sbp ~ age + female + smoker + clinic_urban + (1 | clinic_id), data = clinics)
summary(lmm)
Console output
Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmerModLmerTest] Formula: sbp ~ age + female + smoker + clinic_urban + (1 | clinic_id) Data: clinics REML criterion at convergence: 6896.6 Scaled residuals: Min 1Q Median 3Q Max -2.95086 -0.65976 -0.02077 0.66101 2.98357 Random effects: Groups Name Variance Std.Dev. clinic_id (Intercept) 21.57 4.644 Residual 68.94 8.303 Number of obs: 966, groups: clinic_id, 30 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 90.0114 1.8891 62.2481 47.647 < 2e-16 *** age 0.4076 0.0201 937.0898 20.282 < 2e-16 *** female -2.8050 0.5465 937.4243 -5.133 3.47e-07 *** smokerYes 4.9993 0.6975 935.9141 7.168 1.55e-12 *** clinic_urbanurban 1.8470 1.8851 27.6963 0.980 0.336 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) age female smkrYs age -0.555 female -0.162 -0.006 smokerYes -0.077 0.009 0.008 clnc_rbnrbn -0.663 0.002 0.000 0.006

Reading the output. The first lines (the REML criterion and the scaled residuals) can be ignored for now. The Random effects block gives the clinic variance (21.57, standard deviation 4.644) and the residual variance among patients (68.94, standard deviation 8.303), and the line below it confirms 966 patients in 30 clinics. The Fixed effects block is read like a regression table: age 0.41 mmHg per year, women 2.80 mmHg lower, smokers 5.00 mmHg higher, and urban clinics 1.85 mmHg higher with p = 0.336. The Correlation of Fixed Effects block can be ignored for now.

confint(lmm, parm = "beta_", method = "Wald")   # 95% CIs for the fixed effects
vc <- as.data.frame(VarCorr(lmm))
vc$vcov[1] / sum(vc$vcov)                        # ICC after adjusting for the predictors
Console output
2.5 % 97.5 % (Intercept) 86.3087770 93.7140882 age 0.3682456 0.4470293 female -3.8760676 -1.7338623 smokerYes 3.6322803 6.3663724 clinic_urbanurban -1.8476208 5.5416933 [1] 0.2382825

The 95% confidence intervals agree with the p-values: only the interval for urban clinics (−1.85 to 5.54) includes 0. The adjusted ICC is 0.24. The intervals from method = "Wald" use the normal distribution, while the lmerTest p-values use a t distribution with the degrees of freedom shown, so for urban clinics the two can differ slightly.

par(mfrow = c(1, 3))
plot(fitted(lmm), resid(lmm), main = "Residuals vs fitted"); abline(h = 0, lty = 2)
qqnorm(resid(lmm), main = "Q-Q: residuals"); qqline(resid(lmm))
qqnorm(ranef(lmm)$clinic_id[, 1], main = "Q-Q: clinic effects"); qqline(ranef(lmm)$clinic_id[, 1])
par(mfrow = c(1, 1))
isSingular(lmm)        # FALSE means the clinic variance was estimated without problems
Console output
[1] FALSE

The three plots match Figure 5.5 earlier in this section, and isSingular() returns FALSE, so the clinic variance was estimated without problems.

R Reflect on what you just ran

Use the questions below to interpret the output you produced. Look at your console output and plots before answering.

1. Report the fixed effect for smokerYes with its 95% confidence interval, and explain it in one sentence.

Model answerThe fixed effect is 5.00 mmHg (95% CI 3.63 to 6.37). Smokers have an average systolic blood pressure 5.0 mmHg higher than non-smokers of the same age and gender attending the same clinic, and the interval excludes 0.

2. Use the random-effects block to calculate the adjusted ICC, and explain why it is larger than the ICC of 0.175 from the model with no predictors.

Model answerThe adjusted ICC is 21.57 ÷ (21.57 + 68.94) = 0.24. It is larger than 0.175 because age, gender and smoking explain a good part of the variation among patients within clinics, which shrinks the residual variance from about 104 to 69, while they explain little of the variation between clinics, which stays at about 22. The clinic share of the remaining variation therefore rises.

3. Describe what the three diagnostic plots and isSingular(lmm) show, and state whether the model is acceptable.

Model answerThe residuals-versus-fitted plot is a flat cloud with an even spread, which supports linearity and equal variance. The Q-Q plot of the residuals follows the line, which supports normal residuals. The Q-Q plot of the clinic effects follows the line except for three clinics with higher baselines than expected. isSingular() returns FALSE, so the clinic variance was estimated. The model is acceptable, and the three high clinics are worth noting as a limitation and examining further.
Saved.
Knowledge check: this section

1. What does the term (1 | clinic_id) add to an lmer() model?

The term gives each clinic its own baseline (intercept), with the clinic shifts treated as random effects from a normal distribution.

2. In a linear mixed model of blood pressure, the fixed effect for smoking is 5.0 mmHg. What does it mean?

Fixed effects in a linear mixed model are read as in linear regression: a difference in the average outcome, holding the other predictors constant.

3. A linear mixed model reports a clinic variance of 12 and a residual variance of 36. What is the ICC?

ICC = 12 ÷ (12 + 36) = 12 ÷ 48 = 0.25.

4. Which plot checks the assumption that the cluster effects are normally distributed?

A Q-Q plot of the estimated cluster effects (from ranef()) checks whether they follow a normal distribution.

5. R reports "boundary (singular) fit" for a mixed model. What does this usually mean?

A singular fit usually means that a variance, often the cluster variance, has been estimated at or near zero, because the data show little clustering or the model is too complex for the data.

✎ Reflection

This section introduced the linear mixed model, which gives each cluster its own baseline (a random intercept, written (1 | clinic_id) in lmer()). Its fixed effects are read as differences in the average outcome, as in linear regression, and its random effects are summarised by a cluster variance and a residual variance, whose ratio gives the ICC. The model is checked with a residuals-versus-fitted plot, Q-Q plots of the residuals and of the cluster effects, the number of clusters (about 20 to 30 is a common minimum) and the singular-fit check. A researcher studies body mass index (BMI) among 1,200 employees in 15 workplaces, with age, gender and whether the workplace has an on-site cafeteria as predictors. Explain which terms you would include in an lmer() model and which of them are fixed or random, which predictor is at the workplace level, how you would interpret the cafeteria coefficient, and which assumption checks might be a concern with 15 workplaces.

Model answerThe model would be lmer(bmi ~ age + female + cafeteria + (1 | workplace_id), data = ...). Age, gender and cafeteria are fixed effects, and (1 | workplace_id) is a random intercept that gives each workplace its own baseline BMI. The cafeteria variable is a workplace-level predictor, because it takes one value for every employee in a workplace, so its information comes from the 15 workplaces. Its coefficient would be the difference in average BMI between employees of workplaces with and without a cafeteria, holding age and gender constant, and its confidence interval would be wide because there are only 15 workplaces. With 15 workplaces, the number of clusters is below the usual minimum of about 20 to 30, so the workplace variance would be estimated imprecisely and the model could be singular, and the Q-Q plot of 15 workplace effects would be hard to judge. I would still check the residual plots, report the number of workplaces as a limitation and consider whether more workplaces could be recruited.
✓ Reflection saved!
● Complete the quiz and reflection to continue.
Section 3 of 4

Binary Outcomes: Mixed Models and GEE

⏱ Estimated time: 35 minutes
Lesson 5 · Section 3

Binary Outcomes: Mixed Models and GEE

Two ways to allow for clustering when the outcome is yes or no.

The outcome

Specialist referral in 30 clinics

246
of 966 patients referred
25%
overall referral rate
30
clinics, with different referral rates

Referral is binary, so each method is a version of logistic regression that allows for the clinics.

Approach 1

The GLMM: logistic regression with a random intercept

glmer(referred ~ age + smoker + clinic_urban + (1 | clinic_id), family = binomial, data = clinics)

What it estimates

The odds ratios compare two patients in the same clinic, or in two clinics with the same baseline, so they are clinic-specific (conditional) odds ratios.

Clinic variation

The clinic shifts on the log-odds scale have a variance of 0.30 (standard deviation 0.55).

Approach 2

GEE: a population-averaged model with corrected standard errors

geeglm(referred ~ age + smoker + clinic_urban, id = clinic_id, family = binomial, corstr = "exchangeable", data = clinics)

What it estimates

The odds ratios describe the population as a whole, so they are population-averaged (marginal) odds ratios.

Working correlation

An exchangeable structure assumes that any two patients in a clinic are equally correlated, and here the correlation is estimated at 0.04.

Comparing the models

Smoking barely changes; urban clinic does

Odds ratios with 95% confidence intervals from three models. Smoking: 1.72, 1.67 and 1.64 with similar intervals. Urban clinic: 2.16 from glm, 2.31 from glmer and 2.04 from geeglm, with wider intervals for the last two.
The standard errors for urban clinic on the log scale are 0.17 for glm, 0.28 for the GLMM and 0.29 for GEE.
Two questions

Clinic-specific or population-averaged?

GLMM (clinic-specific)

This compares patients in the same clinic, or in clinics with the same baseline. Its odds ratios are usually a little further from 1. It suits questions about individual clinics.

GEE (population-averaged)

This compares referral across the whole population, averaged over all clinics. It suits questions about policy for the population.

For a continuous outcome, as in Section 2, the two kinds of estimate coincide.

Assumptions and checks

What to check in a GLMM and in GEE

CheckGLMM: glmer()GEE: geeglm()
Binary outcome and enough eventsYes (about 10 per predictor)Yes
Enough clustersAbout 20 to 30About 30 to 40 or more
Cluster effectsAssumed normal (log-odds scale)No assumption about their shape
WarningsRead convergence warningsCheck the cluster id and sort the data
Choosing a method

GLMM or GEE?

QuestionPoints to
Is the question about individual clusters?GLMM
Is the question about the whole population?GEE
Are there fewer than about 30 clusters?GLMM
Is the variation between clusters of interest?GLMM
Is the outcome continuous?Linear mixed model (Section 2)

Both methods agree here: urban clinics refer more patients.

Carry forward

What to take into the next section

  • A GLMM adds a random intercept to a logistic regression, and GEE fits a population-averaged logistic regression with standard errors that allow for clustering.
  • The GLMM gives clinic-specific odds ratios and GEE gives population-averaged odds ratios.
  • Both methods widen the standard errors of cluster-level predictors, and GEE needs more clusters.

Introduction and Overview

This section extends logistic regression to clustered data. The outcome is referred, which records whether each of the 966 patients in phaa_clinics.csv was referred to a specialist (246 were, or 25.5%). The predictors are age, smoking and urban clinic. The section fits a generalized linear mixed model (GLMM) and a generalized estimating equations (GEE) model, compares them with an ordinary logistic regression, and explains when each approach is preferred.

Learning Objectives

  • Fit a logistic GLMM with glmer() and a logistic GEE model with geeglm() in R.
  • Interpret odds ratios from both models and compare their standard errors with an ordinary logistic regression.
  • Distinguish clinic-specific (conditional) from population-averaged (marginal) odds ratios in plain language.
  • Choose between a GLMM and GEE using the research question and the number of clusters.
  • List the checks for each approach.

Two Ways to Allow for Clustering

Logistic regression can allow for clustering in two ways, which correspond to the two families of methods introduced in Section 1 and listed in Table 5.3. A GLMM gives each clinic its own baseline log-odds, and GEE fits one model for the whole population and calculates standard errors that allow for the clustering. The two subsections below describe each approach in turn.

The Generalized Linear Mixed Model

A generalized linear mixed model (GLMM) adds random effects to a generalized linear model. For a binary outcome, it is a logistic regression in which each clinic has its own baseline log-odds, with the clinic shifts drawn from a normal distribution. It is fitted with glmer() from lme4, which takes the same family = binomial argument as glm(). Equation 5.4 writes the GLMM for referral in the clinic data.

Logistic GLMM for referral
\[ \ln\frac{\color{#0B7B6B}{p_{ij}}}{1 - \color{#0B7B6B}{p_{ij}}} = \color{#6D28D9}{\beta_0} + \color{#C2410C}{u_j} + \color{#1D4ED8}{\beta_1}\,\text{age}_{ij} + \color{#1D4ED8}{\beta_2}\,\text{smoker}_{ij} + \color{#1D4ED8}{\beta_3}\,\text{urban}_{j} \]Eq 5.4
The log-odds of referral for patient i in clinic j equals the overall intercept plus a clinic shift plus the fixed effects. Exponentiating a fixed effect gives an odds ratio for patients in the same clinic, or in clinics with the same baseline.

Generalized Estimating Equations

Generalized estimating equations (GEE) fit a logistic regression without clinic effects. A working correlation describes how alike patients in the same clinic are assumed to be, and GEE uses it to weight the data. The standard errors are then calculated from the variation between clinics, which allows for the clustering (Zeger & Liang, 1986). The analyst chooses the working correlation structure. For people in clinics, schools or households, the usual choice is "exchangeable", which assumes that any two members of a cluster are equally correlated. The corrected (sandwich) standard errors remain valid even if the working correlation is not exactly right, provided there are enough clusters. GEE is fitted with geeglm() from geepack, which needs the cluster identifier in the id argument and the data sorted so that each cluster's rows are together.

The three flip cards below define terms that recur in the rest of this section: the clinic-specific odds ratio that a GLMM estimates, the population-averaged odds ratio that GEE estimates, and the working correlation that GEE uses.

Clinic-Specific Odds RatioClick to explore
Population-Averaged Odds RatioClick to explore
Working CorrelationClick to explore

The two approaches therefore estimate different quantities: the GLMM gives clinic-specific odds ratios and GEE gives population-averaged odds ratios. The next part fits both models to the referral data and compares their results with an ordinary logistic regression.

Comparing the Results

Activity 5.3 below fits an ordinary logistic regression (glm()), a GLMM (glmer()) and GEE (geeglm()) with the same three predictors. Table 5.6 shows the odds ratios for smoking (a patient-level predictor) and urban clinic (a clinic-level predictor), with their standard errors on the log-odds scale.

Table 5.6. Odds ratios for smoking and urban clinic, with standard errors on the log-odds scale, from three models of specialist referral.

ModelOR for smokingSE (smoking)OR for urban clinicSE (urban clinic)
Ordinary logistic regression, glm()1.720.1882.160.168
GLMM, glmer()1.670.1942.310.279
GEE, geeglm()1.640.1862.040.289

Figure 5.6 plots the same odds ratios with their 95% confidence intervals.

Odds ratios with 95% confidence intervals from glm, glmer and geeglm. Smoking: 1.72, 1.67 and 1.64 with similar intervals. Urban clinic: 2.16, 2.31 and 2.04, with the glmer and geeglm intervals much wider than the glm interval.
Figure 5.6. The left panel shows the odds ratios for smoking and the right panel the odds ratios for urban clinics, each with a 95% confidence interval. Allowing for clinics barely changes the smoking result but widens the interval for urban clinics.

For smoking, the three models give similar odds ratios and standard errors, because smokers and non-smokers are found in every clinic and the comparison can be made within clinics. For urban clinic, the standard error rises from 0.168 in the ordinary logistic regression to 0.279 in the GLMM and 0.289 in GEE, an increase of about 70%. The odds ratio remains clearly above 1 in all three models: patients in urban clinics have about twice the odds of referral of patients in rural clinics of the same age and smoking status (GLMM OR 2.31, p = 0.003; GEE OR 2.04, p = 0.014).

The GEE odds ratio for urban clinics (2.04) is a little lower than the glm() odds ratio (2.16) because the exchangeable working correlation gives each patient in a large clinic slightly less weight. Refitting GEE with corstr = "independence" reproduces the glm() odds ratio of 2.16, with the corrected standard error of 0.289.

Worked Example 5.3 explains why the GLMM odds ratio for urban clinics in Table 5.6 is further from 1 than the GEE odds ratio.

Worked Example 5.3: Why the GLMM Odds Ratio Is Further From 1

Suppose that a patient in an urban clinic has 2.3 times the odds of referral of a similar patient in a rural clinic with the same baseline. Clinics also differ in their baseline odds, so the population of urban patients mixes clinics with high and low baselines, and so does the population of rural patients. Averaging the probability of referral over this mix pulls the population comparison toward 1, so the population-averaged odds ratio is smaller than the clinic-specific odds ratio. Averaging over clinics is one reason the GEE odds ratio (2.04 here) is smaller than the GLMM odds ratio (2.31), and the weighting used by GEE and chance variation also contribute in this sample. The larger the variation between clinics, the larger the gap. For a continuous outcome with an identity link, averaging does not change the difference, so the two kinds of estimate are the same.

The GLMM and GEE agree that urban clinics have about twice the odds of referral, and the validity of each result rests on assumptions that the next part sets out.

Assumptions and How to Check Them

The GLMM and GEE share some assumptions and differ in others, including the number of clusters each needs. Table 5.7 lists the assumptions of each approach with how to check them, and Box 5.2 describes a step that GEE requires before the model is fitted.

Table 5.7. Assumptions of the logistic GLMM and GEE, with how to check them.

AssumptionGLMMGEEHow to check it
Binary outcome with enough eventsRequiredRequiredtable(outcome); about 10 events per predictor (246 events here).
Independent clustersRequiredRequiredJudged from the study design.
Enough clustersAbout 20 to 30About 30 to 40 or morenlevels(clinic_id); 30 clinics here, at the lower limit for GEE.
Normal cluster effects (log-odds scale)AssumedNot assumedQ-Q plot of ranef(model).
Model fitted without problemsRead any convergence warningsCorrect id, data sorted by clusterWarnings in the console; the "Number of clusters" line in the GEE output.

⚠ Box 5.2: Sort the data before using geeglm()

geeglm() treats consecutive rows with the same id as one cluster. If the rows of a clinic are scattered through the file, the function splits that clinic into several clusters and the standard errors are wrong without any warning. Sorting the data by the cluster identifier (and by time for repeated measures) before fitting avoids the problem. The line Number of clusters: 30 in the output confirms that the clinics were read correctly.

With 30 clinics, the referral data meet the GLMM requirement for the number of clusters and sit at the lower limit for GEE. The checks alone therefore leave the choice between the two models open, and the next part turns to the research question.

Choosing Between a GLMM and GEE

The choice depends mainly on the research question (Subramanian & O'Malley, 2010). A GLMM suits questions about individual clusters, estimates how much the clusters vary, and works with fewer clusters. GEE suits questions about the population as a whole, makes no assumption about the distribution of cluster effects, and needs more clusters for reliable standard errors (Hubbard et al., 2010). In practice, many epidemiological papers report one approach as the main analysis and the other as a sensitivity analysis. When they agree, as they do for urban clinics here, the conclusion is stronger.

This section has fitted a GLMM and GEE to a binary outcome and shown that both widen the standard error for the clinic-level predictor while reaching the same conclusion for urban clinics. Activity 5.3 fits the three models in Table 5.6 in R and compares them, and the knowledge check that follows tests the main ideas of the section. Section 4 applies a linear mixed model and GEE to repeated measurements of the same people.

R Activity 5.3: a GLMM and GEE for specialist referral

This activity fits an ordinary logistic regression, a GLMM and a GEE model to the referral outcome and compares them.

library(lme4); library(geepack)   # glmer() for the GLMM, geeglm() for GEE
clinics <- read.csv("phaa_clinics.csv")   # file must be in the working directory
clinics$clinic_id <- factor(clinics$clinic_id)
clinics$smoker <- factor(clinics$smoker, levels = c("No", "Yes"))
clinics$clinic_urban <- factor(clinics$clinic_urban, levels = c("rural", "urban"))
table(clinics$referred)                  # 1 = referred to a specialist
naive <- glm(referred ~ age + smoker + clinic_urban, family = binomial, data = clinics)
glmm  <- glmer(referred ~ age + smoker + clinic_urban + (1 | clinic_id),
               family = binomial, data = clinics)
summary(glmm)
Console output
0 1 720 246 Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [glmerMod] Family: binomial ( logit ) Formula: referred ~ age + smoker + clinic_urban + (1 | clinic_id) Data: clinics AIC BIC logLik deviance df.resid 1019.6 1044.0 -504.8 1009.6 961 Scaled residuals: Min 1Q Median 3Q Max -1.7792 -0.6038 -0.4221 0.7081 3.6911 Random effects: Groups Name Variance Std.Dev. clinic_id (Intercept) 0.3023 0.5498 Number of obs: 966, groups: clinic_id, 30 Fixed effects: Estimate Std. Error z value Pr(>|z|) (Intercept) -3.943363 0.425525 -9.267 < 2e-16 *** age 0.040244 0.006203 6.488 8.7e-11 *** smokerYes 0.511383 0.193869 2.638 0.00835 ** clinic_urbanurban 0.835563 0.278530 3.000 0.00270 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) age smkrYs age -0.833 smokerYes -0.135 0.036 clnc_rbnrbn -0.498 0.056 0.025

Reading the output. The table shows 246 referrals among 966 patients. In the GLMM, the Random effects block gives a clinic variance of 0.30 on the log-odds scale (standard deviation 0.55), and the Fixed effects block is read like a logistic regression: each coefficient is a log odds ratio for patients in the same clinic or, for urban clinic, in clinics with the same baseline. The urban-clinic coefficient is 0.836 (SE 0.279, p = 0.003), which is an odds ratio of e0.836 = 2.31.

clinics <- clinics[order(clinics$clinic_id), ]   # GEE needs each clinic's rows together
gee <- geeglm(referred ~ age + smoker + clinic_urban, id = clinic_id,
              family = binomial, corstr = "exchangeable", data = clinics)
summary(gee)
options(digits = 7)   # summary() of a GEE model lowers the printed digits; this restores the default
Console output
Call: geeglm(formula = referred ~ age + smoker + clinic_urban, family = binomial, data = clinics, id = clinic_id, corstr = "exchangeable") Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) -3.669981 0.369228 98.796 < 2e-16 *** age 0.038053 0.004334 77.081 < 2e-16 *** smokerYes 0.492949 0.185834 7.036 0.00799 ** clinic_urbanurban 0.711418 0.288681 6.073 0.01373 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation structure = exchangeable Estimated Scale Parameters: Estimate Std.err (Intercept) 0.9896 0.1111 Link = identity Estimated Correlation Parameters: Estimate Std.err alpha 0.04219 0.0157 Number of clusters: 30 Maximum cluster size: 45

The GEE output gives population-averaged coefficients with corrected (sandwich) standard errors in the Std.err column and Wald tests in the next two columns. The urban-clinic coefficient is 0.711 (SE 0.289, p = 0.014), an odds ratio of 2.04. The estimated exchangeable correlation (alpha) is 0.04, and the output confirms 30 clusters. The Estimated Scale Parameters block can be ignored.

# Odds ratios and standard errors for urban clinics from the three models
se <- function(fit) summary(fit)$coefficients["clinic_urbanurban", 2]
b  <- c(naive = coef(naive)[["clinic_urbanurban"]], glmm = fixef(glmm)[["clinic_urbanurban"]],
        gee = coef(gee)[["clinic_urbanurban"]])
round(cbind(OR = exp(b), SE = c(se(naive), se(glmm), se(gee))), 3)
Console output
OR SE naive 2.165 0.168 glmm 2.306 0.279 gee 2.037 0.289

The comparison table shows that the odds ratio for urban clinics is about 2 in all three models, while its standard error rises from 0.168 when clinics are ignored to 0.279 (GLMM) and 0.289 (GEE).

R Reflect on what you just ran

Use the questions below to interpret the output you produced. Look at your console output before answering.

1. Report the odds ratio for urban clinics from the GLMM and from GEE, and explain the difference between what the two odds ratios describe.

Model answerThe GLMM odds ratio is 2.31 and the GEE odds ratio is 2.04. The GLMM odds ratio is clinic-specific: it compares the odds of referral for patients of the same age and smoking status in an urban and a rural clinic with the same baseline. The GEE odds ratio is population-averaged: it compares the proportion referred among all urban-clinic patients with the proportion among all rural-clinic patients of the same age and smoking status, expressed as an odds ratio. The GLMM value is further from 1, as is usual when clusters vary.

2. Compare the standard errors for clinic_urbanurban and for smokerYes across the three models, and explain why one changes much more than the other.

Model answerFor urban clinic, the standard error is 0.168 in the ordinary logistic regression, 0.279 in the GLMM and 0.289 in GEE. For smoking, it is 0.188, 0.194 and 0.186. Urban clinic is a clinic-level predictor, so its information comes from only 30 clinics, and allowing for clustering increases its standard error by about 70%. Smoking varies among patients within every clinic, so its comparison can be made within clinics and its standard error hardly changes.

3. A provincial planner asks whether patients of urban clinics are more likely to be referred across the province. Which model's result would you give, and what check would you mention about it?

Model answerThe GEE result answers a question about the whole population: across the province, patients of urban clinics have about twice the odds of referral of patients of rural clinics (OR 2.04, p = 0.014), adjusting for age and smoking. I would mention that GEE standard errors need enough clusters, usually about 30 to 40 or more, and that with 30 clinics this analysis is at the lower limit, so the GLMM result (OR 2.31, p = 0.003), which reaches the same conclusion, is a useful sensitivity check.
Saved.
Knowledge check: this section

1. Which R function fits a logistic regression with a random intercept for each clinic?

glmer() from lme4 fits generalized linear mixed models; with family = binomial and a term such as (1 | clinic_id) it is a logistic regression with a random intercept.

2. What does GEE do to account for clustering?

GEE fits the regression without cluster effects, which gives population-averaged estimates, and calculates sandwich standard errors from the variation between clusters.

3. A GLMM gives an odds ratio of 2.3 and GEE gives 2.0 for the same predictor. What is the most likely explanation?

For binary outcomes, clinic-specific (conditional) odds ratios are usually further from 1 than population-averaged (marginal) ones, because averaging over clusters with different baselines pulls the comparison toward 1.

4. Why do GEE standard errors need a reasonably large number of clusters?

The corrected (sandwich) standard errors are calculated from the clusters, so with few clusters (fewer than about 30 to 40) they can be too small.

5. Which working correlation does this lesson use for patients clustered in clinics?

Exchangeable assumes every pair of patients in a clinic is equally correlated, which suits people grouped in clinics, schools or households.

✎ Reflection

This section compared two ways of allowing for clustering with a binary outcome. A generalized linear mixed model (GLMM, fitted with glmer()) adds a random intercept for each cluster to a logistic regression and gives clinic-specific (conditional) odds ratios. Generalized estimating equations (GEE, fitted with geeglm()) fit a logistic regression with no cluster effects, calculate standard errors that allow for the correlation within clusters, and give population-averaged (marginal) odds ratios; they need more clusters, about 30 to 40 or more, and data sorted by cluster. Both methods widen the standard errors of cluster-level predictors compared with an ordinary logistic regression. A research team studies whether a smoking-cessation counselling program offered by some pharmacies is associated with quitting, using records of 2,400 customers from 60 pharmacies (30 with the program and 30 without). The team wants to tell the provincial health ministry how quit rates would differ if every pharmacy offered the program. Explain which method you would use as the main analysis and why, what the program variable is (cluster level or person level), what an ordinary logistic regression would get wrong, and which other method you would report as a sensitivity analysis.

Model answerThe program is a pharmacy-level (cluster-level) variable, because every customer of a pharmacy either has access to it or does not, so the information about it comes from the 60 pharmacies. The ministry's question is about the population of customers across the province, so GEE, which gives a population-averaged odds ratio, is the natural main analysis: geeglm(quit ~ program + age + ..., id = pharmacy_id, family = binomial, corstr = "exchangeable"), with the data sorted by pharmacy. With 60 pharmacies there are enough clusters for the corrected standard errors. An ordinary logistic regression would treat the 2,400 customers as independent and would give a standard error for the program that is too small, so the confidence interval would be too narrow and the program could look effective when the evidence is weak. A GLMM with a random intercept for pharmacy would be a useful sensitivity analysis; its odds ratio would be pharmacy-specific (cluster-specific) and probably a little further from 1, and it would also show how much quit rates vary between pharmacies.
✓ Reflection saved!
● Complete the quiz and reflection to continue.
Section 4 of 4

Repeated Measures Over Time

⏱ Estimated time: 45 minutes
Lesson 5 · Section 4

Repeated Measures Over Time

When the same people are measured several times, each person becomes a cluster of visits.

Running example

A wellness trial with four visits

200
adults, 100 per arm
4
visits: 0, 6, 12 and 18 months
2
outcomes: blood pressure and adherence

Each person is a cluster of visits, so the visits from one person are not independent.

Data layout

Long format: one row per visit

idarmvisit (months)sbp_mmhgadherent
R0001control01320
R0001control61310
R0001control121410
R0001control181260
R0002control01271

Mixed models and GEE need long format, with an id column that identifies the person (the cluster).

Look first

Individual lines and arm averages

Panel A: blood pressure over four visits for 24 participants, one line each. Panel B: mean blood pressure by arm, falling from 130.4 to 127.0 in the control arm and from 128.9 to 123.1 in the intervention arm.
Each person stays in their own band (A), and the intervention arm falls more steeply (B).
Missed visits

Mixed models use every visit that was attended

96
of 800 measurements missing (12%)
121
of 200 people attended all four visits
704
measurements used by the mixed model

The results are valid if missing a visit depends only on information in the model (missing at random).

The model

Does blood pressure fall faster in the intervention arm?

lmer(sbp_mmhg ~ arm * visit + (1 | id), data = visits)

TermWhat it measures
(1 | id)A baseline level for each person (random intercept)
arminterventionThe difference between the arms at the start (month 0)
visitThe change per month in the control arm
armintervention:visitThe extra change per month in the intervention arm
Reading the output

An extra fall of about 2.7 mmHg over 18 months

TermEstimateSEp-value
Intervention arm at month 0−1.32 mmHg1.070.22
Change per month, control arm−0.160 mmHg0.043< 0.001
Extra change per month, intervention arm−0.149 mmHg0.0610.014
−2.7
mmHg extra over 18 months (18 × −0.149)
0.56
ICC: how alike each person's visits are
A repeated binary outcome

Adherence over time

Proportion adherent at months 0, 6, 12 and 18. Control: 0.43, 0.41, 0.40, 0.37. Intervention: 0.42, 0.55, 0.69, 0.73.
Adherence rises in the intervention arm and drifts down in the control arm, and the model is fitted with geeglm(adherent ~ arm * visit, id = id, family = binomial, corstr = "exchangeable").
Reading the output

Adherence rises faster in the intervention arm

TermOdds ratio95% CI
Intervention arm at month 01.000.60 to 1.69
Change per month, control arm0.990.96 to 1.02
Extra change per month, intervention arm1.0931.041 to 1.147

Over one six-month interval, 1.0936 ≈ 1.7. The working correlation between a person's visits is 0.07, and there are 200 clusters.

Assumptions and checks

What to check in a repeated-measures analysis

CheckHow
Long format with a person identifierhead(data); sort by id and visit for GEE
Time coded and shaped correctlyPlot the averages by visit; straight lines support a linear trend
Residuals for the mixed modelThe residual and Q-Q plots from Section 2
Missed visitsCount them by visit and compare people who missed visits with those who did not
Lesson summary

Data structure, model and main check

SituationModel and R functionMain check
Any clustered dataICC and design effectHow strong is the clustering, and how much information is lost?
Continuous outcomeLinear mixed model, lmer()Residual plots, Q-Q of cluster effects, number of clusters
Binary outcome, cluster-specificGLMM, glmer()Events, convergence warnings
Binary outcome, population-averagedGEE, geeglm()About 30 to 40 clusters or more, sorted data
Repeated measuresThe same models with the person as the clusterLong format, time trend, missed visits

Introduction and Overview

This section applies the methods of the lesson to repeated measures, in which the same people are measured at several times. The example is phaa_repeated.csv, a simulated wellness trial in which 200 adults were assigned to a control arm or an intervention arm (100 each) and were due to attend visits at 0, 6, 12 and 18 months. At each visit the study recorded systolic blood pressure (sbp_mmhg) and whether the person was adhering to the program (adherent = 1). The section covers the layout of repeated-measures data, missed visits, a linear mixed model for the change in blood pressure, and GEE for the repeated binary outcome.

Learning Objectives

  • Recognise repeated-measures data and explain why each person is a cluster of observations.
  • Distinguish long from wide format and identify the variables a repeated-measures model needs.
  • Explain how a mixed model handles missed visits and what the missing-at-random condition means.
  • Fit and interpret a linear mixed model with an interaction between study arm and time.
  • Fit and interpret a GEE model for a repeated binary outcome, and list the checks for repeated-measures analyses.

Repeated Measures as Clustered Data

When the same person is measured more than once, the measurements share that person's genes, habits, circumstances and baseline health, so they are more alike than measurements from different people. Each person is therefore a cluster of visits, and the methods from the earlier sections apply with the person in place of the clinic. In this trial, the ICC for blood pressure is 0.56: more than half of the variation is a stable difference between people, much more than the clustering of patients in clinics.

The two subsections below describe how the trial data are laid out for these methods and what a plot of the data shows about the dependence.

Long and Wide Format

Repeated-measures data can be stored in wide format, with one row per person and a column for each visit, or in long format, with one row per visit and columns for the person identifier, the visit time and the outcome. Mixed models and GEE in R need long format. The trial file is already in long format, with 800 rows (200 people × 4 visits). Table 5.8 shows the first five rows of the file.

Table 5.8. The first five rows of phaa_repeated.csv in long format, with one row per visit.

idarmvisit (months)sbp_mmhgadherent
R0001control01320
R0001control61310
R0001control121410
R0001control181260
R0002control01271

The id column identifies the person (the cluster), visit records months since the start, and each person appears once for each visit. In wide format, the same person R0001 would be a single row with columns such as sbp_0, sbp_6, sbp_12 and sbp_18. The base R function reshape() or tidyr::pivot_longer() converts wide data to long format.

Looking at the Data

Before any model is fitted, a plot of the measurements over time shows how strongly each person’s visits are linked and what shape the trend takes. Figure 5.7 shows the blood pressure of individual participants at each visit and the average in each arm.

Panel A: systolic blood pressure at four visits for 24 participants, one line each, coloured by arm. Panel B: mean blood pressure by arm at months 0, 6, 12 and 18, falling from 130.4 to 127.0 in the control arm and from 128.9 to 123.1 in the intervention arm.
Figure 5.7. Panel A shows one line for each of 24 randomly chosen participants, and panel B shows the average blood pressure in each arm at each visit. Each person tends to stay in their own band, and the intervention arm falls more steeply than the control arm.

The spaghetti plot (panel A) shows the dependence directly: people start at very different levels and tend to stay there. The averages (panel B) fall in both arms, from 130.4 to 127.0 mmHg in the control arm and from 128.9 to 123.1 mmHg in the intervention arm, and the roughly straight lines suggest that a linear trend over time is reasonable.

The trial data are therefore in the long format that the models need, and Figure 5.7 shows both the dependence within people and a roughly linear trend. Some planned visits were missed, however, and the next part considers how missed visits affect the analysis.

Missed Visits

Of the 800 planned measurements of blood pressure, 96 (12%) are missing, and only 121 of the 200 participants attended all four visits. Dropping everyone with a missed visit (a complete-case analysis) would discard 79 people and could bias the results if the people who missed visits differ from those who attended every one. For a regression model, a complete-case analysis still gives unbiased coefficients when the chance that a record is incomplete depends only on the predictors in the model and not on the outcome, although it loses precision; this exception refines the statement in Lesson 2 that listwise deletion is biased unless the data are MCAR. In repeated-measures data a missed visit often depends on the person’s earlier outcome values, and then the exception does not apply. A mixed model uses every measurement that was taken, here 704 measurements from all 200 people.

Box 5.3 recalls the three missing-data mechanisms from Lesson 2 and ends with a retrieval question, and Box 5.4 restates the missing-at-random condition in plain terms for the trial.

Box 5.3: Recall: Lesson 2, Section 1 (Missing Data)

Lesson 2, Section 1 set out Rubin’s three missing-data mechanisms. Data are missing completely at random (MCAR) when the chance that a value is missing is unrelated to any value, observed or not. They are missing at random (MAR) when that chance depends only on values that were observed, and missing not at random (MNAR) when it depends on the missing value itself. Students who took HSCI 230 first met these mechanisms in Lesson 8, Section 2 (Attrition and Nonresponse Bias) and Lesson 11, Section 2 (Statistical Inference and Model Issues).

Retrieval question. A participant in the wellness trial misses the 12-month visit. Give one reason for the missed visit that would make the missing blood pressure MAR, and one that would make it MNAR.

AnswerThe value is MAR if the chance of missing the visit depends only on recorded information, for example if people whose blood pressure was high at the 6-month visit were more likely to miss the next one. It is MNAR if the chance depends on the unrecorded value itself, for example if the person stayed away because their blood pressure was high at the time of the 12-month visit, in a way that the earlier visits could not predict.

Box 5.4: Missing at random, in plain terms

A mixed model fitted to all available visits gives valid results when the chance of missing a visit depends only on information the model already uses, such as the person's arm or their blood pressure at earlier visits. This condition is called missing at random (MAR). It would fail if, for example, people stayed away from a visit because their blood pressure was high on that day, in a way the earlier visits could not predict. MAR cannot be proved from the data, but it is a much weaker assumption than the complete-case assumption that people who attended every visit are typical of everyone. Comparing people who missed visits with those who did not at their first visit shows which measured characteristics predict a missed visit, so that they can be included in the model. The comparison cannot confirm MAR, because MAR concerns the values that were never recorded.

When the missing-at-random condition holds, a mixed model fitted to all 704 available measurements gives valid results. The next part fits that model to the blood pressure data.

A Mixed Model for Change Over Time

The model for blood pressure includes a random intercept for each person and fixed effects for arm, visit and their interaction. Equation 5.5 writes the model.

Random-intercept model for change over time
\[ \text{SBP}_{it} = \color{#6D28D9}{\beta_0} + \color{#C2410C}{u_i} + \color{#1D4ED8}{\beta_1}\,\text{arm}_i + \color{#1D4ED8}{\beta_2}\,\text{visit}_{t} + \color{#0B7B6B}{\beta_3}\,(\text{arm}_i \times \text{visit}_{t}) + \varepsilon_{it} \]Eq 5.5
Blood pressure for person i at time t equals an overall intercept, a shift for person i, a difference between arms at the start and a change per month in the control arm, plus the extra change per month in the intervention arm (β3), which is the effect of the program on the trend.

Table 5.9 gives the estimates of the fixed effects for arm, visit and their interaction, with an interpretation of each.

Table 5.9. Fixed effects from the random-intercept model for blood pressure over time.

TermEstimateSEp-valueInterpretation
Intervention arm (month 0)−1.32 mmHg1.070.22The arms start at similar levels, as expected after randomization.
Visit (control arm, per month)−0.160 mmHg0.043< 0.001Blood pressure falls by about 0.16 mmHg per month in the control arm.
Arm × visit (per month)−0.149 mmHg0.0610.014Blood pressure falls an extra 0.15 mmHg per month in the intervention arm.

Over the 18 months of the trial, the extra fall in the intervention arm is 18 × −0.149 = −2.7 mmHg. The random effects give a between-person variance of 34.80 and a within-person variance of 27.58, so the ICC is 34.80 ÷ (34.80 + 27.58) = 0.56.

Worked Example 5.4 combines the intercept of the model with the estimates in Table 5.9 to predict blood pressure at month 18 in each arm.

Worked Example 5.4: Predicted Blood Pressure at Month 18

For a typical person in the control arm, the model predicts 130.02 − 0.160 × 18 = 127.1 mmHg at month 18. For a typical person in the intervention arm, it predicts 130.02 − 1.32 + (−0.160 − 0.149) × 18 = 123.1 mmHg. Both predictions are close to the observed averages of 127.0 and 123.1 mmHg.

The accordion below describes two extensions of this model, random slopes for time and a correlation that depends on the time between visits, both of which are beyond the scope of this lesson.

Going further: random slopes for time and correlation over time

The random-intercept model gives each person their own level but the same slope over time within each arm. A random slope, (visit | id), would let each person have their own rate of change. In this trial that model produces a convergence warning, because four visits per person give little information about individual slopes. Visits close together in time are also often more alike than visits far apart, a pattern that models with an autoregressive correlation structure (such as AR(1)) describe. Both extensions appear in published longitudinal studies and are beyond the scope of this lesson.

The mixed model shows an extra fall of about 2.7 mmHg over 18 months in the intervention arm. The trial also recorded a binary outcome at each visit, and the next part models it with GEE.

GEE for a Repeated Binary Outcome

Adherence is recorded as yes or no at each visit. The share adherent falls slightly in the control arm (0.43 at month 0 to 0.37 at month 18) and rises steadily in the intervention arm (0.42 to 0.73). Figure 5.8 shows these proportions at each visit.

Proportion adherent by arm at months 0, 6, 12 and 18. Control: 0.43, 0.41, 0.40, 0.37. Intervention: 0.42, 0.55, 0.69, 0.73.
Figure 5.8. The lines show the proportion of participants adhering to the program at each visit in each arm. Adherence drifts down in the control arm and rises in the intervention arm.

A GEE model with the person as the cluster gives population-averaged odds ratios for the trial, which matches the question of how the program changes adherence across participants. With 200 people there are plenty of clusters for the corrected standard errors. Table 5.10 gives the odds ratios from this model.

Table 5.10. Odds ratios from the GEE model for adherence over time.

TermOdds ratio95% CIInterpretation
Intervention arm (month 0)1.000.60 to 1.69The arms start with the same odds of adherence.
Visit (control arm, per month)0.990.96 to 1.02Adherence in the control arm barely changes.
Arm × visit (per month)1.0931.041 to 1.147The monthly odds ratio is 0.99 × 1.093 = 1.08 in the intervention arm, a rise of about 8% per month.

Over one six-month interval between visits, the interaction corresponds to an odds ratio of 1.0936 ≈ 1.70. The estimated exchangeable correlation between a person's visits is 0.07, much weaker than the correlation for blood pressure, but GEE allows for it in the standard errors.

GEE needs a stronger condition about missed visits than the mixed model. Its results are valid when the chance of missing a visit depends only on the predictors in the GEE model, here arm and visit. If missed visits depend on earlier outcomes, such as earlier blood pressure or adherence, a mixed model, or a weighted form of GEE that is beyond this lesson, is preferred.

GEE therefore shows a steady rise in adherence in the intervention arm, provided that missed visits depend only on arm and visit. The next part collects the checks that apply to both repeated-measures models.

Assumptions and How to Check Them

A repeated-measures analysis needs some of the checks from the earlier sections together with checks of its own, which concern the layout of the data, the coding of time and the missed visits. Table 5.11 lists them for the mixed model and for GEE.

Table 5.11. Assumptions and requirements of a repeated-measures analysis, with how to check them.

Assumption or requirementHow to check it
Data in long format with a person identifierhead(data); for GEE, sort by id and then by visit.
Time coded correctly and the trend shape reasonablePlot the averages by visit; roughly straight lines support a linear trend in visit.
Residual assumptions for the mixed modelResiduals versus fitted values and Q-Q plots, as in Section 2.
Missed visits: missing at random for the mixed model; related only to the model's predictors for GEECount missing values by visit and compare people who missed visits with those who did not. This shows which variables to include; it cannot confirm MAR.
Enough clusters for GEEThe Number of clusters line of the output (200 here).
Independent peopleJudged from the study design (one person per household, for example).

These checks complete the analysis of the wellness trial. The final part sets the four examples of the lesson side by side.

Bringing the Lesson Together

The four sections of this lesson apply one approach. The cluster variable is identified first, the clustering is measured with the ICC, and a model is chosen that allows for it: a linear mixed model for a continuous outcome, a GLMM or GEE for a binary outcome, and the same models with the person as the cluster for repeated measures. Table 5.12 summarises the four examples with the cluster, the model and R function, and the main result of each.

Table 5.12. The four examples of the lesson, with the cluster, model and main result of each.

ExampleClusterModel and R functionMain result
Blood pressure in clinicsClinicLinear mixed model, lmer()Urban minus rural 1.85 mmHg (−1.85 to 5.54); ICC 0.24
Referral in clinicsClinicGLMM, glmer(), and GEE, geeglm()Urban OR 2.31 (GLMM) and 2.04 (GEE)
Blood pressure over timePersonLinear mixed model, lmer()Extra fall of 0.149 mmHg per month in the intervention arm
Adherence over timePersonGEE, geeglm()Interaction OR 1.093 per month (1.041 to 1.147)

Activity 5.4 fits the mixed model for blood pressure over time in R, and Activity 5.5 continues from it with the GEE model for adherence. The knowledge check that follows the two activities tests the main ideas of this section, and the final page of the lesson gathers the key takeaways and the final assessment.

R Activity 5.4: a mixed model for blood pressure over time

This activity fits the linear mixed model for blood pressure in the wellness trial. Download phaa_repeated.csv and save it in your working directory. The next activity continues from this one.

library(lmerTest)   # loads lme4 and adds p-values
visits <- read.csv("phaa_repeated.csv")   # file must be in the working directory
visits$id  <- factor(visits$id)
visits$arm <- factor(visits$arm, levels = c("control", "intervention"))
head(visits, 8)                                      # long format: one row per visit
table(visit = visits$visit, missing = is.na(visits$sbp_mmhg))
round(tapply(visits$sbp_mmhg, list(visits$arm, visits$visit), mean, na.rm = TRUE), 1)
Console output
id arm age female visit sbp_mmhg adherent 1 R0001 control 54 1 0 132 0 2 R0001 control 54 1 6 131 0 3 R0001 control 54 1 12 141 0 4 R0001 control 54 1 18 126 0 5 R0002 control 59 1 0 127 1 6 R0002 control 59 1 6 125 0 7 R0002 control 59 1 12 131 1 8 R0002 control 59 1 18 143 1 missing visit FALSE TRUE 0 180 20 6 176 24 12 180 20 18 168 32 0 6 12 18 control 130.4 128.9 128.1 127.0 intervention 128.9 126.8 125.4 123.1

The first eight rows show the long format: person R0001 has four rows, one per visit, followed by person R0002. The table of missing values shows between 20 and 32 missed measurements at each visit (96 in all), and the averages by arm and visit match Figure 5.7.

lmm_t <- lmer(sbp_mmhg ~ arm * visit + (1 | id), data = visits)
summary(lmm_t)
Console output
Linear mixed model fit by REML. t-tests use Satterthwaite's method [ lmerModLmerTest] Formula: sbp_mmhg ~ arm * visit + (1 | id) Data: visits REML criterion at convergence: 4672.1 Scaled residuals: Min 1Q Median 3Q Max -2.7904 -0.5636 0.0015 0.5875 3.3301 Random effects: Groups Name Variance Std.Dev. id (Intercept) 34.80 5.899 Residual 27.58 5.252 Number of obs: 704, groups: id, 200 Fixed effects: Estimate Std. Error df t value Pr(>|t|) (Intercept) 130.02227 0.75761 335.39051 171.622 < 2e-16 *** armintervention -1.31505 1.06859 332.64253 -1.231 0.219329 visit -0.15981 0.04314 514.38871 -3.704 0.000235 *** armintervention:visit -0.14925 0.06055 514.62646 -2.465 0.014032 * --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation of Fixed Effects: (Intr) armntr visit armintrvntn -0.709 visit -0.502 0.356 armntrvntn: 0.358 -0.498 -0.712

Reading the output. The line Number of obs: 704, groups: id, 200 confirms that the model used all 704 available measurements from all 200 people. In the Fixed effects block, armintervention (−1.32, p = 0.22) is the difference between the arms at month 0, visit (−0.160, p < 0.001) is the change per month in the control arm, and armintervention:visit (−0.149, p = 0.014) is the extra change per month in the intervention arm.

vc <- as.data.frame(VarCorr(lmm_t))
vc$vcov[1] / sum(vc$vcov)                            # how alike one person's visits are
18 * fixef(lmm_t)[["armintervention:visit"]]          # extra change over 18 months
Console output
[1] 0.5578058 [1] -2.686584

The ICC of 0.56 shows that a person's visits are strongly alike, and the extra fall in blood pressure in the intervention arm over 18 months is 2.7 mmHg.

R Reflect on what you just ran

Use the questions below to interpret the output you produced. Look at your console output before answering.

1. Report the estimate, standard error and p-value for armintervention:visit, and explain in one or two sentences what it says about the program.

Model answerThe interaction is −0.149 mmHg per month (SE 0.061, p = 0.014). Blood pressure fell 0.149 mmHg per month faster in the intervention arm than in the control arm, which adds up to an extra fall of about 2.7 mmHg over the 18 months, and the p-value indicates that this difference in trends is unlikely to be due to chance alone.

2. How many measurements did the model use, and why is this better than analysing only the 121 people who attended all four visits? What condition must hold for the result to be valid?

Model answerThe model used 704 measurements from all 200 people. Analysing only the 121 people with complete data would discard 79 people and could bias the result if people who missed visits differ from those who attended all of them. The mixed model result is valid if the chance of missing a visit depends only on information in the model, such as arm and earlier blood pressure readings (missing at random).

3. Report the ICC from lmm_t and compare it with the ICC for patients in clinics (0.175). What does the difference mean?

Model answerThe ICC is 0.56, compared with 0.175 for patients in clinics. More than half of the variation in blood pressure in the trial is a stable difference between people, so the four visits of one person are much more alike than the patients of one clinic. Ignoring this dependence would be a serious error in a repeated-measures analysis.
Saved.
R Activity 5.5: GEE for adherence over time

This activity continues from the previous one and uses the visits data frame created there. It fits a GEE model for the repeated binary outcome, adherence.

library(geepack)
adh <- visits[!is.na(visits$adherent), ]             # visits with adherence recorded
adh <- adh[order(adh$id, adh$visit), ]               # each person's rows together
round(tapply(adh$adherent, list(adh$arm, adh$visit), mean), 2)   # share adherent
gee_a <- geeglm(adherent ~ arm * visit, id = id, family = binomial,
                corstr = "exchangeable", data = adh)
summary(gee_a)
options(digits = 7)   # restore the default number of printed digits
Console output
0 6 12 18 control 0.43 0.41 0.40 0.37 intervention 0.42 0.55 0.69 0.73 Call: geeglm(formula = adherent ~ arm * visit, family = binomial, data = adh, id = id, corstr = "exchangeable") Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) -0.278541 0.198606 1.967 0.160773 armintervention 0.005248 0.264926 0.000 0.984197 visit -0.012340 0.016809 0.539 0.462865 armintervention:visit 0.088695 0.024784 12.807 0.000345 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Correlation structure = exchangeable Estimated Scale Parameters: Estimate Std.err (Intercept) 1 0.0292 Link = identity Estimated Correlation Parameters: Estimate Std.err alpha 0.06572 0.03557 Number of clusters: 200 Maximum cluster size: 4

Reading the output. The table of shares adherent matches Figure 5.8. In the GEE output, armintervention (0.005) shows no difference between the arms at month 0, visit (−0.012) shows little change in the control arm, and armintervention:visit (0.089, p < 0.001) shows that adherence rises faster in the intervention arm. The estimated correlation between a person's visits (alpha) is 0.07, and there are 200 clusters.

round(exp(cbind(OR = coef(gee_a), confint.default(gee_a))), 3)   # odds ratios, 95% CIs
Console output
OR 2.5 % 97.5 % (Intercept) 0.757 0.513 1.117 armintervention 1.005 0.598 1.690 visit 0.988 0.956 1.021 armintervention:visit 1.093 1.041 1.147

The odds ratio for the interaction is 1.093 per month (95% CI 1.041 to 1.147). The function confint.default() gives Wald confidence intervals from the corrected standard errors.

R Reflect on what you just ran

Use the questions below to interpret the output you produced. Look at your console output before answering.

1. Report the odds ratio and 95% confidence interval for armintervention:visit, and explain what it means for adherence in the two arms.

Model answerThe odds ratio is 1.093 per month (95% CI 1.041 to 1.147). In the intervention arm the odds of adherence rise by about 8% per month (0.99 × 1.093 = 1.08), while in the control arm they fall by about 1% per month (OR 0.99). Over one six-month interval this is an odds ratio of 1.0936 ≈ 1.70, so the program is associated with a steady rise in adherence.

2. Why is GEE, with id = id, a reasonable choice for this outcome, and why must the data be sorted before fitting it?

Model answerEach person contributes up to four yes-or-no adherence measurements, so the visits are clustered within people, and GEE allows for that clustering in the standard errors while giving population-averaged odds ratios that answer a question about the trial population. With 200 people there are enough clusters. geeglm() treats consecutive rows with the same id as one cluster, so the rows must be sorted by person (and by visit) or a person's visits could be split into several clusters, giving wrong standard errors.

3. The estimated correlation (alpha) between a person's adherence measurements is 0.07, while the ICC for blood pressure was 0.56. Would it be acceptable to analyse adherence with an ordinary glm() because the correlation is small? Explain.

Model answerIt would not be the right choice. Even a small correlation among repeated measurements breaks the independence assumption, and the study design guarantees that the measurements are grouped within people. GEE allows for the correlation at little cost, and its standard errors remain valid whether the correlation is small or large. Reporting the GEE result, perhaps with the ordinary glm() result as a comparison, is the appropriate approach.
Saved.
Knowledge check: this section

1. In repeated-measures data, what is the cluster?

Each person contributes several measurements that share their characteristics, so the person is the cluster and the visits are the observations within it.

2. Which data layout do lmer() and geeglm() need for repeated measures?

Long format, with one row per visit and columns for the person identifier, the time and the outcome, is required by mixed models and GEE in R.

3. In the model lmer(sbp_mmhg ~ arm * visit + (1 | id)), which term tests whether blood pressure changes at a different rate in the two arms?

The interaction term is the difference in the slope over time between the arms, which answers whether the program changes the trend.

4. Why can a mixed model use people who missed some visits?

A mixed model fitted to all available measurements gives valid results when the data are missing at random, meaning that missingness depends only on observed information such as arm or earlier values.

5. A GEE interaction odds ratio is 1.093 per month. What is the corresponding odds ratio over a six-month interval?

Odds ratios multiply over time: 1.0936 ≈ 1.70.

✎ Reflection

This section treated repeated measures as clustered data in which each person is a cluster of visits. The data must be in long format (one row per visit) with a person identifier. A linear mixed model such as lmer(outcome ~ arm * visit + (1 | id)) gives each person a random intercept, its arm:visit interaction tests whether the outcome changes at a different rate in the two arms, and it uses every visit that was attended, which is valid if missed visits are missing at random (missingness depends only on information in the model). GEE with id as the cluster gives population-averaged odds ratios for a repeated binary outcome. A study follows 300 older adults, half of whom join a weekly walking group, and records their walking speed (in metres per second) and whether they had a fall in the previous three months, every three months for one year (five visits). About 15% of visits are missed, more often among people who were frailer at their first visit. Describe the data layout you would need, the model you would fit for walking speed and the term that answers the main question, the model you would fit for falls, and what the pattern of missed visits means for the analysis.

Model answerThe data need to be in long format, with one row per person per visit (up to five rows per person) and columns for the person identifier, group (walking group or not), months since the start, walking speed and fall (yes or no). For walking speed, a continuous outcome, I would fit lmer(speed ~ group * month + frailty + (1 | id)), which gives each person their own baseline speed; the group:month interaction answers the main question, whether walking speed changes at a different rate in the walking group, and I would check the residual plots and whether the average speeds by visit look roughly linear. For falls, a repeated binary outcome, I would fit GEE with geeglm(fall ~ group * month + frailty, id = id, family = binomial, corstr = "exchangeable") after sorting the data by person and visit; with 300 people there are plenty of clusters, and the interaction odds ratio would describe how the odds of a fall change over time in the walking group compared with the other group. Because frailer people at the first visit miss more visits, a complete-case analysis would keep a healthier group and could bias the results, while the mixed model uses all visits and is valid if missingness depends on observed information such as baseline frailty, so I would include baseline frailty in the models (the first-visit walking speed is already part of the outcome data that the mixed model uses) and compare people who missed visits with those who did not. For the GEE model of falls, baseline frailty also needs to be a predictor, because GEE is valid only when missed visits depend on the predictors in the model.
✓ Reflection saved!
● Complete the quiz and reflection to continue.
Final Assessment

Lesson 5: Final Assessment

15 questions • 100% required to pass

Bringing It All Together

This lesson addressed the assumption of independence that every model in the previous lessons shares. Section 1 showed that observations are dependent when they share a cluster, such as patients in the same clinic or visits by the same person. The intracluster correlation coefficient measured the clustering in the clinic data (0.175 for systolic blood pressure), the design effect of 6.45 showed that the 966 patients carried about as much information as 150 independent patients, and an ordinary regression that ignored the clinics gave a standard error for urban clinics of 0.74 mmHg, compared with 1.96 mmHg from a model that allowed for them.

Section 2 fitted a linear mixed model with a random intercept for each clinic using lmer(). Its fixed effects were read like the coefficients of a linear regression, its random effects gave a clinic variance and a residual variance whose ratio was the adjusted ICC of 0.24, and it was checked with residual plots, a Q-Q plot of the clinic effects, the number of clinics and a test for a singular fit. Section 3 turned to a binary outcome, specialist referral, and compared a logistic mixed model fitted with glmer(), which gives cluster-specific odds ratios, with GEE fitted with geeglm(), which gives population-averaged odds ratios with sandwich standard errors. Section 4 treated each person in a wellness trial as a cluster of visits, arranged the data in long format, and used a mixed model for blood pressure and GEE for adherence, each with an arm-by-visit interaction that answered the main question of the trial and each using every visit that was attended.

Across the lesson, the same steps were followed each time: the cluster variable was identified, the clustering was measured, a model that allows for it was chosen to suit the outcome type and the research question, and the model's assumptions were checked before the results were reported. The final assessment asks for these steps to be applied to new studies.

Key Takeaways from this lesson

  • Observations are dependent when they share a cluster, and the cluster variable (clinic, school, neighbourhood or person) is identified from the study design before any model is fitted.
  • The ICC is the share of the outcome's variation that lies between clusters, and the design effect, 1 + (average cluster size − 1) × ICC, shows how much information the clustering removes.
  • Ignoring clustering usually makes standard errors too small, and the problem is largest for cluster-level predictors such as an urban clinic or a school-wide program; for person-level predictors, allowing for clustering changes the standard error little and can make it smaller.
  • A linear mixed model with a random intercept, lmer(y ~ x + (1 | cluster)), gives each cluster its own baseline, and its fixed effects are read like ordinary regression coefficients.
  • A mixed model is checked with residual plots, a Q-Q plot of the cluster effects, enough clusters (about 20 to 30 or more) and a test for a singular fit.
  • For a binary outcome, a GLMM (glmer()) gives cluster-specific odds ratios and GEE (geeglm()) gives population-averaged odds ratios, which are usually a little closer to 1.
  • Repeated measures are analysed in long format, the group-by-time interaction tests whether the groups change at different rates, and a mixed model that uses every attended visit is valid when missed visits are missing at random.

The final assessment covers all four sections. All 15 questions must be answered correctly (100%), and the final reflection completed, to finish the lesson.

Reflection

A study evaluates a school-based physical activity program in 40 elementary schools, 20 of which were randomly assigned to run the program. About 50 students are surveyed in each school (2,000 students in all). The outcomes are the minutes of physical activity each student reports per day and whether the student meets the guideline of 60 minutes per day (yes or no). The predictors are the program (a school-level variable), the student's age and gender. A model with no predictors gives an ICC of 0.04 for minutes of activity. This lesson showed that observations in the same cluster are dependent, that the design effect is 1 + (average cluster size − 1) × ICC and the effective sample size is the sample size divided by the design effect, that a linear mixed model with a random intercept is fitted with lmer(y ~ x + (1 | cluster)) and checked with residual plots, a Q-Q plot of the cluster effects, the number of clusters and isSingular(), and that a binary outcome can be modelled with a GLMM (glmer(), cluster-specific odds ratios) or GEE (geeglm(), population-averaged odds ratios). Identify the cluster, calculate the design effect and the effective sample size, and explain what would happen to the standard error for the program if the schools were ignored. Then state the model and R code you would use for each outcome, what the program coefficient or odds ratio would mean, and the main checks you would carry out.

Model answerThe cluster is the school, because students in the same school share teachers, facilities and the program itself. With about 50 students per school and an ICC of 0.04, the design effect is 1 + 49 × 0.04 = 2.96, so the 2,000 students carry about as much information as 2,000 / 2.96 ≈ 676 independent students. The program is a school-level predictor, so an ordinary regression that ignored the schools would treat 2,000 students as independent and give a standard error for the program that is far too small, with a confidence interval that is too narrow and a p-value that is too small. For minutes of activity, a continuous outcome, I would fit a linear mixed model with a random intercept for each school: lmer(minutes ~ program + age + gender + (1 | school_id), data = students) after loading lmerTest. The program coefficient would be the difference in average minutes of activity between students of the same age and gender in program and non-program schools; for example, a coefficient of 8 would mean about 8 more minutes per day. I would check the residuals versus fitted values and the Q-Q plot of the residuals, a Q-Q plot of the school effects from ranef(), that 40 schools is enough (it is above the rough minimum of 20 to 30), and that isSingular() is FALSE. For meeting the guideline, a binary outcome, I would fit either a GLMM, glmer(meets ~ program + age + gender + (1 | school_id), family = binomial), whose odds ratio compares a student in a program school with a student of the same age and gender in a non-program school with the same baseline, or GEE, geeglm(meets ~ program + age + gender, id = school_id, family = binomial, corstr = "exchangeable") after sorting by school, whose odds ratio compares all students in program schools with all students in non-program schools. Because the question is whether the program works across schools, the population-averaged GEE odds ratio answers it directly, and 40 schools is enough for the sandwich standard errors; I would expect the GLMM odds ratio to be a little further from 1. For both models I would check that there are enough students meeting the guideline for the number of predictors.

Minimum 20 characters required.

✓ Reflection saved

Final Knowledge Assessment

Final assessment: the 15 questions

1. Which of these studies produces clustered data?

Students in the same school share teachers, facilities and policies, so the school is a cluster and the students' outcomes are not independent.

2. In a study in which each participant is measured at four visits, what is the cluster?

Each participant forms a cluster of visits, because measurements on the same person tend to be more alike than measurements on different people.

3. A mixed model with no predictors gives a between-clinic variance of 10 and a within-clinic variance of 90. What is the ICC?

The ICC is the between-cluster variance divided by the total variance: 10 / (10 + 90) = 0.10.

4. Clusters contain 41 people on average and the ICC is 0.05. What is the design effect?

The design effect is 1 + (41 − 1) × 0.05 = 1 + 2.0 = 3.0.

5. A clustered study of 1,200 people has a design effect of 4. What is its effective sample size?

The effective sample size is the sample size divided by the design effect: 1,200 / 4 = 300.

6. A study of patients in 30 clinics ignores the clinics and treats every patient as independent. What usually happens to the standard error for a clinic-level predictor such as an urban location?

Ignoring clustering overstates the amount of independent information, so the standard error of a cluster-level predictor is too small and its p-value is too small.

7. In lmer(sbp ~ age + smoker + (1 | clinic_id)), what does (1 | clinic_id) do?

The term (1 | clinic_id) is a random intercept: each clinic is shifted up or down from the overall average, and the model estimates the variance of these shifts.

8. In a linear mixed model of blood pressure, the fixed effect for female is −2.8 mmHg. What does it mean?

A fixed effect is read like an ordinary regression coefficient, and in a mixed model it compares people who share a cluster, so the clinic shift cancels out of the comparison.

9. Which output would show that a mixed model's random-effects part is more complex than the data can support?

A singular fit means that a random-effect variance has been estimated as zero (or a correlation as ±1), which usually means the random-effects structure should be simplified.

10. Which plot checks the assumption that the clinic effects of a mixed model are normally distributed?

The clinic effects are estimated by ranef(), and a Q-Q plot of them with points near the line supports the normality assumption.

11. What does GEE do to account for clustering in a logistic regression?

GEE fits a logistic regression without cluster effects, weights the data using a working correlation, and calculates sandwich standard errors from the variation between clusters.

12. A GLMM gives an odds ratio of 2.3 for urban clinics and GEE gives 2.0. What is the most likely explanation?

For a binary outcome the cluster-specific (GLMM) and population-averaged (GEE) odds ratios answer different questions, and the population-averaged odds ratio is usually closer to 1.

13. A GEE analysis has only 8 clusters. What is the main concern?

Sandwich standard errors are estimated from the clusters, so with few clusters they are unstable and tend to be too small; at least 30 to 40 clusters is a common guide.

14. A dataset has one row per person with the columns sbp_0, sbp_6, sbp_12 and sbp_18. What must be done before fitting lmer()?

Mixed models and GEE need long format, in which each visit is a row and the person identifier links the rows that belong to the same person.

15. A trial fits lmer(sbp ~ arm * visit + (1 | id)). Some participants missed visits, and missingness depends on their earlier blood pressure. Which statement is correct?

A mixed model uses all the visits that were attended, and its results are valid when missingness depends only on observed information. The term that tests whether the arms change at different rates is the arm:visit interaction.

✦ Before submitting: pass every section knowledge check (100%) and complete every reflection.

🏆 Congratulations!

This lesson completes the regression sequence of the course. The models from Lessons 3 and 4 can now be applied when people are grouped in clinics, schools or neighbourhoods, or are measured more than once.

You have successfully completed this lesson: Modelling Dependent Data.

Your responses have been downloaded automatically.