dm 'log;clear;output;clear;'; options ps=50 ls=70 pageno=1; goptions reset=global border ftext=swiss gunit=cm htext=0.4 htitle=0.5; goptions display noprompt; *********************************************************************; ** **; ** AUTHOR: Chris Bilder **; ** COURSE: STAT 4043 **; ** DATE: 11-17-00 **; ** UPDATE: **; ** PURPOSE: Do loess example in NKNW p.426-7 **; ** NOTES: **; ** **; *********************************************************************; title1 'Chris Bilder, STAT 4043'; data set1; infile 'c:/chris/osu/stat4043/NKNW_data/ch10ta06.dat'; input X1 X2 Y; format X1 d4.0; format X2 d4.1; format Y d4.2; run; *Create scatter plot of the data; proc g3d data=set1; title2 '3D scatter plot'; scatter X1*X2=Y / shape='cube' grid color='blue'; run; data pred_points; do X1 = 28 to 78 by 2; do X2 = 1 to 10 by 0.5; output; end; end; format X1 d4.0; format X2 d4.1; run; title2 "Loess method"; proc loess data=set1; model Y = X1 X2 / degree=1 details smooth=0.5 clm scale=SD std t residual; score data=pred_points / clm ; ods Output ScoreResults=ScoreOut OutputStatistics=StatOut; run; title2 'Predicted value at X_h1=30 and X_h2=3'; proc print data=scoreout; where X1=30 and X2=3; run; proc g3d data=scoreout; format X1 f2.1; title2 '3D surface, q=0.5'; plot X1*X2=p_y / grid xticknum=5 yticknum=6 zmin=-100 zmax=500 zticknum=7; run; proc gcontour data=scoreout; plot X2*X1=p_y / haxis=axis1 vaxis=axis2 legend=legend1 grid levels=50 100 150 200 250 300 350 400 autolabel=(reveal); title2 'Contour plot'; axis1 label = ('X1') length=9 order = (30 to 75 by 9); axis2 label=(a=90 'X2') length=9 order = (3 to 10 by 1); legend1 label=('Predicted Y') down=2; run; ***************************************; ** Usual multiple linear regression **; ***************************************; data extra; input X1 X2 Y; datalines; 30 3 . ; run; data set1; set set1 extra; run; title2 'Multiple linear regression model'; proc reg data=set1 outest=beta_hat; model Y = X1 X2 / clm; output out=out_set1 p=predict; symbol1 v=dot h=.2 cv=blue; *plot residual.*X1 / nostat vref=0 cvref=red; *plot residual.*X2 / nostat vref=0 cvref=red; *plot residual.*predicted. / nostat vref=0; *plot student.*predicted. / nostat vref=0 cvref=red; run; proc score data=pred_points score=beta_hat out=Y_hat predict type=parms; var X1 X2; run; proc g3d data=Y_hat; title2 '3D surface'; plot X1*X2=model1 / grid xticknum=5 yticknum=6 zmin=-100 zmax=500 zticknum=7; run; *Create contour plot; proc gcontour data=Y_hat; plot X2*X1=model1 / haxis=axis1 vaxis=axis2 legend=legend1 grid levels=50 100 150 200 250 300 350 400 autolabel=(reveal); title2 'Contour plot'; axis1 label = ('X1') length=9 order = (30 to 75 by 9); axis2 label=(a=90 'X2') length=9 order = (3 to 10 by 1); legend1 label=('Predicted Y') down=2; run; data residual2; set out_set1; residual2 = Y - predict; run; data residual3; merge residual2 statout(rename=(residual=residual1)); residual2 = abs(residual2); residual1 = abs(residual1); run; proc gplot data=residual3; plot residual2*obs residual1*obs / overlay haxis=axis1 vaxis=axis2 legend=legend1 autohref lhref=35; title2 'Compare absolute value of residuals'; axis1 label = ('Observation'); axis2 label=(a=90 '|Residual|'); symbol1 v=dot h=.4 cv=blue; symbol2 v=square h=.4 cv=red; legend1 label = none position = (bottom center outside) across = 2 down = 1 mode = reserve frame offset =(0.5, 0.5) value = (j=l h=0.3 'Multiple linear reg. residual' 'Loess residual'); run; quit;