### SAS Textbook Examples Multilevel Analysis Techniques and Applications by Joop Hox Chapter 2: The Basic Two-Level Regression Model: Introduction

The data set used in this chapter is popular.

Table 2.1 on page 17.

Part 1: Intercept only.

data pop;
set popular;
run;
proc mixed data = pop ;
model popular = /solution;
random intercept / subject = school type = un;
run;
The Mixed Procedure
Model Information
Data Set                     WORK.POP
Dependent Variable           POPULAR
Covariance Structure         Unstructured
Subject Effect               SCHOOL
Estimation Method            REML
Residual Variance Method     Profile
Fixed Effects SE Method      Model-Based
Degrees of Freedom Method    Containment

Dimensions
Covariance Parameters             2
Columns in X                      1
Columns in Z Per Subject          1
Subjects                        100
Max Obs Per Subject              26
Observations Used              2000
Observations Not Used             0
Total Observations             2000

Covariance Parameter Estimates
Cov Parm     Subject    Estimate
UN(1,1)      SCHOOL       0.8798
Residual                  0.6387

Fit Statistics
-2 Res Log Likelihood          5115.6
AIC (smaller is better)        5119.6
AICC (smaller is better)       5119.6
BIC (smaller is better)        5124.8
Null Model Likelihood Ratio Test
DF    Chi-Square      Pr > ChiSq
1       1379.30          <.0001

Solution for Fixed Effects
Standard
Effect       Estimate       Error      DF    t Value    Pr > |t|
Intercept      5.3076     0.09550      99      55.58      <.0001

Part 2: The variable sex is included as a random effect and teacher experience (texp) as fixed effect (M1).

proc mixed data = pop;
model popular = sex texp /solution;
random intercept sex /subject = school type=un ;
run;
The Mixed Procedure

Covariance Parameter Estimates
Cov Parm     Subject    Estimate
UN(1,1)      SCHOOL       0.4116
UN(2,1)      SCHOOL      0.02089
UN(2,2)      SCHOOL       0.2733
Residual                  0.3925

Fit Statistics
-2 Res Log Likelihood          4275.9
AIC (smaller is better)        4283.9

Fit Statistics
AICC (smaller is better)       4283.9
BIC (smaller is better)        4294.3

Null Model Likelihood Ratio Test
DF    Chi-Square      Pr > ChiSq
3       1279.56          <.0001

Solution for Fixed Effects
Standard
Effect       Estimate       Error      DF    t Value    Pr > |t|
Intercept      3.3400      0.1608      98      20.77      <.0001
SEX            0.8431     0.05969      99      14.13      <.0001
TEXP           0.1084     0.01022    1800      10.61      <.0001

Type 3 Tests of Fixed Effects
Num     Den
Effect         DF      DF    F Value    Pr > F
SEX             1      99     199.54    <.0001
TEXP            1    1800     112.51    <.0001
Table 2.2. on page 20.
Part 1: M1, as shown above.
Part 2: Cross-level interaction with teacher experience (M2). We first created a variable genxexp an interaction term of sex and terp. The variable genxexp is included as a fixed effect.
data pop1;
set pop;
genxexp = sex*texp;
run;
proc mixed data = pop1 ;
model popular = sex texp genxexp/solution;
random intercept sex /subject = school type=un ;
run;
The Mixed Procedure
Covariance Parameter Estimates
Cov Parm     Subject    Estimate
UN(1,1)      SCHOOL       0.4120
UN(2,1)      SCHOOL      0.02343
UN(2,2)      SCHOOL       0.2264
Residual                  0.3924

Fit Statistics
-2 Res Log Likelihood          4268.4
AIC (smaller is better)        4276.4
AICC (smaller is better)       4276.5
BIC (smaller is better)        4286.9

Null Model Likelihood Ratio Test
DF    Chi-Square      Pr > ChiSq
3       1274.81          <.0001

Solution for Fixed Effects
Standard
Effect       Estimate       Error      DF    t Value    Pr > |t|
Intercept      3.3135      0.1610      98      20.58      <.0001
SEX            1.3296      0.1330      98       9.99      <.0001
TEXP           0.1102     0.01023    1800      10.77      <.0001
genxexp      -0.03403    0.008457    1800      -4.02      <.0001

Type 3 Tests of Fixed Effects
Num     Den
Effect         DF      DF    F Value    Pr > F
SEX             1      98      99.86    <.0001
TEXP            1    1800     116.07    <.0001
genxexp         1    1800      16.20    <.0001

Table 2.3 on page 21.

Part 1: M1 from Table 2.2.
Part 2: Standardized variables. We first standardized all the variables

proc stdize data = pop out=pops;
var popular sex texp;
run;
proc mixed data = pops;
model popular = sex texp / solution ;
random intercept sex /subject = school type=un ;
run;
Covariance Parameter Estimates
Cov Parm     Subject    Estimate
UN(1,1)      SCHOOL       0.3305
UN(2,1)      SCHOOL      0.05122
UN(2,2)      SCHOOL      0.04545
Residual                  0.2612

Fit Statistics
-2 Res Log Likelihood          3460.0
AIC (smaller is better)        3468.0

Fit Statistics
AICC (smaller is better)       3468.0
BIC (smaller is better)        3478.4

Null Model Likelihood Ratio Test
DF    Chi-Square      Pr > ChiSq
3       1279.56          <.0001

Solution for Fixed Effects
Standard
Effect       Estimate       Error      DF    t Value    Pr > |t|
Intercept    -0.00978     0.05867      98      -0.17      0.8680
SEX            0.3439     0.02434      99      14.13      <.0001
TEXP           0.5791     0.05459    1800      10.61      <.0001

