Building models, checking assumptions, and clustered data

In 14.1 you fitted generalized linear models by maximum likelihood and read a likelihood ratio test. That gave you the machinery. This lesson is about judgement: which model to fit, how to tell whether it fits, and what to do when the data refuse to behave. The largest failures in applied medical statistics are rarely arithmetic. They come from forcing one model onto data whose design it does not match.

Decide what the model is for

Before you name a single variable, settle one question: are you estimating an effect or predicting an outcome? The answer changes how you build the model.

• Estimation. You want the effect of one exposure, say metformin adherence on HbA1c, with everything else held fair. You include confounders, keep the model small enough to interpret, and you report a coefficient with a confidence interval.
• Prediction. You want an accurate forecast for a new patient, say the risk of readmission. You include any variable that improves prediction, and you judge the model on how well it predicts data it has not seen, not on the meaning of one coefficient.

The same dataset gives two different model-building strategies depending on which goal you pick. Choosing variables to estimate an effect is mostly about confounding. Choosing variables to predict is mostly about out-of-sample accuracy.

Confounders versus effect modifiers

For estimation, two roles get confused. A confounder is associated with both the exposure and the outcome and distorts the crude association. Age confounds the link between exercise and blood pressure: older patients exercise less and have higher pressure. You adjust for a confounder by putting it in the model, which removes its distortion.

An effect modifier is different. It changes the size of the exposure effect across its own levels. If a drug lowers HbA1c more in newly diagnosed patients than in long-standing ones, baseline duration modifies the treatment effect. You do not adjust an effect modifier away. You report the exposure effect within each level, or you add an interaction term and read the two effects separately. Treating a modifier as a plain confounder hides a finding that matters for who should get the treatment.

The assumptions, and how to check them

The linear model from 5.2 rests on four assumptions. The generalized linear model keeps two of them and replaces the other two with a chosen distribution and link. Check them in order.

1. Correct mean structure (linearity). The linear predictor describes the mean (after the link). Plot residuals against fitted values. A flat band is fine. A curve means a term is mis-specified, often a continuous variable that needs a quadratic or a log.
2. Independence. Each observation carries its own information. This one comes from the design, not a plot. Repeated measures or grouped sampling break it.
3. Constant variance. For a linear model the residual spread should not change with the fitted value. A fan or funnel shape signals non-constant variance. For a GLM the variance is tied to the mean by the family, so you check for over-dispersion instead.
4. Distribution of residuals. A linear model assumes normal residuals: read a normal quantile plot. A GLM sets the distribution through the family (binomial for binary, Poisson for counts), so here you check that the family and link are the right choice.

When an assumption fails: three moves

A broken assumption does not mean you abandon regression. You have three responses, in rough order of how much they change the model.

1. Transform the outcome. When the outcome is skewed or its variance grows with its mean, a log transform often restores linearity and constant variance at the same time. Costs and lengths of stay are the usual candidates.
2. Change the GLM family or link. A count that piles up at zero belongs in a Poisson or negative-binomial model, not a linear one. A binary outcome belongs in logistic regression. Matching the family to the outcome type usually fixes the variance assumption by construction.
3. Keep the model, repair the standard errors. When the mean model is right but the variance or independence assumption is off, leave the coefficients alone and compute sandwich standard errors (also called cluster-robust standard errors). These estimate the spread from the residuals in the data rather than from the assumed distribution. The estimate does not move; the standard error, and so the confidence interval, does.

Worked example 1: sandwich standard errors for clustered data

A diabetes-care study covers 192 patients across 16 Klinik Kesihatan sites. Eight clinics adopted a structured education programme; eight did not. The programme is a clinic-level exposure: every patient in a clinic shares the same value. The outcome is HbA1c. An ordinary regression treats all 192 patients as independent and reports a small standard error. But patients in one clinic resemble each other, so the real information is closer to 16 clinics than 192 patients. Cluster-robust standard errors correct for that.

set.seed(8)
n_clinic <- 16
per      <- 12
clinic   <- rep(1:n_clinic, each = per)
program  <- rbinom(n_clinic, 1, 0.5)[clinic]
u        <- rnorm(n_clinic, 0, 0.8)[clinic]
hba1c    <- 8 - 0.4 * program + u + rnorm(n_clinic * per, 0, 0.7)
dia      <- data.frame(hba1c, program, clinic = factor(clinic))

m <- # fit lm(hba1c ~ program)

