Example 41.3: Plotting the Likelihood
The data for this example are from Hemmerle and Hartley (1973) and
are also used as an example for the VARCOMP procedure. The response
variable consists of measurements from an oven experiment, and the
model contains a fixed effect A and random effects
B and A*B.
The SAS code is as follows:
data hh;
input a b y @@;
datalines;
1 1 237 1 1 254 1 1 246
1 2 178 1 2 179
2 1 208 2 1 178 2 1 187
2 2 146 2 2 145 2 2 141
3 1 186 3 1 183
3 2 142 3 2 125 3 2 136
;
ods output ParmSearch=parms;
proc mixed data=hh asycov mmeq mmeqsol covtest;
class a b;
model y = a / outp=predicted;
random b a*b;
lsmeans a;
parms (17 to 20 by .1) (.3 to .4 by .005) (1.0);
run;
proc print data=predicted;
run;
The ASYCOV option in the PROC statement requests the
asymptotic variance matrix of the covariance parameter
estimates. This matrix is the observed inverse Fisher
information matrix, which equals 2H-1,
where H is the Hessian matrix of the objective function
evaluated at the final covariance parameter estimates.
The MMEQ and MMEQSOL options in the PROC statement request
that the mixed model equations and their solution be
displayed.
The OUTP= option in the MODEL statement produces the
data set predicted, containing the predicted values. Least-squares
means (LSMEANS) are requested for A.
The PARMS and ODS statements are used to construct
a data set containing the likelihood surface.
The results from this analysis are shown in Output 41.3.1.
Output 41.3.1: Plotting the Likelihood
Model Information |
Data Set |
WORK.HH |
Dependent Variable |
y |
Covariance Structure |
Variance Components |
Estimation Method |
REML |
Residual Variance Method |
Profile |
Fixed Effects SE Method |
Model-Based |
Degrees of Freedom Method |
Containment |
|
The "Model Information" table lists details about this variance
components model.
Class Level Information |
Class |
Levels |
Values |
a |
3 |
1 2 3 |
b |
2 |
1 2 |
|
The "Class Level Information" table lists the levels for
A and B.
Dimensions |
Covariance Parameters |
3 |
Columns in X |
4 |
Columns in Z |
8 |
Subjects |
1 |
Max Obs Per Subject |
16 |
Observations Used |
16 |
Observations Not Used |
0 |
Total Observations |
16 |
|
The "Dimensions" table reveals that X is 16 ×4
and Z is 16 ×8. Since there are no SUBJECT= effects,
PROC MIXED considers the data effectively to be from one subject
with 16 observations.
Parameter Search |
CovP1 |
CovP2 |
CovP3 |
Variance |
Res Log Like |
-2 Res Log Like |
17.0000 |
0.3000 |
1.0000 |
80.1400 |
-52.4699 |
104.9399 |
17.0000 |
0.3050 |
1.0000 |
80.0466 |
-52.4697 |
104.9393 |
17.0000 |
0.3100 |
1.0000 |
79.9545 |
-52.4694 |
104.9388 |
17.0000 |
0.3150 |
1.0000 |
79.8637 |
-52.4692 |
104.9384 |
17.0000 |
0.3200 |
1.0000 |
79.7742 |
-52.4691 |
104.9381 |
17.0000 |
0.3250 |
1.0000 |
79.6859 |
-52.4690 |
104.9379 |
17.0000 |
0.3300 |
1.0000 |
79.5988 |
-52.4689 |
104.9378 |
17.0000 |
0.3350 |
1.0000 |
79.5129 |
-52.4689 |
104.9377 |
17.0000 |
0.3400 |
1.0000 |
79.4282 |
-52.4689 |
104.9377 |
17.0000 |
0.3450 |
1.0000 |
79.3447 |
-52.4689 |
104.9378 |
... | ... | ... | ... | ... | ... |
20.0000 |
0.3550 |
1.0000 |
78.2003 |
-52.4683 |
104.9366 |
20.0000 |
0.3600 |
1.0000 |
78.1201 |
-52.4684 |
104.9368 |
20.0000 |
0.3650 |
1.0000 |
78.0409 |
-52.4685 |
104.9370 |
20.0000 |
0.3700 |
1.0000 |
77.9628 |
-52.4687 |
104.9373 |
20.0000 |
0.3750 |
1.0000 |
77.8857 |
-52.4689 |
104.9377 |
20.0000 |
0.3800 |
1.0000 |
77.8096 |
-52.4691 |
104.9382 |
20.0000 |
0.3850 |
1.0000 |
77.7345 |
-52.4693 |
104.9387 |
20.0000 |
0.3900 |
1.0000 |
77.6603 |
-52.4696 |
104.9392 |
20.0000 |
0.3950 |
1.0000 |
77.5871 |
-52.4699 |
104.9399 |
20.0000 |
0.4000 |
1.0000 |
77.5148 |
-52.4703 |
104.9406 |
|
Only a portion of the "Parameter Search" table is shown because the
full listing has 651 rows.
Iteration History |
Iteration |
Evaluations |
-2 Res Log Like |
Criterion |
1 |
2 |
104.93416367 |
0.00000000 |
Convergence criteria met. |
|
Convergence is quick because PROC MIXED starts from the best
value from the grid search.
Covariance Parameter Estimates |
Cov Parm |
Estimate |
Standard Error |
Z Value |
Pr Z |
b |
1464.36 |
2098.01 |
0.70 |
0.2426 |
a*b |
26.9581 |
59.6570 |
0.45 |
0.3257 |
Residual |
78.8426 |
35.3512 |
2.23 |
0.0129 |
|
The preceding table lists the variance components estimates. Note that B is
much more variable than A*B.
Asymptotic Covariance Matrix of Estimates |
Row |
Cov Parm |
CovP1 |
CovP2 |
CovP3 |
1 |
b |
4401640 |
1.2831 |
-273.32 |
2 |
a*b |
1.2831 |
3558.96 |
-502.84 |
3 |
Residual |
-273.32 |
-502.84 |
1249.71 |
|
The asymptotic covariance matrix also reflects the large variability
of B relative to A*B.
Fit Statistics |
Res Log Likelihood |
-52.5 |
Akaike's Information Criterion |
-55.5 |
Schwarz's Bayesian Criterion |
-53.5 |
-2 Res Log Likelihood |
104.9 |
PARMS Model Likelihood Ratio Test |
DF |
Chi-Square |
Pr > ChiSq |
2 |
0.00 |
1.0000 |
|
The PARMS likelihood ratio test (LRT) compares the best
model from the grid search with the final fitted model.
Since these models are nearly the same, the LRT is
not significant.
Mixed Model Equations |
Row |
Effect |
a |
b |
Col1 |
Col2 |
Col3 |
Col4 |
Col5 |
Col6 |
Col7 |
Col8 |
Col9 |
Col10 |
Col11 |
Col12 |
Col13 |
1 |
Intercept |
|
|
0.2029 |
0.06342 |
0.07610 |
0.06342 |
0.1015 |
0.1015 |
0.03805 |
0.02537 |
0.03805 |
0.03805 |
0.02537 |
0.03805 |
36.4143 |
2 |
a |
1 |
|
0.06342 |
0.06342 |
|
|
0.03805 |
0.02537 |
0.03805 |
0.02537 |
|
|
|
|
13.8757 |
3 |
a |
2 |
|
0.07610 |
|
0.07610 |
|
0.03805 |
0.03805 |
|
|
0.03805 |
0.03805 |
|
|
12.7469 |
4 |
a |
3 |
|
0.06342 |
|
|
0.06342 |
0.02537 |
0.03805 |
|
|
|
|
0.02537 |
0.03805 |
9.7917 |
5 |
b |
|
1 |
0.1015 |
0.03805 |
0.03805 |
0.02537 |
0.1022 |
|
0.03805 |
|
0.03805 |
|
0.02537 |
|
21.2956 |
6 |
b |
|
2 |
0.1015 |
0.02537 |
0.03805 |
0.03805 |
|
0.1022 |
|
0.02537 |
|
0.03805 |
|
0.03805 |
15.1187 |
7 |
a*b |
1 |
1 |
0.03805 |
0.03805 |
|
|
0.03805 |
|
0.07515 |
|
|
|
|
|
9.3477 |
8 |
a*b |
1 |
2 |
0.02537 |
0.02537 |
|
|
|
0.02537 |
|
0.06246 |
|
|
|
|
4.5280 |
9 |
a*b |
2 |
1 |
0.03805 |
|
0.03805 |
|
0.03805 |
|
|
|
0.07515 |
|
|
|
7.2676 |
10 |
a*b |
2 |
2 |
0.03805 |
|
0.03805 |
|
|
0.03805 |
|
|
|
0.07515 |
|
|
5.4793 |
11 |
a*b |
3 |
1 |
0.02537 |
|
|
0.02537 |
0.02537 |
|
|
|
|
|
0.06246 |
|
4.6802 |
12 |
a*b |
3 |
2 |
0.03805 |
|
|
0.03805 |
|
0.03805 |
|
|
|
|
|
0.07515 |
5.1115 |
|
The mixed model equations are analogous to the normal equations in
the standard linear model. For this example, rows 1 -4
correspond to the fixed effects, rows 5 -12 correspond to the
random effects, and Col13 corresponds to the dependent
variable.
Mixed Model Equations Solution |
Row |
Effect |
a |
b |
Col1 |
Col2 |
Col3 |
Col4 |
Col5 |
Col6 |
Col7 |
Col8 |
Col9 |
Col10 |
Col11 |
Col12 |
Col13 |
1 |
Intercept |
|
|
761.84 |
-29.7718 |
-29.6578 |
|
-731.14 |
-733.22 |
-0.4680 |
0.4680 |
-0.5257 |
0.5257 |
-12.4663 |
-14.4918 |
159.61 |
2 |
a |
1 |
|
-29.7718 |
59.5436 |
29.7718 |
|
-2.0764 |
2.0764 |
-14.0239 |
-12.9342 |
1.0514 |
-1.0514 |
12.9342 |
14.0239 |
53.2049 |
3 |
a |
2 |
|
-29.6578 |
29.7718 |
56.2773 |
|
-1.0382 |
1.0382 |
0.4680 |
-0.4680 |
-12.9534 |
-14.0048 |
12.4663 |
14.4918 |
7.8856 |
4 |
a |
3 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
5 |
b |
|
1 |
-731.14 |
-2.0764 |
-1.0382 |
|
741.63 |
722.73 |
-4.2598 |
4.2598 |
-4.7855 |
4.7855 |
-4.2598 |
4.2598 |
26.8837 |
6 |
b |
|
2 |
-733.22 |
2.0764 |
1.0382 |
|
722.73 |
741.63 |
4.2598 |
-4.2598 |
4.7855 |
-4.7855 |
4.2598 |
-4.2598 |
-26.8837 |
7 |
a*b |
1 |
1 |
-0.4680 |
-14.0239 |
0.4680 |
|
-4.2598 |
4.2598 |
22.8027 |
4.1555 |
2.1570 |
-2.1570 |
1.9200 |
-1.9200 |
3.0198 |
8 |
a*b |
1 |
2 |
0.4680 |
-12.9342 |
-0.4680 |
|
4.2598 |
-4.2598 |
4.1555 |
22.8027 |
-2.1570 |
2.1570 |
-1.9200 |
1.9200 |
-3.0198 |
9 |
a*b |
2 |
1 |
-0.5257 |
1.0514 |
-12.9534 |
|
-4.7855 |
4.7855 |
2.1570 |
-2.1570 |
22.5560 |
4.4021 |
2.1570 |
-2.1570 |
-1.7134 |
10 |
a*b |
2 |
2 |
0.5257 |
-1.0514 |
-14.0048 |
|
4.7855 |
-4.7855 |
-2.1570 |
2.1570 |
4.4021 |
22.5560 |
-2.1570 |
2.1570 |
1.7134 |
11 |
a*b |
3 |
1 |
-12.4663 |
12.9342 |
12.4663 |
|
-4.2598 |
4.2598 |
1.9200 |
-1.9200 |
2.1570 |
-2.1570 |
22.8027 |
4.1555 |
-0.8115 |
12 |
a*b |
3 |
2 |
-14.4918 |
14.0239 |
14.4918 |
|
4.2598 |
-4.2598 |
-1.9200 |
1.9200 |
-2.1570 |
2.1570 |
4.1555 |
22.8027 |
0.8115 |
|
This solution matrix results from sweeping all but the last row of
the mixed model equations matrix. The final column contains a
solution vector for the fixed and random effects. The first four
rows correspond to fixed effects and the last eight to random
effects.
Type 3 Tests of Fixed Effects |
Effect |
Num DF |
Den DF |
F Value |
Pr > F |
a |
2 |
2 |
28.00 |
0.0345 |
|
The A factor is significant at the 5% level.
Least Squares Means |
Effect |
a |
Estimate |
Standard Error |
DF |
t Value |
Pr > |t| |
a |
1 |
212.82 |
27.6014 |
2 |
7.71 |
0.0164 |
a |
2 |
167.50 |
27.5463 |
2 |
6.08 |
0.0260 |
a |
3 |
159.61 |
27.6014 |
2 |
5.78 |
0.0286 |
|
The significance of A appears to be from the difference between
its first level and its other two levels.
Obs |
a |
b |
y |
Pred |
StdErrPred |
DF |
Alpha |
Lower |
Upper |
Resid |
1 |
1 |
1 |
237 |
242.723 |
4.72563 |
10 |
0.05 |
232.193 |
253.252 |
-5.7228 |
2 |
1 |
1 |
254 |
242.723 |
4.72563 |
10 |
0.05 |
232.193 |
253.252 |
11.2772 |
3 |
1 |
1 |
246 |
242.723 |
4.72563 |
10 |
0.05 |
232.193 |
253.252 |
3.2772 |
4 |
1 |
2 |
178 |
182.916 |
5.52589 |
10 |
0.05 |
170.603 |
195.228 |
-4.9159 |
5 |
1 |
2 |
179 |
182.916 |
5.52589 |
10 |
0.05 |
170.603 |
195.228 |
-3.9159 |
6 |
2 |
1 |
208 |
192.670 |
4.70076 |
10 |
0.05 |
182.196 |
203.144 |
15.3297 |
7 |
2 |
1 |
178 |
192.670 |
4.70076 |
10 |
0.05 |
182.196 |
203.144 |
-14.6703 |
8 |
2 |
1 |
187 |
192.670 |
4.70076 |
10 |
0.05 |
182.196 |
203.144 |
-5.6703 |
9 |
2 |
2 |
146 |
142.330 |
4.70076 |
10 |
0.05 |
131.856 |
152.804 |
3.6703 |
10 |
2 |
2 |
145 |
142.330 |
4.70076 |
10 |
0.05 |
131.856 |
152.804 |
2.6703 |
11 |
2 |
2 |
141 |
142.330 |
4.70076 |
10 |
0.05 |
131.856 |
152.804 |
-1.3297 |
12 |
3 |
1 |
186 |
185.687 |
5.52589 |
10 |
0.05 |
173.374 |
197.999 |
0.3134 |
13 |
3 |
1 |
183 |
185.687 |
5.52589 |
10 |
0.05 |
173.374 |
197.999 |
-2.6866 |
14 |
3 |
2 |
142 |
133.542 |
4.72563 |
10 |
0.05 |
123.013 |
144.072 |
8.4578 |
15 |
3 |
2 |
125 |
133.542 |
4.72563 |
10 |
0.05 |
123.013 |
144.072 |
-8.5422 |
16 |
3 |
2 |
136 |
133.542 |
4.72563 |
10 |
0.05 |
123.013 |
144.072 |
2.4578 |
|
The preceding output lists the predicted values from the model.
These values are the sum of the fixed effects estimates and the
empirical best linear unbiased predictors (EBLUPs) of the random
effects. It is often useful to plot predicted values and residuals
to assess the adequacy of the model, using another SAS procedure to
generate plots and diagnostic measures.
To plot the likelihood surface using the G3D procedure from
SAS/GRAPH software, use the following source:
proc g3d data=parms;
plot CovP1*CovP2 = ResLogLike
/ ctop=red cbottom=blue caxis=black;
run;
The results from this plot are shown in Output 41.3.2.
The peak of the surface is the REML estimates for the
B and A*B variance components.
Output 41.3.2: Plot of Likelihood Surface
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.