Type 3 Tests of Fixed Effects
Num     Den
Effect         DF      DF    F Value    Pr > F
SEX             1      99     199.54    <.0001
TEXP            1    1800     112.51    <.0001

Figure 2.1 on page 23. Return back to model M1. The outp option requests the predicted values, residuals and other statistics to be saved in a temporary file called test.

proc mixed data = pop;
model popular = sex texp /solution outp=test;
random intercept sex / subject = school type = un;
run;
proc stdize data= test out=fig2;
var resid;
run;
proc univariate data = fig2 noprint;
var resid;
qqplot /href = 0 vref = 0 ;
run;

Figure 2.2 Level 1 standardized residuals plotted against predicted popularity using the fixed part of the model.

data pred;
set test;
eqpred =  3.3400 + 0.8431* sex + .1084*texp;
run;

proc stdize data = pred out = pred1;
var resid;
run;

axis1 order = (-3.7, -2.8, -1.9, -.9, 0, .9, 1.9, 2.8, 3.7) label = (a=90 "std(const)");
axis2 order = (3.5 to 7.1 by .9) label = ("pred.val.") minor = none;
symbol value = star color = blue;
proc gplot data = pred1;
plot resid*eqpred /vref=0 vaxs = axis1 haxis = axis2;
run;
quit;



Figure 2.3 Level 2 residuals plotted against predicted popularity. We need to use the ods output statement to output the level 2 residuals. The output data file pop2 contains both the residuals for the intercept and variable sex. The data step after proc mixed separates them apart. Then we need to get the average prediction at level 2 and this is done using proc sql. Then we are ready to plot them one at a time. The final step is to put these two plots together and this is done using proc greplay.

goptions nodisplay;
proc greplay igout=mycat tc=tempcat nofs;
delete _all_;
run;
proc mixed data = pop;
model popular = sex texp /solution ;
random intercept sex /solution subject = school type = un;
ods output SolutionR = pop2;
run;

data int sex;
set pop2;
school = subject + 0;
keep school estimate stderrpred;
if effect = "Intercept" then output int;
if effect = "SEX" then output sex;
run;
data pintall;
merge pred int;
by school;
run;
proc sql;
create table pall1 as
select distinct mean(eqpred) as p, estimate as int
from pintall
group by school;
quit;
symbol v = star w = 3 c=blue;
axis1 order = (-1.6 to 1.6 by .4) minor =none label=(a=90 "(const)");
axis2 order = (3.8 to 7 by .8) minor = none label=("pred. val." );
proc gplot data = pall1 gout=mycat ;
format int p 4.1 ;
plot int*p = 1 /vref=0 vaxis=axis1 haxis=axis2;
run;
quit;
data psex;
merge pred sex;
by school;
run;
proc sql;
create table psex2 as
select distinct mean(eqpred) as p, estimate as sex
from psex
group by school;
quit;
symbol v = star w = 3 c=blue;
axis1 order = (-1.2 to 1.2 by .3) label=(a = 90 "(sex)");
axis2 order = (3.8 to 7 by .8) minor = none label=("pred. val.");

proc gplot data = psex2 gout=mycat ;
format sex p 4.1 ;
plot sex*p = 1 /vref=0 vaxis=axis1 haxis=axis2;
run;
quit;
goptions display;
proc greplay nofs igout=mycat tc=sashelp.templt template=v2;
treplay 1:gplot 2:gplot1 ;
run;
quit;



Figure 2.4 Error bar plot of level 2 residuals.

goptions nodisplay;
proc greplay igout=mycat2 tc=tempcat nofs;
delete _all_;
run;

proc sort data =int;
by estimate;
run;

data popint1;
set int;
id = _n_;
est= estimate; id = id; output;
est = estimate + 1.39*stderrpred; id = id; output;
est = estimate - 1.39*stderrpred; id = id; output;
run;

symbol i=hiloT;
axis1 order =(-1.6, -1, -.5, 0, .5, 1, 1.6, 2.1) label=(a=90 "const") minor=none;
axis2 order = (1 to 101 by 25) label = ("Rank") minor=none;
proc gplot data = popint1 gout=mycat2;
format est 4.1;
plot ( est )*id  /overlay vref=0 vaxis=axis1 haxis=axis2;
run;
quit;
proc sort data = sex;
by estimate;
run;

data popsex1;
set sex;
id = _n_;
est = estimate; id = id; output;
est = estimate + 1.39*stderrpred; id = id; output;
est = estimate - 1.39*stderrpred; id = id; output;
run;

symbol i=hiloT;
axis1 order =(-1.6, -1, -.5, 0, .5, 1, 1.6, 2.1) label=(a=90 "const") minor=none;
axis2 order = (1 to 101 by 25) label = ("Rank") minor=none;
proc gplot data = popsex1 gout = mycat2;
format est 4.1;
plot ( est )*id  /overlay vref=0 vaxis=axis1 haxis=axis2;
run;
quit;
goptions display;
proc greplay nofs igout=mycat2 tc=sashelp.templt template=v2;
treplay 1:gplot 2:gplot1 ;
run;
quit;

Figure 2.6. The predicted values from outp = test statement is from the random effect model. We can take away the fixed effect of variable texp.

data test1;
set test;
p = pred - .1084*texp;
run;
proc sort data = test1;
by school;
run;

axis1 order =(1.8 to 6.6 by 1.2 ) minor=none label=(a=90 "Predicted Popularity");
axis2 order = (0 to 1 by 1) minor = none label =("Pupil Sex");
symbol i = join r=100  c = black;
proc gplot data = test1;
plot p*sex=school /overlay nolegend haxis=axis2 vaxis=axis1;
run;
quit; 

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.