The MIXED Procedure

Example 65.3 Plotting the Likelihood

The data for this example are from Hemmerle and Hartley (1973) and are also used for an example in 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 statements are 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 MIXED statement requests the asymptotic variance matrix of the covariance parameter estimates. This matrix is the observed inverse Fisher information matrix, which equals $2\mb{H} ^{-1}$, where $\mb{H} $ is the Hessian matrix of the objective function evaluated at the final covariance parameter estimates. The MMEQ and MMEQSOL options in the PROC MIXED 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 65.3.1Output 65.3.13.

The "Model Information" table in Output 65.3.1 lists details about this variance components model.

Output 65.3.1: Model Information

The Mixed Procedure

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 "Class Level Information" table in Output 65.3.2 lists the levels for A and B.

Output 65.3.2: Class Level Information

Class Level Information
Class Levels Values
a 3 1 2 3
b 2 1 2



The "Dimensions" table in Output 65.3.3 reveals that $\mb{X}$ is $16 \times 4$ and $\mb{Z}$ is $16 \times 8$. Since there are no SUBJECT= effects, PROC MIXED considers the data to be effectively from one subject with 16 observations.

Output 65.3.3: Model Dimensions and Number of Observations

Dimensions
Covariance Parameters 3
Columns in X 4
Columns in Z 8
Subjects 1
Max Obs per Subject 16

Number of Observations
Number of Observations Read 16
Number of Observations Used 16
Number of Observations Not Used 0



Only a portion of the "Parameter Search" table is shown in Output 65.3.4 because the full listing has 651 rows.

Output 65.3.4: Selected Results of Parameter Search

The Mixed Procedure

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



As Output 65.3.5 shows, convergence occurs quickly because PROC MIXED starts from the best value from the grid search.

Output 65.3.5: Iteration History and Convergence Status

Iteration History
Iteration Evaluations -2 Res Log Like Criterion
1 2 104.93416367 0.00000000

Convergence criteria met.



The "Covariance Parameter Estimates" table in Output 65.3.6 lists the variance components estimates. Note that B is much more variable than A*B.

Output 65.3.6: Estimated Covariance Parameters

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 asymptotic covariance matrix in Output 65.3.7 also reflects the large variability of B relative to A*B.

Output 65.3.7: Asymptotic Covariance Matrix of Covariance Parameters

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



As Output 65.3.8 shows, 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.

Output 65.3.8: Fit Statistics and Likelihood Ratio Test

Fit Statistics
-2 Res Log Likelihood 104.9
AIC (Smaller is Better) 110.9
AICC (Smaller is Better) 113.6
BIC (Smaller is Better) 107.0

PARMS Model Likelihood Ratio Test
DF Chi-Square Pr > ChiSq
2 0.00 1.0000



The mixed model equations are analogous to the normal equations in the standard linear model. As Output 65.3.9 shows, 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.

Output 65.3.9: Mixed Model Equations

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 solution matrix in Output 65.3.10 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 correspond to random effects.

Output 65.3.10: Solutions of the Mixed Model Equations

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



The A factor is significant at the 5% level (Output 65.3.11).

Output 65.3.11: Tests of Fixed Effects

Type 3 Tests of Fixed Effects
Effect Num DF Den DF F Value Pr > F
a 2 2 28.00 0.0345



Output 65.3.12 shows that the significance of A appears to be from the difference between its first level and its other two levels.

Output 65.3.12: Least Squares Means for A Effect

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



Output 65.3.13 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.

Output 65.3.13: Predicted Values

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



To plot the likelihood surface by using ODS Graphics, use the following statements:

proc template;
   define statgraph surface;
      begingraph;
         layout overlay3d;
            surfaceplotparm x=CovP1 y=CovP2 z=ResLogLike;
         endlayout;
      endgraph;
   end;
run;
proc sgrender data=parms template=surface;
run;

The results from this plot are shown in Output 65.3.14. The peak of the surface is the REML estimates for the B and A*B variance components.

Output 65.3.14: Plot of Likelihood Surface

Plot of Likelihood Surface