Introduction
You can get the datasets and SAS commands for this workshop on the web page:
Thelab examples come from the book, “Linear Mixed Models: A Practical Guide Using Statistical Software”, by Brady T. West, Kathleen B. Welch, and Andrzej T. Gałecki, Chapman & Hall/CRC, 2006 (referred to as WWG).The model and figure numbers in the labs correspond to those in WWG.
Go to
to get more detailed analysis for each example in the book, along with the appropriate syntax in SAS (using SAS 9.1.3), SPSS, Stata, R, and HLMsoftware programs for each analysis presented in the book.
In this workshop we will be using four actual datasets plus the hypothetical banana data.
- rat_pup.dat This data set was provided by Jose Pinheiro and Doug Bates in their book, “Mixed-Effects Models in S and S-Plus” (2000). In the analysis of this data set we compare the birth weights of 322 rat pups in 27 litters (litter size varies from 2 to 18 pups) born to mothers treated with High, Low, and Control doses of a chemical.
- classroom.csv This data set was from a study by of Instructional Improvement (SII; Hill, Rowan, and Ball, 2004) and includes data from 1190 first grade students from 312 classrooms in 107 schools.
- ratbrain.datThis data set was reported by Douglas et. al., 2004 and shows measurement of the effects of two different chemicals on three different regions of the brain in five adult male rats.
- autism.csv This data set was derived from a study of 156 children with Autism Spectrum Disorder (Oti, Anderson, Lord, 2006). Measurements were made at five basic ages for each child: 2, 3, 5, 9, and 13 years. Not all children were measured at all time points.
Lab Example 2
Two-Level Clustered Data
The Rat Pup Data
(Chapter 3 in WWG)
This is a two-level clustered data set, in which the clusters are litters, and individual pups are the units of analysis. This analysis is displayed in more detail than the analyses for later examples.
In this data set, 30 pregnant rats were treated with one of three doses of a drug (High, Low, or Control), and the birth weights of therat pupsthat were born were measured.There were originally 10 litters for each of the drug doses. However, 3 litters in the high dose did not survive, resulting in 27 litters. There were 322 rat pups in the final study.
The table below shows a portion of the Rat Pup data. We have Level 1 covariates (e.g. Sex) that are specific for each rat pup and Level 2 covariates (e.g. Treatment and Litsize) that are constant for all rat pups within a given litter. The response variable, Weight, varies for each rat pup.
Portion of the Rat Pup Data Set
LITTER / TREATMENT / LITSIZE / PUP_ID / WEIGHT / SEX1 / Control / 12 / 1 / 6.60 / Male
1 / Control / 12 / 2 / 7.40 / Male
1 / Control / 12 / 3 / 7.15 / Male
1 / Control / 12 / 4 / 7.24 / Male
1 / Control / 12 / 5 / 7.10 / Male
1 / Control / 12 / 6 / 6.04 / Male
1 / Control / 12 / 7 / 6.98 / Male
1 / Control / 12 / 8 / 7.05 / Male
1 / Control / 12 / 9 / 6.95 / Female
…
11 / Low / 16 / 132 / 5.65 / Male
11 / Low / 16 / 133 / 5.78 / Male
…
21 / High / 14 / 258 / 5.09 / Male
21 / High / 14 / 259 / 5.57 / Male
Data Exploration
We first examine the data using descriptive statistics and take a graphical look at the data using Boxplots.
data ratpup;
infile "rat_pup.dat" firstobs=2 dlm="09"X;
input pup_id weight sex $litter litsize treatment $;
if treatment = "High" then treat = 1;
if treatment = "Low" then treat = 2;
if treatment = "Control" then treat = 3;
run;
proc format;
value trtfmt 1 = "High"
2 = "Low"
3 = "Control";
run;
In the data step we set up the new variable, Treat, so that the highest value (=3) is coded for Control, so that Control will be the “reference level” for Treat.
title "Summary statistics for weight by treatment and sex";
proc means data=ratpup maxdec=2;
class treat sex;
var weight;
format treat trtfmt.;
run;
Analysis Variable : weightTreat / Sex / N Obs / N / Mean / Std Dev / Minimum / Maximum
High / Female / 32 / 32 / 5.85 / 0.60 / 4.48 / 7.68
Male / 33 / 33 / 5.92 / 0.69 / 5.01 / 7.70
Low / Female / 65 / 65 / 5.84 / 0.45 / 4.75 / 7.73
Male / 61 / 61 / 6.03 / 0.38 / 5.25 / 7.13
Control / Female / 54 / 54 / 6.12 / 0.69 / 3.68 / 7.57
Male / 77 / 77 / 6.47 / 0.75 / 4.57 / 8.33
There is a tendency for the mean Weight to be lower for the High and Low treatment groups, compared to the Control. We also see that the mean for males is somewhat higher than for females within each level of Treat.
We now look at boxplots of Weight for each combination of treatment and sex.
We now graphically examine the relationship between Litter Size and Weight within each Treatment. We first create a new variable called RANKLIT that ranks the litters by size.
/*Create Ranklit*/
proc sort data=ratpup;
by litsize litter;
run;
data ratpup2;
set ratpup;
by litsize litter;
if first.litter then ranklit+1;
label ranklit="New Litter ID (Ordered by Size)";
run;
There appears to be a relationship between litter size and birthweight, so that larger litters tend to have smaller pups, within each treatment and sex, except for Low Males.
Models
We consider several competing models for the Rat Pup Data. We use a top-down model fitting strategy, where we first set up a “loaded” model that includes all the candidate fixed effects, then we decide on an appropriate covariance structure for the random effects, and finally we reduce the fixed effects portion of the model to arrive at our final model.
The table below shows a summary of each of the models we fit. You can use the model numbers to keep track of what we're doing as we move through the examples.
Term / Variable / General Notation / Model3.1 / 3.2A / 3.2B / 3.3
Fixed Effects / Intercept / 0 / / / /
TREAT1
(High vs. Control) / 1 / / / /
TREAT2
(Low vs. Control) / 2 / / / /
SEX1
(Female vs. Male) / 3 / / / /
LITSIZE / 4 / / / /
TREAT1 × SEX1 / 5 / / /
TREAT2 × SEX1 / 6 / / /
Random Effects / Litter
(j) / Intercept / uj / / / /
Residuals / Rat Pup (pup i in litter j) / εij / / / /
Covariance
Parameters (θD)
for D matrix / Litter Level / Variance of intercepts / σ2litter / / / /
Covariance Parameters (θR)
for Ri matrix / Rat Pup Level / Variances of Residuals / σ2high
σ2low
σ2control / σ2 / σ2high
σ2low
σ2control / σ2high/low,
σ2control / σ2high/low,
σ2control
Model Specification
This model includes the fixed effects of Treatment and Sex and their interaction. It also includes the fixed effect of LITSIZE. There is a random intercept (uj ) for each litter and a random error (εij) for each rat pup.
Hypothesis Tests
A summary of the hypotheses tested in the analysis of the Rat Pup data is shown in the following table. Note that different test types and sometimes different estimation methods are considered for the different hypothesis tests. We use Likelihood Ratio tests (LRT) with REML estimation to help us decide on the appropriate covariance structure; we use an F-test with REML estimation to test the Fixed Effects of the Treatment by Sex interaction; and we use an LRT test with ML estimation to decide whether to keep the Fixed Effects of Treatment in the model.
Hypothesis Specification / Hypothesis TestLabel / Null
(H0) / Alternative
(HA) / Test / Models Compared / Estimation Method / Asymptotic / Approximate Dist. of Test Statistic under H0
Nested Model
(H0) / Reference Model
(HA)
3.1 / Drop u0j
(σ2litter = 0) / Retain u0j
(σ2litter > 0) / LRT / Model
3.1A / Model
3.1 / REML / 0.5χ20 + 0.5χ21
3.2 / Homogeneous Residual Variance
(σ2high = σ2low = σ2control= σ2) / Residual Variances are not all equal / LRT / Model
3.1 / Model 3.2A / REML / χ22
3.3 / Grouped
Heterogeneous
Residual Variance
(σ2high = σ2low) / σ2high ≠ σ2low / LRT / Model
3.2B / Model 3.2A / REML / χ21
3.4 / Homogeneous Residual Variance
(σ2high/low = σ2control= σ2) / GroupedHeterogeneous Residual Variance
(σ2high/low ≠
σ2control) / LRT / Model 3.1 / Model 3.2B / REML / χ21
3.5 / Drop TREATMENT × SEX effects
(5 = 6 = 0) / 5 ≠ 0, or
6 ≠ 0 / Type III
F-test / N/A / Model 3.2B / REML / F(2,194)1
3.6 / Drop TREATMENT effects
(1 = 2 = 0) / 1 ≠ 0, or
2 ≠ 0 / LRT / Model 3.3A / Model 3.3 / ML / χ22
Type III F-test / N/A / REML / F(2,24.3)
Step 1: Fit a model with a “loaded” mean structure (Model 3.1).
The SAS commands used to fit Model 3.1 to the Rat Pup data using Proc Mixedare shown below.
title "Model 3.1";
procmixed data=ratpup order=internal covtest ;
class treat sex litter ;
model weight = treat sex litsize treat*sex /
solution ddfm = sat ;
random intercept / subject=litter ;
format treat trtfmt. ;
run ;
Annotated output from running these commands is shown below:
Model 3.1
The Mixed Procedure
Model Information
Data Set WORK.RATPUP
Dependent Variable weight
Covariance Structure Variance Components
Subject Effect litter
Estimation Method REML
Residual Variance Method Profile
Fixed Effects SE Method Model-Based
Degrees of Freedom Method Satterthwaite
Class Level Information
Class Levels Values
treat 3 High Low Control
sex 2 Female Male
litter 27 1 2 3 4 5 6 7 8 9 10 11 12 13
14 15 16 17 18 19 20 21 22 23
24 25 26 27
Dimensions
Covariance Parameters 2
Columns in X 13
Columns in Z Per Subject 1
Subjects 27
Max Obs Per Subject 18
Number of Observations
Number of Observations Read 322
Number of Observations Used 322
Number of Observations Not Used 0
Iteration History
Iteration Evaluations -2 Res Log Like Criterion
0 1 490.50994069
1 3 401.17987994 0.00076832
2 1 401.10574867 0.00001581
3 1 401.10432307 0.00000001
Convergence criteria met.
Covariance Parameter Estimates
Standard Z
Cov Parm Subject Estimate Error Value Pr Z
Intercept litter 0.09649 0.03279 2.94 0.0016
Residual 0.1635 0.01351 12.10 <.0001
Fit Statistics
-2 Res Log Likelihood 401.1
AIC (smaller is better) 405.1
AICC (smaller is better) 405.1
BIC (smaller is better) 407.7
Solution for Fixed Effects
Standard
Effect sex treat Estimate Error DF t Value Pr > |t|
Intercept 8.3233 0.2733 32.9 30.45 <.0001
treat High -0.9061 0.1915 30.8 -4.73 <.0001
treat Low -0.4670 0.1582 28.2 -2.95 0.0063
treat Control 0 . . . .
sex Female -0.4117 0.07315 295 -5.63 <.0001
sex Male 0 . . . .
litsize -0.1284 0.01875 31.8 -6.85 <.0001
treat*sex Female High 0.1070 0.1318 304 0.81 0.4173
treat*sex Male High 0 . . . .
treat*sex Female Low 0.08387 0.1057 299 0.79 0.4281
treat*sex Male Low 0 . . . .
treat*sex Female Control 0 . . . .
treat*sex Male Control 0 . . . .
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
treat 2 24.3 11.49 0.0003
sex 1 303 46.99 <.0001
litsize 1 31.8 46.87 <.0001
treat*sex 2 302 0.47 0.6282
Step 2: Select a structure for the random effects
We want to decide if we should include the random intercept for each litter in our model, so we need to test Hypothesis 3.1 shown below:
H0: σ2litter = 0 (drop the random litter effects)
HA:σ2litter > 0 (keep the random litter effects)
We now fit a nested model, excluding the random intercept for each litter (the Random statement is excluded from our SAS code).
title "Model 3.1A";
proc mixed data=ratpup order=internal covtest ;
class treat sex litter ;
model weight = treat sex treat*sex litsize /
solution ddfm = sat ;
*random intercept / subject=litter;
format treat trtfmt. ;
run;
The fit statistics that result after fittingModel 3.1A are shown below:
Fit Statistics
-2 Res Log Likelihood 490.5
AIC (smaller is better) 492.5
AICC (smaller is better) 492.5
BIC (smaller is better) 496.3
We can now compare the two models to decide if we wish to keep the random litter effects in the model (i.e. if σ2litter. = 0). We use a likelihood ratio test comparing the -2 Res log likelihood for Model 3.1 to that of Model 3.1A .
Because we are testing a null hypothesis that is on the boundary of the parameter space, we must use a p-value based ona mixture of chi-square distributions to get the appropriate p-value for this test. We use 0.5χ20 + 0.5χ21.
However, we can ignore the χ20portion, and we simply end up with a p-value that isone-half the p-value for χ21 as shown below.
data _null_;
lrtstat = 490.5 - 401.1;
df = 1;
pvalue = 0.5*(1 - probchi(lrtstat,df));
format pvalue 10.8;
run;
title "P-value for variance of random litter effects";
proc print data=test1;
run;
The results of this test are displayed in the output window:
Obs lrtstat df pvalue
1 89.4 1 0.00000000
Because this result is significant, we will keep the random intercept for each litter in this model and for all subsequent models.
We could stop our analysis at this point and be satisfied with the results, or we could decide to drop the treatment by sex interaction term from the model. However, we will first think about whether we really have equal variances for all of the residuals for all of the treatments.
Step 3: Select a covariance structure for the residuals
One of the very nice features of SAS Proc Mixed is that it allows you to have unequal residual variances for different groups of subjects. We will investigate this option in this part of the lab exercise.
Based on our examination of the boxplots of Weight for different litters by levels of Treatment and sex, we noted that the variance in the Control treated litters tended to be larger than that observed in the High and Low treated litters..
We can build this heterogeneity of variance into the SAS code, using the Repeated statement (which is used to specify the residual covariance structure), as shown in the syntax for Model 3.2A below. Notice that we are only using the repeated statement because we wish to modify the residual variance from being constant for all treatments, to being different for different treatments.
title "Model 3.2A";
proc mixed data=ratpup order=internal;
class treat litter sex;
model weight = treat sex treat*sex litsize /
solution ddfm=sat ;
random int / subject=litter;
repeated / group=treat;
format treat trtfmt.;
run;
The residual variance estimates for each of the treatments based on this model is shown below. Notice that the residual variance for Control is about twice as high as for the High and Low treatments, and that the variance for the High and Low treatments is very similar (.1084 and .0846, respectively).
Covariance Parameter Estimates
Cov Parm Subject Group Estimate
Intercept litter 0.09827
Residual treat High 0.1084
Residual treat Low 0.08459
Residual treat Control 0.2650
The fit statistics for this model are shown below:
Fit Statistics for Model 3.2A
-2 Res Log Likelihood 359.9
AIC (smaller is better) 367.9
AICC (smaller is better) 368.0
BIC (smaller is better) 373.1
We can now proceed to test the following hypothesis:
H0: σ2high = σ2low = σ2control = σ2 (the residual variances are all equal)
HA: At least one pair of residual variances is not equal
Because we are no longer on the boundary of a parameter space, the REML-based likelihood ratio test that we calculate has 2 df, which correspond to the two additional covariance parameters in Model 3.2A vs. Model 3.1.
title "P-value for Hypothesis 3.2";
data _null_;
lrtstat = 401.1-359.9;
df=2;
pvalue = 1 - probchi(lrtstat,df);
format pvalue 10.8;
put lrtstat=df=pvalue=;
run;
The results of this hypothesis test are displayed in the SAS log . The result is significant, so we decide that Model 3.2A is preferable to Model 3.1:
lrtstat=41.2 df=2 pvalue=0.000000
We now further refine the residual covariance structure by deciding if we need to have separate residuals variances for each level of treatment, or if we can pool the residual variance for the High and Low treatment groups. To do this, we fit Model 3.2B, but first we create a new variable, called TRTGRP, combining the High and Low treatment groups.
title "RATPUP3 dataset";
data ratpup3;
set ratpup2;
if treatment in ("High", "Low") then TRTGRP = 1;
if treatment = "Control" then TRTGRP = 2;
run;
proc format;
value tgrpfmt 1="High/Low"
2="Control";
run;
title "Model 3.2B";
proc mixed data=ratpup3 order=internal;
class treat litter sex trtgrp;
model weight = treat sex treat*sex litsize /
solution ddfm=sat;
random int / subject=litter;
repeated / group=trtgrp;
format treat trtfmt. trtgrp tgrpfmt.;
run;
The covariance parameter estimates for Model 3.2B are shown below:
Covariance Parameter Estimates
Cov Parm Subject Group Estimate
Intercept litter 0.09895
Residual TRTGRP High/Low 0.09242
Residual TRTGRP Control 0.2650
The fit statistics for this model are:
Fit Statistics for Model 3.2B
-2 Res Log Likelihood 361.1
AIC (smaller is better) 367.1
AICC (smaller is better) 367.2
BIC (smaller is better) 371.0
We can now test the hypotheses shown below:
H0: σ2high = σ2low
HA: σ2high ≠ σ2low
We use a REML-based likelihood ratio test. The SAS code for this hypothesis test is shown below:
title "P-value for Hypothesis 3.3";
data _null_;
lrtstat = 361.1-359.9;
df = 1;
pvalue = 1-probchi(lrtstat, df);
format pvalue 10.8;
put lrtstat = df = pvalue =;
run;
The results of this test are not significant (p=0.27), so we decide to use the more simple model, Model 3.2B as the preferred model at this stage, and we pool the residual variance for the High and Low treatment groups.
lrtstat=1.2 df=1 pvalue=0.27332168
Hypothesis 3.4 involves comparing Model 3.2B to Model 3.1.
H0: σ2high/low = σ2control= σ2
HA: σ2high/low ≠ σ2control
We do not show syntax for Hypothesis 3.4 here. Instead, we simply work with Model 3.2Bwith unequal variances for the treated and control litters as our preferred model.
Step 4: Reduce the model by removing nonsignificant fixed effects.
We now test Hypothesis 3.5 to decide whether we can remove the Fixed Effects associated with the Treatment by Sex interaction. For this test we use an approximate F-test, with the Satterthwaite degrees of freedom, from our current model, Model 3.2B.
The Type 3 F-test results for Fixed Effects for Model 3.2B are shown below. We decide to drop the Treatment by Sex interaction (p=.7289) from the model.
Model 3.2B
Type 3 Tests of Fixed Effects
Num Den
Effect DF DF F Value Pr > F
treat 2 24.4 11.18 0.0004
sex 1 296 59.17 <.0001
treat*sex 2 194 0.32 0.7289
litsize 1 31.2 49.33 <.0001
We can now test Hypothesis 3.6. We choose to do this using a likelihood ratio test. Because the two models that we wish to compare differ only in the Fixed Effects that they include, we use Maximum Likelihood (method=ml) to fit both models, and then compare them.
title "Model 3.3 using ML";
proc mixed data=ratpup3 order=internal method=ml;
class treat litter sex trtgrp;
model weight = treat sex litsize /solution ddfm=sat ;
random int / subject=litter;
repeated / group=trtgrp;
format treat trtfmt.;
run ;
title "Model 3.3A using ML";
proc mixed data=ratpup3 order=internal method=ml;
class treat litter sex trtgrp;
model weight = sex litsize /solution ddfm=sat ;
random int / subject=litter;
repeated / group=trtgrp;
format treat trtfmt.;
run ;
The fit statistics for Model 3.3 and Model 3.3A, both fit using ML estimation are shown below: