BIO226 Lab Session 8: Generalized Linear Mixed Effects Models (GLMMs)

BIO226 Lab Session 8: Generalized Linear Mixed Effects Models(GLMMs) Professor Brent Coull TA: Shira Mitchell May 3, 2012

Key Points of GLMM 1. GLMMs extend the approach of linear mixed effects models to categorical data. 2. GLMMs assume heterogeneity across individuals in a subset of regression coefficients (e.g. intercepts and slopes). 3. While Marginal Models (GEEs) focus on inferences about populations, GLMMs focus on inferences about individuals. 4. Regression parameters from GLMMs have ‘subject specific’ interpretations in terms of changes in the transformed mean response for a specific individual.

Specification of GLMM The GLMM can be considered in 2 steps: 1. Assume conditional distribution of each Yij, given individual-specific effects bi, belongs to exponential family with conditional mean g(E[Yij |bi]) = X′ij β + Z′ijbi, where g(.) is known link function and Zij is known design vector (a subset of Xij) linking random effects bi to Yij.

Specification of GLMM 2. The bi are assumed to vary independently from one individual to another and bi ∼ N(0,G), where G is covariance matrix for random effects. Note: additional assumption of “conditional independence”, i.e. given bi, the responses Yi1, Yi2, ..., Yini are assumed to be mutually independent.

GLMM ExampleLongitudinal Binary Response to Depression Medication Cross-classification of responses on depression at 3 times (N=Normal, A=Abnormal) Reprinted by Agresti (2002) with permission from original source (Koch et al., 1977, Biometrics) The data stored in ‘depress.txt’ have already been converted into long form and contains the following 6 variables: ID, Y (0=Abnormal , 1=Normal ) Severe (0=mild, 1=severe), Drug (0=standard, 1=new), Time, and Drug*Time.

SAS Code DATA depress; INFILE ‘depress.txt’; INPUT id y severe drug time dt; RUN; DATA depress; SET depress; t=time; \* create categorical time variable*\ RUN; PROC PRINT DATA=depress(WHERE=(id=65 or id=101)); RUN;

SAS Output 1. drug*time Question of interest: Do patient-specific changes in probability of normal differ between the two treatments?

Marginal Model for Depression Data To obtain initial parameter values for the GLMM, we fit the following Marginal Model (GEE): logit{Pr(Yij = 1)} = ηij = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej where: • Yij = 0 subject i is abnormal in period j; 1 subject i is normal in period j • severei= 0 mild depression, initial diagnosis; 1 severe depression, initial diagnosis • drugi= 0 standard; 1 new drug • timej= 0 if baseline; 1 if time 1; 2 if time 2 and we assume: • Yij ∼ Bernoulli (eηij/(1+eηij)) • Var(Yij) = E(Yij)(1 − E(Yij)), note that Pr(Yij = 1) = E(Yij) because Yij is binary. • log OR(Yij,Yik) = αjk

SAS Code PROC GENMOD DESCENDING DATA=depress; CLASS id t; MODEL y=severe drug time dt / DIST=binomial LINK=logit; REPEATED SUBJECT=id / WITHINSUBJECT=t LOGOR=fullclust; RUN;

SAS Output • Log Odds Ratio • Parameter Information • Parameter Group • Alpha1 (1, 2) • Alpha2 (1, 3) • Alpha3 (2, 3)

GLMM for Depression Data • Consider the following GLMM: logit{Pr(Yij = 1| bi1)} = ηij = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej+ bi1 where bi1 is a random intercept that allows a different baseline probability of normal (vs abnormal) for each subject. and we assume: • Yij|bi1 ∼ Bernoulli (e ηij/(1+eηij)) which implies that Var(Yij|bi1) = E(Yij|bi1)(1 − E(Yij| bi1)). Note: E(Yij|bi1) = Pr(Yij = 1| bi1) because Yij is binary. • Given bi1, the responses Yi0, Yi1, Yi2, are mutually independent. • The bi1 are assumed to vary independently from one individual to another and bi1 ∼ N(0, σ2b).

