Fitting GLMMs with SAS PROC GLIMMIX

From a clustered binary outcome to SAS statements

For a binary response measured on multiple individuals within each cluster, the model must represent both the event probability and the dependence among observations in the same cluster. PROC GLIMMIX expresses these choices through separate statements: the response and fixed effects go in MODEL, categorical variables are identified in CLASS, and the random effects go in RANDOM. The distribution and link are options on MODEL.

The workflow below uses a hypothetical study in which patients are nested within clinics, the outcome is whether a patient improved, and the main predictor is treatment group. It is a random-intercept logistic GLMM. The goal here is to translate that model into a clear SAS specification and check that the fitted model corresponds to the intended analysis.

Model definition: Let $Y_{ij}=1$ if patient $i$ in clinic $j$ improves and $Y_{ij}=0$ otherwise. Conditional on a clinic random intercept $b_j$, assume $Y_{ij}$ has a Bernoulli distribution with probability $p_{ij}$, where $$ \operatorname{logit}(p_{ij})=\log\left(\frac{p_{ij}}{1-p_{ij}}\right)=\beta_0+\beta_1\text{Treatment}_{ij}+\beta_2\text{AgeC}_{ij}+b_j,\qquad b_j\sim N(0,\sigma_b^2). $$ Here Treatment is 1 for the active treatment and 0 for control; AgeC is age centered at 50 years. The model assumes independent random effects across clinics and conditionally independent outcomes given the random effects and covariates. It also assumes the specified linear relationship on the log-odds scale and a normal random-intercept distribution.

The coefficient $\beta_1$ describes the treatment difference in log odds for patients of the same age and clinic random effect. Thus $\exp(\beta_1)$ is a conditional odds ratio, not a population-averaged odds ratio. The random intercept allows the baseline log odds to vary across clinics.

A complete PROC GLIMMIX specification

In SAS, the CLASS statement declares variables that represent categories. The MODEL statement identifies the response and fixed effects; DIST=BINARY and LINK=LOGIT specify the Bernoulli-logit model. The RANDOM statement identifies a random intercept and the variable defining which observations share it. The example requests a Laplace approximation for likelihood fitting and fixed-effect solutions in the output.

Worked Example: Fitting a Clinic Random-Intercept Model

Suppose a hypothetical teaching data set has one record per patient, with clinic, treatment group, age in years, and improvement coded 1 for yes and 0 for no. A small excerpt is shown below; the complete analysis would use all eligible patients.

Clinic	Treatment	Age (years)	Improved
01	Control	48	0
01	Active	55	1
02	Control	51	0
02	Active	46	1

The response is Improved, the fixed effects are Treatment and centered age, and the random effect is the clinic-specific intercept. Centering age gives AgeC = age minus 50, so the intercept corresponds to the control group at age 50 when the random intercept is zero.

One way to prepare the variables and fit the model is:

SAS workflow:
data patients;
  input clinic $ treatment $ age improved;
  agec = age - 50;
  datalines;
01 Control 48 0
01 Active 55 1
02 Control 51 0
02 Active 46 1
;
run;

proc glimmix data=patients method=laplace;
  class clinic treatment(ref='Control') / param=ref;
  model improved(event='1') = treatment agec /
    dist=binary link=logit solution;
  random intercept / subject=clinic;
run;

The excerpt is only to illustrate variable roles and syntax; four records are not adequate for reliable estimation of a clinic-level variance. In a real analysis, the data would have enough clinics and observations to support the model. Before fitting, check that each row represents one patient, that clinic identifiers are correct, and that Improved contains only the intended 0 and 1 values. With PARAM=REF and Control as the reference, the treatment coefficient compares Active with Control.

Event coding and the direction of the comparison

The event option makes the modeled outcome level explicit. For a 0/1 response, EVENT='1' makes the modeled probability the probability of improvement. The response coding and option must agree: reversing the event changes the direction of the log odds and reverses the interpretation of the treatment contrast. For a categorical response, also check the formatted levels and SAS's ordering rules rather than assuming that a default level is the desired event.

Worked Example: Checking Which Outcome Is Modeled

Imagine a hypothetical data check finds 68 records coded Improved=1 and 92 coded Improved=0. The scientific question concerns improvement, so the event must be 1. With EVENT='1', the model targets $P(\text{Improved}=1\mid\text{predictors},b_j)$. If the event were instead 0, the modeled probability would be non-improvement.

The consequences are visible algebraically. Since $P(Y=0)=1-P(Y=1)$, $$ \operatorname{logit}\{P(Y=0)\}=-\operatorname{logit}\{P(Y=1)\}. $$ Thus, with the same predictors and parameterization, coefficients for the complementary event have the opposite sign. For example, a treatment log-odds estimate of 0.60 for improvement corresponds to a log-odds estimate of -0.60 for non-improvement, and the corresponding odds ratios are $\exp(0.60)=1.82$ and $\exp(-0.60)=0.549$, respectively. The data check establishes that both levels occur; it does not by itself establish that the chosen coding is scientifically correct.

Reading fixed effects on the model scale

