/****************************************************************/ /* NAME: IMLREG1 */ /* TITLE: Regression with IML */ /* DATA: THERAPY data */ /****************************************************************/ data THERAPY; input NAME $ SEX $ PERSTEST THERAPY INTEXT; If SEX = 'F' then SX = 1; /* Dummy variable for sex */ else SX = 0; cards; John M 26 32 3 Susan F 24 40 4 Mary F 22 44 8 Paul M 33 44 4 Jenny F 27 48 6 Rick M 36 52 4 Cathy F 30 56 10 Robert M 38 56 4 Lisa F 30 60 12 Tina F 34 68 15 ; proc iml; /*---------REGRESSION ROUTINE------------*/ /* given X, and y, this fits y = X b + e */ /* by least squares. */ start reg; n = nrow(X); /* number of observations */ k = ncol(X); /* number of variables */ xpx = t(X)*X; /* cross-products, X'X */ xpy = t(X)*y; xpxi = inv(xpx); /* inverse crossproducts */ b = xpxi*xpy; /* parameter estimates */ yhat = X*b; /* predicted values */ resid = y-yhat; /* residuals */ sse = t(resid)*resid; /* sum of squared errors */ dfe = n-k; /* degrees of freedom error */ mse = sse/dfe; /* mean squared error */ rmse = sqrt(mse); /* root mean squared error */ page; covb = xpxi#mse; /* covariance of estimates */ stdb = sqrt(vecdiag(covb)); /* standard errors */ t = b/stdb; /* ttest for estimates=0 */ probt = 1-probf(t#t,1,dfe); /* significance probability */ print vnames b stdb t[f=7.3] probt[f=7.4]; s = diag(1/stdb); corrb = s*covb*s; /* correlation of estimates */ print "Covariance of Estimates", covb[r=vnames c=vnames], "Correlation of Estimates",corrb[r=vnames c=vnames f=7.3]; if nrow(alpha)= 0 then return; /* is alpha specified? */ H = X*xpxi*t(X); /* hat matrix */ vresid = (I(n)-H)*mse; /* covariance of residuals */ vpred = H#mse; /* covariance of yhat values */ hat= vecdiag(H); /* hat leverage values */ alpha=.05; t = tinv(1-alpha/2,dfe); lowerm=yhat-t#sqrt(hat*mse); /* lower conf. limit for mean */ upperm=yhat+t#sqrt(hat*mse); /* upper */ lower =yhat-t#sqrt(hat*mse+mse); /* lower conf. limit for indiv */ upper =yhat+t#sqrt(hat*mse+mse); /* upper */ print "Predicted Values, Residuals, Leverage, and Conf. Limits" , id y yhat resid[f=6.2] hat[f=6.2] lowerm[f=6.2] upperm[f=6.2] /* lower[f=6.2] upper[f=6.2] */ ; finish; /*---routine to test a linear combination of the estimates---*/ /* given L, this routine tests hypothesis that L b = 0. */ start test; dfh=nrow(L); /* degrees of freedom of test */ Lb=L*b; /* lin. comb of coeffs tested */ vlb=L*xpxi*t(L); /* variance of LB */ q=t(Lb)*inv(vlb)*Lb /dfh; /* Mean square for H0: LB=0 */ f=q/mse; prob=1-probf(f,dfh,dfe); print l[c=vnames], lb f dfh dfe prob; finish; /* Read the SAS data set THERAPY into IML variables */ use therapy; read all var{THERAPY} into y; read all var{PERSTEST INTEXT SX} into X[r=NAME c=VNAMES]; X = J( nrow(X),1) || X; /* add column of 1s */ vnames = 'Intercept' // t(vnames); id = name; alpha=.05; /* Error rate for confidence intervals */ run reg; print "Test Coef for X1"; L= {0 1 0 0}; run test; * print "TEST Coef for X1, X2"; L= {0 1 0 0,0 0 1 0} ; *run test; quit;