PROC MIXED: Plotting the Likelihood :: SAS/STAT(R) 9.22 User's Guide

The MIXED Procedure

Example 56.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 $\text{[math]}$ , where $\text{[math]}$ 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 56.3.1–Output 56.3.13.

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

Output 56.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 56.3.2 lists the levels for A and B.

Output 56.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 56.3.3 reveals that $\text{[math]}$ is $\text{[math]}$ and $\text{[math]}$ is $\text{[math]}$ . Since there are no SUBJECT= effects, PROC MIXED considers the data to be effectively from one subject with 16 observations.

Output 56.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 56.3.4 because the full listing has 651 rows.

Output 56.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 56.3.5 shows, convergence occurs quickly because PROC MIXED starts from the best value from the grid search.

Output 56.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 56.3.6 lists the variance components estimates. Note that B is much more variable than A*B.

Output 56.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 56.3.7 also reflects the large variability of B relative to A*B.

Output 56.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 56.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 56.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 56.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 56.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 56.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 56.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 56.3.11).

Output 56.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 56.3.12 shows that the significance of A appears to be from the difference between its first level and its other two levels.

Output 56.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 56.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 56.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 56.3.14. The peak of the surface is the REML estimates for the B and A*B variance components.

Output 56.3.14 Plot of Likelihood Surface

Top of Page