The MIXED Procedure

Example 78.3 Plotting the Likelihood

(View the complete code for this example.)

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 , where 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 78.3.1Output 78.3.13.

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

Output 78.3.1: Model Information

The Mixed Procedure

Model Information
Data SetWORK.HH
Dependent Variabley
Covariance StructureVariance Components
Estimation MethodREML
Residual Variance MethodProfile
Fixed Effects SE MethodModel-Based
Degrees of Freedom MethodContainment


The "Class Level Information" table in Output 78.3.2 lists the levels for A and B.

Output 78.3.2: Class Level Information

Class Level Information
ClassLevelsValues
a31 2 3
b21 2


The "Dimensions" table in Output 78.3.3 reveals that is and is . Since there are no SUBJECT= effects, PROC MIXED considers the data to be effectively from one subject with 16 observations.

Output 78.3.3: Model Dimensions and Number of Observations

Dimensions
Covariance Parameters3
Columns in X4
Columns in Z8
Subjects1
Max Obs per Subject16

Number of Observations
Number of Observations Read16
Number of Observations Used16
Number of Observations Not Used0


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

Output 78.3.4: Selected Results of Parameter Search

The Mixed Procedure

CovP1CovP2CovP3VarianceRes Log Like-2 Res Log Like
17.00000.30001.000080.1400-52.4699104.9399
17.00000.30501.000080.0466-52.4697104.9393
17.00000.31001.000079.9545-52.4694104.9388
17.00000.31501.000079.8637-52.4692104.9384
17.00000.32001.000079.7742-52.4691104.9381
17.00000.32501.000079.6859-52.4690104.9379
17.00000.33001.000079.5988-52.4689104.9378
17.00000.33501.000079.5129-52.4689104.9377
17.00000.34001.000079.4282-52.4689104.9377
17.00000.34501.000079.3447-52.4689104.9378
......
......
......
20.00000.35501.000078.2003-52.4683104.9366
20.00000.36001.000078.1201-52.4684104.9368
20.00000.36501.000078.0409-52.4685104.9370
20.00000.37001.000077.9628-52.4687104.9373
20.00000.37501.000077.8857-52.4689104.9377
20.00000.38001.000077.8096-52.4691104.9382
20.00000.38501.000077.7345-52.4693104.9387
20.00000.39001.000077.6603-52.4696104.9392
20.00000.39501.000077.5871-52.4699104.9399
20.00000.40001.000077.5148-52.4703104.9406


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

Output 78.3.5: Iteration History and Convergence Status

Iteration History
IterationEvaluations-2 Res Log LikeCriterion
12104.934163670.00000000

Convergence criteria met.


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

Output 78.3.6: Estimated Covariance Parameters

Covariance Parameter Estimates
Cov ParmEstimateStandard
Error
Z ValuePr > Z
b1464.362098.010.700.2426
a*b26.958159.65700.450.3257
Residual78.842635.35122.230.0129


The asymptotic covariance matrix in Output 78.3.7 also reflects the large variability of B relative to A*B.

Output 78.3.7: Asymptotic Covariance Matrix of Covariance Parameters

Asymptotic Covariance Matrix of Estimates
RowCov ParmCovP1CovP2CovP3
1b44016401.2831-273.32
2a*b1.28313558.96-502.84
3Residual-273.32-502.841249.71


As Output 78.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 78.3.8: Fit Statistics and Likelihood Ratio Test

Fit Statistics
-2 Res Log Likelihood104.9
AIC (Smaller is Better)110.9
AICC (Smaller is Better)113.6
BIC (Smaller is Better)107.0

PARMS Model Likelihood Ratio Test
DFChi-SquarePr > ChiSq
20.001.0000


The mixed model equations are analogous to the normal equations in the standard linear model. As Output 78.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 78.3.9: Mixed Model Equations