NLMIXED in SAS PROC NLMIXED DATA=depress QPOINTS=20; PARAMS beta1=-0.03 beta2=-1.3 beta3=-0.05 beta4=0.48 beta5=1.01 sigma=0.07; eta = beta1 + beta2*severe + beta3*drug + beta4*time + beta5*dt + b1; p = (exp(eta)/(1 + exp(eta)); MODEL y ~ BINARY(p); RANDOM b1 ~ NORMAL(0, sigma*sigma) SUBJECT = id; ESTIMATE ’treatment effect, time 1’ beta3 + beta5; ESTIMATE ’treatment effect, time 2’ beta3 + 2*beta5; ESTIMATE ’time trend standard treatment’ beta4; ESTIMATE ’time trend new treatment’ beta4 + beta5; RUN;

NLMIXED in SAS • PARAMS statement: lists all parameters (fixed effects and covariance for random effects) and their initial values (default initial value is 1). • Program statements: defines linear predictor eta (includes fixed and random effects) and relates mean response (p) to linear predictor (eta). • MODEL statement: specifies response variable and conditional distribution of response given random effects (e.g. BINARY). • RANDOM effects ∼ distribution SUBJECT=variable: defines random effects (RANDOM) and variable that determines clustering of observations within an individual (SUBJECT). Note: PROC NLMIXED does not have a CLASS statement, therefore, it is critical that the dataset is sorted by ID prior to analysis.

Estimate Statements Treatment effect, time 1 logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej+ bi1 For drugi = 0 and timej = 1, logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β4 + bi1 For drug=1 and time=1, logit{Pr(Yi’j = 1| bi’1)} = β1 + β2severei’ + β3 + β4 + β5 + bi’1 Thus, difference = β3 + β5 assuming bi1 = bi’1 and severei = severei’

Estimate Statements Treatment effect, time 2 logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej+ bi1 For drug=0 and time=2, logit{Pr(Yij = 1| bi1)} = β1 + β2severei + 2β4 + bi1. For drug=1 and time=2, logit{Pr(Yi’j = 1| bi’1)} = β1 + β2severei’ + β3 + 2β4 + 2β5 + bi1. Thus, the difference = β3 + 2β5 assuming bi1 = bi’1 and severei = severei’

Estimate Statements logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej+ bi1 Time Trend, Standard Treatment logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β4timeij + bi1. Time Trend, New Treatment logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β3 + (β4 +β5)timeij + bi1.

NL Mixed Output

Not significant at baseline (RCT) NL Mixed Output logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej+bi1

NL Mixed Output logit{Pr(Yij = 1| bi1)} = β1 + β2severei + β3drugi + β4timej + β5drugi ∗ timej+bi1

Conclusions • Research question: are patient-specific changes in probability of normal different between the two treatments over time? This corresponds to a testing H0: β5 = 0 • β5 = 1.0184 (p-value<.0001). Thus, we reject H0 of no treatment effect and conclude that there are greater patient-specific changes in probability of normal for the new treatment. • The estimated odds ratio of normal comparing a patient on the new treatment to a patient on the standard treatment with the same random intercept and severity of initial diagnosis is 2.61 (1.93, 3.52) [e0.9587(e.659,e1.258)] for time 1, and 7.22 (4.28, 12.19) [e1.977(e1.453, e2.501)] for time 2.

Conclusions, continued • We estimate that the odds of normal for a subject on standard treatment increases by a factor of 1.62 (e0.483) for each time period. We estimate that the odds of normal for a subject on the new treatment increases by a factor of 4.49 (e1.501) for each time period. • The odds of normal of a subject with an initial diagnosis of severe depression are 0.27 (e−1.315) times the odds of normal of a subject with mild depression and the same random intercept (i.e., a lower odds of normal). • There appears to be little heterogeneity among subjects (σb = 0.06583). Approximately 95% of patients in the standard group with an initial diagnosis of mild depression are expected to have a baseline (time=0) log odds of normal between -0.1568 and 0.1011 (−0.02795 ± 1.96 × 0.06583) or baseline probability of normal between 0.461 = e-0.1568/(1+e-0.1568) and 0.525 = e0.1011/(1+e0.1011). (Lecture 20 slide 23)

Conclusions, continued Note that, when we interpret the parameter estimates from the mixed model, we interpret them at the patient level. When we report odds ratios comparing two patients, we assume that they have the same random intercepts (i.e. the same baseline propensity for normal).

BIO226 Lab Session 8: Generalized Linear Mixed Effects Models (GLMMs)