m <- lm(hba1c ~ program, data = dia)
round(summary(m)$coefficients["program", ], 4)
# Estimate Std. Error t value Pr(>|t|)
# the model-based s.e. ignores clustering and is too small;
# for the cluster-robust s.e. in R use sandwich::vcovCL(m, cluster = dia$clinic)
# with lmtest::coeftest(), or a mixed model with lme4::lmer(hba1c ~ program + (1|clinic))

import numpy as np, pandas as pd
import statsmodels.formula.api as smf
rng = np.random.default_rng(8)
n_clinic, per = 16, 12
clinic  = np.repeat(np.arange(n_clinic), per)
program = rng.integers(0, 2, n_clinic)
u       = rng.normal(0, 0.8, n_clinic)
hba1c   = 8.0 - 0.4*program[clinic] + u[clinic] + rng.normal(0, 0.7, n_clinic*per)
dia = pd.DataFrame({"hba1c": hba1c, "program": program[clinic], "clinic": clinic})

m  = # ordinary OLS
mc = # OLS with cov_type="cluster"

m  = smf.ols("hba1c ~ program", dia).fit()
mc = smf.ols("hba1c ~ program", dia).fit(
        cov_type="cluster", cov_kwds={"groups": dia["clinic"]})
print(round(m.params["program"], 3),  round(m.bse["program"], 4),  round(m.pvalues["program"], 4))
print(round(mc.params["program"], 3), round(mc.bse["program"], 4), round(mc.pvalues["program"], 4))
# -0.37 0.1394 0.0086   ordinary: looks clearly significant
# -0.37 0.3151 0.24     cluster-robust: standard error doubles, p = 0.24

The coefficient is the same in both rows. The cluster-robust standard error is more than twice the ordinary one, and a result that looked convincing (p = 0.009) becomes weak (p = 0.24). The ordinary model overstated the evidence because it counted 192 near-duplicates as 192 independent facts.

Overfitting and events per variable

The opposite failure is a model that is too large for the data. With many predictors and few observations, the model fits the noise in your particular sample, standard errors balloon, and the fit falls apart on new data. The guard is to count how much outcome information each parameter has to lean on.

For logistic and survival models the count that matters is events, not the total sample. The working rule is at least 10 events per variable:

If 40 of 600 patients died and you want 8 predictors, you have 40 events for 8 parameters, an EPV of 5. That is thin: halve the predictors or accept that the estimates are unstable. For linear regression the rougher version is at least 10 to 15 observations per predictor.

Clustered data: why independence breaks

Clustered data is the most common way the independence assumption fails. It arises in three shapes: repeated measures on the same patient over time, several measurements on one patient (different teeth, both eyes), and patients grouped within a unit (a clinic, a household, a village). In each, observations inside a cluster are more alike than observations across clusters.

The consequence is one-directional and it bites. Correlated observations carry less information than the same number of independent ones, so the effective sample size is smaller than the row count. Ordinary standard errors, computed as if every row were independent, come out too small. Confidence intervals are too narrow and p-values are too small. You saw exactly this in worked example 1. The amount of clustering is summarised by the intra-class correlation, the share of total variance that sits between clusters rather than within them.

Two ways to model clustering

Two model families take clustering seriously, and they answer slightly different questions.

1. Random-effects (multilevel) models. Add a cluster-level term to the linear predictor, assumed to vary randomly between clusters with mean zero. The model splits the variance into between-cluster and within-cluster parts and adjusts both the estimate and its standard error. It answers a subject-specific question: the effect for a given patient or clinic. In R you fit these with lme4::lmer (numerical outcome) or lme4::glmer (binary or count). These need a package beyond base R, so they are described here, not run.
2. Generalized estimating equations (GEE). Specify only the mean and a working correlation structure (commonly exchangeable, meaning every pair in a cluster shares one correlation), then attach sandwich standard errors. GEE answers a population-averaged question: the average effect across the population. It is the natural choice for clustered logistic regression, where random-effects models are harder to fit. In R you would use geepack::geeglm; the runnable version below is in Python.

Two lighter options exist when you do not need individual-level predictors: collapse each cluster to a summary measure and analyse those, or fit the ordinary model and report cluster-robust standard errors as in worked example 1. Match the tool to the question.

Worked example 2: a GEE for clustered binary data

A tuberculosis study traces 179 household contacts of 45 index cases in Selangor. The outcome is a positive tuberculin skin test in the contact. The exposure is whether the index case was smear-positive, a marker of higher infectiousness, recorded at the household level. Contacts in one household share the same index case and home, so their results are clustered. Compare an ordinary logistic regression, which ignores this, with a GEE that does not.