PROC GLIMMIX reports fixed-effect estimates on the scale of the linear predictor. For the logit link, these estimates are log odds. Exponentiating a coefficient gives a conditional odds ratio, holding the other fixed predictors and the cluster random effect constant. A confidence interval for a coefficient can likewise be exponentiated endpoint by endpoint to obtain a confidence interval for its odds ratio.

Worked Example: Converting a Hypothetical Treatment Estimate

Suppose a hypothetical PROC GLIMMIX output reports a treatment estimate of 0.60 with standard error 0.25. For a large-sample Wald calculation, the test statistic for $H_0:\beta_1=0$ is $$ z=\frac{0.60}{0.25}=2.40. $$ A two-sided standard normal p-value is approximately 0.0164. At a 0.05 significance level, this is evidence against a zero conditional treatment coefficient, under the fitted model and its inferential approximation.

The estimated conditional odds ratio is $$ \exp(0.60)=1.82. $$ A 95% Wald interval on the log-odds scale is $0.60\pm1.96(0.25)$, or $0.11$ to $1.09$. Exponentiating gives $\exp(0.11)=1.12$ and $\exp(1.09)=2.97$. Under this model, the estimated odds of improvement for active treatment are 1.82 times the odds for control, conditional on age and the clinic random intercept; the approximate 95% interval is 1.12 to 2.97. This is not a statement that the probability of improvement is 1.82 times as large.

Justification for the conditional interpretation

The interpretation follows directly from the model. For two patients in the same clinic and of the same age, let the first receive active treatment and the second control. Their linear predictors differ by $\beta_1$, because age and the shared $b_j$ cancel. Since odds equal the exponential of log odds, their conditional odds ratio is

\frac{p_{\text{Active}}/(1-p_{\text{Active}})} {p_{\text{Control}}/(1-p_{\text{Control}})} =\frac{\exp(\beta_0+\beta_1+\beta_2\text{AgeC}+b_j)} {\exp(\beta_0+\beta_2\text{AgeC}+b_j)} =\exp(\beta_1).

This cancellation relies on the model having the same treatment coefficient across clinics and no treatment-by-clinic random slope. If treatment effects vary by clinic, the random-intercept model does not describe that extra variation.

Checks, alternatives, and limitations

A successful run is not sufficient evidence that the intended model was fitted. Check that the response has the intended event level, that reference categories and contrast coding are as planned, and that the number of observations and clusters match the analysis data. Review convergence messages, the fixed-effect estimates and standard errors, and the estimated random-intercept variance. A variance estimate near zero can indicate little detectable between-clinic variation, but can also occur when the data provide limited information about that variance.

Warning signs include complete or near separation, very sparse outcome levels within many clinics, extremely large estimates or standard errors, a non-positive-definite Hessian, and failure to converge. These can indicate weak information, over-complex fixed effects, coding problems, or an inadequate model. Recheck the data and model specification before relying on Wald tests. With few clusters, ordinary large-sample approximations can be unreliable; report the fitting method and treat inferential results cautiously.

The Laplace approximation is a likelihood-based approach to integrating over random effects. It can be computationally more demanding than pseudo-likelihood methods, and estimates may differ when data are sparse or cluster sizes are small. Check whether reasonable alternative fitting settings materially change the conclusions. The earlier material on GLMM software contrasts this model-based approach with alternatives: a logistic model that omits the random effect ignores within-clinic dependence, while a marginal model targets a different, population-averaged effect. Neither is automatically interchangeable with the conditional GLMM.

A random intercept assumes a common distribution of clinic baselines and a common treatment effect across clinics. It does not correct confounding, selection bias, or incorrect measurement, and it does not make an observational treatment comparison causal. Conclusions also depend on the chosen link, covariate specification, and random-effects assumptions. State these limits when reporting results.

Common pitfall

Common mistake: “The treatment coefficient is 0.60, so treatment increases the probability of improvement by 0.60.” This is incorrect because the coefficient is on the conditional log-odds scale, not the probability scale. Exponentiating it gives a conditional odds ratio of 1.82; probability differences depend on the starting probability and the random effect. Also verify the modeled event before assigning a direction to the result.

Check Your Understanding

Answer each question, then compare with the solutions.

In a model with patients nested in clinics, which statement identifies the clinic random intercept in PROC GLIMMIX?
If improvement is coded 1 and non-improvement 0, what does EVENT='1' select?
A hypothetical treatment estimate is -0.30. What is its conditional odds ratio, rounded to three significant figures?
Name two checks to make if PROC GLIMMIX reports a convergence warning.
Why is a random-intercept model not automatically equivalent to a population-averaged logistic model?

Solutions. (1) random intercept / subject=clinic; specifies a separate random intercept for each clinic. (2) It selects the probability of improvement as the modeled event. (3) $\exp(-0.30)=0.741$, so the estimated conditional odds are about 0.741 times those for the reference treatment, at the same covariate values and random intercept. (4) Check response and predictor coding, reference levels, sparse or separated outcome patterns, the data and cluster counts, and the convergence or Hessian diagnostics. (5) The random-intercept model is conditional on the random effect; averaging probabilities over its distribution produces a generally different marginal effect.