Mixed Model Equations
RowEffectabCol1Col2Col3Col4Col5Col6Col7Col8Col9Col10Col11Col12Col13
1Intercept  0.20290.063420.076100.063420.10150.10150.038050.025370.038050.038050.025370.0380536.4143
2a1 0.063420.06342  0.038050.025370.038050.02537    13.8757
3a2 0.07610 0.07610 0.038050.03805  0.038050.03805  12.7469
4a3 0.06342  0.063420.025370.03805    0.025370.038059.7917
5b 10.10150.038050.038050.025370.1022 0.03805 0.03805 0.02537 21.2956
6b 20.10150.025370.038050.03805 0.1022 0.02537 0.03805 0.0380515.1187
7a*b110.038050.03805  0.03805 0.07515     9.3477
8a*b120.025370.02537   0.02537 0.06246    4.5280
9a*b210.03805 0.03805 0.03805   0.07515   7.2676
10a*b220.03805 0.03805  0.03805   0.07515  5.4793
11a*b310.02537  0.025370.02537     0.06246 4.6802
12a*b320.03805  0.03805 0.03805     0.075155.1115


The solution matrix in Output 78.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 78.3.10: Solutions of the Mixed Model Equations

Mixed Model Equations Solution
RowEffectabCol1Col2Col3Col4Col5Col6Col7Col8Col9Col10Col11Col12Col13
1Intercept  761.84-29.7718-29.6578 -731.14-733.22-0.46800.4680-0.52570.5257-12.4663-14.4918159.61
2a1 -29.771859.543629.7718 -2.07642.0764-14.0239-12.93421.0514-1.051412.934214.023953.2049
3a2 -29.657829.771856.2773 -1.03821.03820.4680-0.4680-12.9534-14.004812.466314.49187.8856
4a3              
5b 1-731.14-2.0764-1.0382 741.63722.73-4.25984.2598-4.78554.7855-4.25984.259826.8837
6b 2-733.222.07641.0382 722.73741.634.2598-4.25984.7855-4.78554.2598-4.2598-26.8837
7a*b11-0.4680-14.02390.4680 -4.25984.259822.80274.15552.1570-2.15701.9200-1.92003.0198
8a*b120.4680-12.9342-0.4680 4.2598-4.25984.155522.8027-2.15702.1570-1.92001.9200-3.0198
9a*b21-0.52571.0514-12.9534 -4.78554.78552.1570-2.157022.55604.40212.1570-2.1570-1.7134
10a*b220.5257-1.0514-14.0048 4.7855-4.7855-2.15702.15704.402122.5560-2.15702.15701.7134
11a*b31-12.466312.934212.4663 -4.25984.25981.9200-1.92002.1570-2.157022.80274.1555-0.8115
12a*b32-14.491814.023914.4918 4.2598-4.2598-1.92001.9200-2.15702.15704.155522.80270.8115


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

Output 78.3.11: Tests of Fixed Effects

Type 3 Tests of Fixed Effects
EffectNum DFDen DFF ValuePr > F
a2228.000.0345


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

Output 78.3.12: Least Squares Means for A Effect

Least Squares Means
EffectaEstimateStandard
Error
DFt ValuePr > |t|
a1212.8227.601427.710.0164
a2167.5027.546326.080.0260
a3159.6127.601425.780.0286


Output 78.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 78.3.13: Predicted Values

ObsabyPredStdErrPredDFAlphaLowerUpperResid
111237242.7234.72563100.05232.193253.252-5.7228
211254242.7234.72563100.05232.193253.25211.2772
311246242.7234.72563100.05232.193253.2523.2772
412178182.9165.52589100.05170.603195.228-4.9159
512179182.9165.52589100.05170.603195.228-3.9159
621208192.6704.70076100.05182.196203.14415.3297
721178192.6704.70076100.05182.196203.144-14.6703
821187192.6704.70076100.05182.196203.144-5.6703
922146142.3304.70076100.05131.856152.8043.6703
1022145142.3304.70076100.05131.856152.8042.6703
1122141142.3304.70076100.05131.856152.804-1.3297
1231186185.6875.52589100.05173.374197.9990.3134
1331183185.6875.52589100.05173.374197.999-2.6866
1432142133.5424.72563100.05123.013144.0728.4578
1532125133.5424.72563100.05123.013144.072-8.5422
1632136133.5424.72563100.05123.013144.0722.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 78.3.14. The peak of the surface is the REML estimates for the B and A*B variance components.

Output 78.3.14: Plot of Likelihood Surface

Plot of Likelihood Surface