set.seed(3)
n_house <- 45
sizes   <- sample(2:6, n_house, replace = TRUE)
house   <- rep(1:n_house, times = sizes)
smear   <- rbinom(n_house, 1, 0.5)[house]
u       <- rnorm(n_house, 0, 1.4)[house]
p       <- 1 / (1 + exp(-(-0.1 + 0.9 * smear + u)))
pos     <- rbinom(length(p), 1, p)
tb      <- data.frame(pos, smear, house = factor(house))

m <- # glm(pos ~ smear, family = binomial)

m <- glm(pos ~ smear, data = tb, family = binomial)
round(summary(m)$coefficients["smear", ], 4)
exp(coef(m)["smear"])   # odds ratio, ignoring household clustering
# for a GEE in R: geepack::geeglm(pos ~ smear, id = house, family = binomial,
#                                 corstr = "exchangeable")
# for a random-effects logistic model: lme4::glmer(pos ~ smear + (1|house), family = binomial)

import numpy as np, pandas as pd
import statsmodels.api as sm
rng = np.random.default_rng(3)
n_house = 45
sizes = rng.integers(2, 7, n_house)
house = np.repeat(np.arange(n_house), sizes)
smear = rng.integers(0, 2, n_house)[house]
u     = rng.normal(0, 1.4, n_house)[house]
pos   = rng.binomial(1, 1/(1+np.exp(-(-0.1 + 0.9*smear + u))))
tb = pd.DataFrame({"pos": pos, "smear": smear, "house": house})
X  = sm.add_constant(tb["smear"])

glm = # ordinary logistic GLM
gee = # sm.GEE with Exchangeable correlation

glm = sm.GLM(tb["pos"], X, family=sm.families.Binomial()).fit()
gee = sm.GEE(tb["pos"], X, groups=tb["house"],
             family=sm.families.Binomial(),
             cov_struct=sm.cov_struct.Exchangeable()).fit()
print(round(glm.params["smear"], 3), round(glm.bse["smear"], 4), round(np.exp(glm.params["smear"]), 3))
print(round(gee.params["smear"], 3), round(gee.bse["smear"], 4), round(np.exp(gee.params["smear"]), 3))
# 1.481 0.3326 4.399   logistic, clustering ignored: OR 4.40
# 1.360 0.4865 3.896   GEE: OR 3.90, standard error wider

The GEE differs from the ordinary fit on two counts. Its standard error is larger, because it credits the within-household correlation and so trusts the data less. And its coefficient moves a little, because GEE gives less weight to contacts from large households, where each extra contact adds little new information. Ignoring the clustering would have reported a tighter, slightly larger odds ratio than the data support.

Common mistakes

• Treating clustered observations as independent. Repeated measures or grouped patients analysed row by row give standard errors that are too small. Decide the clustering from the design before you fit anything.
• Adjusting an effect modifier away. Putting a modifier in as a plain confounder averages over a real difference in the effect. If you suspect the exposure works differently across a variable, test an interaction rather than only adjusting.
• Cramming predictors with too few events. Nine predictors on 35 events is an EPV of 4. The model overfits and the coefficients swing from sample to sample. Count events per variable before you add a term.
• Reading a likelihood ratio test after asking for cluster-robust standard errors. The log-likelihood does not change when you switch to sandwich standard errors, so a likelihood ratio test still ignores the clustering. Use Wald tests in that setting.
• Transforming when the family is the real fix. Logging a skewed count to force it into a linear model is a workaround. A Poisson or negative-binomial model matches the outcome and keeps the results on a scale you can report.

Tips

• Decide estimation versus prediction first. It dictates which variables belong in the model and how you judge it.
• Read the design for clustering before you open the data. Repeated measures, grouped patients, and multiple sites all break independence.
• Count events per variable for any logistic or survival model. Below 10, trim predictors rather than hope.
• For a population-averaged exposure effect with clustered data, reach for GEE. When the cluster variation is itself of interest, fit a random-effects model.
• If the mean model is sound but the spread is not, sandwich standard errors fix the inference without disturbing the estimate.

The thread through all of this is one idea: match the model to the design and the data. There is no single regression that fits every study. A skewed cost outcome wants a transform or a different family. Patients nested in clinics want a model that knows they are nested. Read the data-generating structure first, then choose the model that respects it.