SPSSX Discussion

Solving for Parameter Estimates in Linear Mixed Model

Classic

List

Threaded

1 message

Ryan

Solving for Parameter Estimates in Linear Mixed Model

This post is intended for those interested in learning how to compute parameter estimates for a linear mixed model using matrix algebra. The data set was taken from this website:

http://ftp.sas.com/samples/A55235

Metal is the fixed effects variable, ingot is the random effects variable and pressure is the dependent
variable.

Hope this is helpful.

Best Wishes,

Ryan
--

***** Create data set *****.
DATA LIST /observation 1-2 ingot 4 metal 6 pressure 8-11(1).
BEGIN DATA
01 1 3 67.0
02 1 2 71.9
03 1 1 72.2
04 2 3 67.5
05 2 2 68.8
06 2 1 66.4
07 3 3 76.0
08 3 2 82.6
09 3 1 74.5
10 4 3 72.7
11 4 2 78.1
12 4 1 67.3
13 5 3 73.1
14 5 2 74.2
15 5 1 73.2
16 6 3 65.8
17 6 2 70.8
18 6 1 68.7
19 7 3 75.6
20 7 2 84.9
21 7 1 69.0
END DATA.

***** random intercept model *****.
MIXED pressure BY metal
/FIXED=metal | SSTYPE(3)
/METHOD=REML
/PRINT=CORB COVB DESCRIPTIVES G LMATRIX R SOLUTION TESTCOV
/RANDOM=INTERCEPT | SUBJECT(ingot) COVTYPE(VC).

***** repeated measures model with CS var-cov residual matrix *****.
MIXED pressure BY metal
/FIXED=metal | SSTYPE(3)
/METHOD=REML
/PRINT=CORB COVB DESCRIPTIVES LMATRIX R SOLUTION TESTCOV
/REPEATED=metal | SUBJECT(ingot) COVTYPE(CS).

***** beta hats for the random intercept model using matrix algebra *****.
matrix.

*****Y vector *****.
COMPUTE Y = {67.0;
71.9;
72.2;
67.5;
68.8;
66.4;
76.0;
82.6;
74.5;
72.7;
78.1;
67.3;
73.1;
74.2;
73.2;
65.8;
70.8;
68.7;
75.6;
84.9;
69.0}.

***** X is a 21 X 4 fixed effects design matrix *****.
COMPUTE X = {1,1,0,0;
1,0,1,0;
1,0,0,1;
1,1,0,0;
1,0,1,0;
1,0,0,1;
1,1,0,0;
1,0,1,0;
1,0,0,1;
1,1,0,0;
1,0,1,0;
1,0,0,1;
1,1,0,0;
1,0,1,0;
1,0,0,1;
1,1,0,0;
1,0,1,0;
1,0,0,1;
1,1,0,0;
1,0,1,0;
1,0,0,1}.

***** Z is a 21 X 7 random effects design matrix *****.
COMPUTE Z = {1,0,0,0,0,0,0;
0,1,0,0,0,0,0;
0,0,1,0,0,0,0;
0,0,0,1,0,0,0;
0,0,0,0,1,0,0;
0,0,0,0,0,1,0;
0,0,0,0,0,0,1;
1,0,0,0,0,0,0;
0,1,0,0,0,0,0;
0,0,1,0,0,0,0;
0,0,0,1,0,0,0;
0,0,0,0,1,0,0;
0,0,0,0,0,1,0;
0,0,0,0,0,0,1;
1,0,0,0,0,0,0;
0,1,0,0,0,0,0;
0,0,1,0,0,0,0;
0,0,0,1,0,0,0;
0,0,0,0,1,0,0;
0,0,0,0,0,1,0;
0,0,0,0,0,0,1}.

***** G is a 7 X 7 random effects variance-covariance matrix ******.
***** Note: 11.4477... is the variance of the intercepts of the 7 ingots ******.
COMPUTE G = {11.447777777778326,0,0,0,0,0,0;
0,11.447777777778326,0,0,0,0,0;
0,0,11.447777777778326,0,0,0,0;
0,0,0,11.447777777778326,0,0,0;
0,0,0,0,11.447777777778326,0,0;
0,0,0,0,0,11.447777777778326,0;
0,0,0,0,0,0,11.447777777778326}.

***** R is a 21 X 21 residual variance covariance matrix ******.
***** Note: 10.37... is the residual the variance in Y ******.
COMPUTE R = {10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574,0;
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10.371587301585574}.

***** variance-covariance matrix *****.
COMPUTE V = Z*G*TRANSPOS(Z)+R.

***** beta hats vector *****.

COMPUTE B = GINV(TRANSPOS(X)*INV(V)*X)*TRANSPOS(X)*INV(V)*Y.

***** print the results ******.
print Y /title "Matrix Y".
print X /title "Matrix X".
print Z /title "Matrix Z".
print G /title "Matrix G".
print R /title "Matrix R".
print V /title "Matrix V".
print B /title "beta hats".

end matrix.

***** beta hats *****.
COMPUTE b0 = 54.29642857 + 16.80357143.
EXECUTE.
COMPUTE b1 = (54.29642857 + 15.88928571) - (54.29642857 + 16.80357143).
EXECUTE.
COMPUTE b2 = (54.29642857 + 21.60357143) - (54.29642857 + 16.80357143) .
EXECUTE.