* ======================================================================= 
*  File:  	matrix 003.SPS .
*  Date:  	09-Dec-2003 .
*  Author:  	Bruce Weaver, weaverb@mcmaster.ca .
*  Notes:	Use matrix language to solve multiple regression equation .
* ======================================================================= .

* The following example is taken from "Statistics The Easy Way" (3rd Ed.)
* ISBN 0-8120-9392-5.

* -------------------------------------------------- .
* First, perform analysis using REGRESSION procedure .
* -------------------------------------------------- .

DATA LIST LIST /c x1 x2 y (4f5.0).
BEGIN DATA.
1  10   4   52
1  12   7   72
1  16   7   78
1  15  12   109
1  15   4   57
1  16  15   128
1   9  12   95
1  10  15   116
END DATA.

REGRESSION
  /MISSING LISTWISE
  /STATISTICS COEFF OUTS CI BCOV R ANOVA
  /CRITERIA=PIN(.05) POUT(.10)
  /NOORIGIN
  /DEPENDENT y
  /METHOD=ENTER x1 x2  .

* ------------------------------------------ .
* Now perform analysis using Matrix Language .
* ------------------------------------------ .

* NOTE:  Previous files (matrix 001 and 002) used COMPUTE
* statements to create matrices.  But in this case, I use
* GET statements to read them in from the working data file.

matrix.

GET y
 /file = *
 /var = y.
GET x
 /file = *
 /var = c x1 x2.
compute n = nrow(y).			/* N = number of rows in Y-matrix .
compute p = ncol(x)-1.			/* p = number of predictors = columns in X matrix less 1.
compute ybar = msum(y)/n.		/* Mean of Y values .
compute dev.y = y-ybar.			/* Matrix of (Y-Ybar) deviation scores.
compute ss.y = t(dev.y)*dev.y.		/* SS(Y) = dev'dev.
compute b = inv(t(x)*x) * t(x)*y.	/* B = INV(X'X)*X'Y .
compute yhat = x*b.			/* Yhat = XB .
compute e = y - yhat.			/* E = Y-Yhat .
compute ss.e = t(e)*e.			/* SSE = E'E .
compute ss.reg = ss.y - ss.e.		/* SS_regression .
compute ms.reg = ss.reg/p.		/* MS regression .
compute ms.e = ss.e/(n-p-1).		/* Variance of residuals .
compute cov.b = inv(t(x)*x)*ms.e.	/* (X'X)^-1*MS_error .
compute var.b = diag(cov.b).		/* Variances of B .
compute se.b = sqrt(var.b).		/* SE of B .
compute f = ms.reg/ms.e.		/* F-statistic .

print x /title = 'Matrix X'.
print y /title = 'Observed scores'.
print ybar /title = 'mean of Y'.
print yhat /title = 'Predicted scores'.
print e /title = 'Errors'.
print b /title = 'Matrix of coefficients (B)'.
print ss.e /title = 'Residual sum of squares'.
print cov.b /title = 'Covariance of B'.
print var.b /title = 'Variance of B'.
print b /title = 'Matrix of coefficients (B)'.
print se.b /title = 'SE of B'.
print ss.reg /title = 'SS Regression'.
print ss.e /title = 'SS Residual'.
print ss.y /title = 'SS Total'.
print ms.reg /title = 'MS Regression'.
print ms.e /title = 'MS Residual'.
print f /title = 'F-ratio'.

end matrix.

* ======================================================================= .