/* Use the NLIN procedure in SAS to do the analysis for problem 2 on assignment 9, Spring 2002 */ data set1; infile 'c:\courses\st511\snew\pnonlin.dat'; input y x1 x2; run; proc print data=set1; run; /* Get starting values */ data set2; set set1; logy = log(y); logx1 = log(x1); logx2 = log(x2); run; proc reg data=set2 outest=set3; model logy = logx1 logx2 ; run; proc print data=set3; run; /* Add an extra case to the data file to get an estimated mean for x1=50 and x2=60 and 95% coonfidence interval for the mean response and a 95% prediction interval */ data temp; input y x1 x2; cards; . 50 60 run; data set1a; set set1 temp; run; /* Use the NLIN procedure to fit the model Here we include formulas for first partial derivatives and use the Gauss-Newton procedure */ proc nlin data=set1a method=gauss convergeparm=.000001 maxiter=50 outest=setp ; parms b0 = 9.591 b1=0.51485 b2=0.29845; model y = b0*(x1**b1)*(x2**b2); der.bo = (x1**b1)*(x2**b2); der.b1 = b0*(x1**b1)*(x2**b2)*log(x1); der.b2 = b0*(x1**b1)*(x2**b2)*log(x2); output out=set4 residual=r predicted=yhat student=student stdi=stdi stdp=stdp l95=l95pred u95=u95pred l95m=l95mean u95m=u95mean; run; proc print data=set4; var y x1 x2 r yhat l95pred u95pred l95mean u95mean; run; /* Make residual plots */ proc rank data=set4 normal=blom out=set4; var r; ranks q; run; /* UNIX users can replace the following goptions line with filename graffile pipe 'lpr -Dpostscript'; goptions gsfmode=replace gsfname=graffile targetdevice=ps300 rotate=landscape; */ goptions device=WIN target=WINPRTC rotate=landscape; axis1 label=(f=swiss h=1.8 "Temperature(C)") value=(f=swiss h=1.8) w=3.0 length= 6.5 in; axis2 label=(f=swiss h=1.8 "Pressure(psi)") value=(f=swiss h=1.8) w=3.0 length= 6.5 in; axis3 label=(f=swiss h=2.0 a=90 r=0 "Residuals") value=(f=swiss h=1.8) w= 3.0 length = 5.0 in; axis4 label=(f=swiss h=1.8 "q") order = -3 to 3 by 1 value=(f=swiss h=1.8) w=3.0 length= 6.5 in; SYMBOL1 V=CIRCLE H=1.7 w=3 l=1 i=none ; SYMBOL2 V=DIAMOND H=1.7 w=3 l=3 i=join ; SYMBOL3 V=square H=1.7 w=3 l=9 i=join ; PROC GPLOT DATA=set4; PLOT r*x1 / vaxis=axis3 haxis=axis1; TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Temperature"; footnote ls=0.01in; RUN; PROC GPLOT DATA=set4; PLOT r*x2 / vaxis=axis3 haxis=axis2; TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Pressure"; footnote ls=0.01in; RUN; PROC GPLOT DATA=set4; PLOT r*q / vaxis=axis3 haxis=axis4; TITLE1 ls=0.01in H=2.0 F=swiss "Normal Probability Plot"; footnote ls=0.01in; RUN; /* Use the estimated parameters and the covariance matrix produced by the outest= option in the nlin procedure to obtain an approximate t-test for beta1=beta2 */ proc print data=setp; run; data setpp; set setp; if(_TYPE_="FINAL" or _TYPE_="COVB"); keep b0 b1 b2; run; proc print data=setpp; run; proc iml; start b1vsb2; use setpp; read all into x; use set1; read all into w; nn=nrow(w); np=ncol(x); nr=nrow(x); b = x[1, ]; vcov = x[2:nr, ]; c={0 1 -1}; t=c*t(b)/sqrt(c*vcov*t(c)); at=abs(t); df=nn-np; pvalue=(1-probt(at,df))*2; print t df pvalue; finish b1vsb2; run b1vsb2; /* Another way to test the hypothesis that beta1=beta2 is to fit a model with the contraints that beta1=beta2 and perform an approximate f-test f=(SSE(restricted)-SSE(unrestricted))/(df)/(MSE(unrestricted)) */ data set1; set set1; x12=x1*x2; run; proc nlin data=set1 method=gauss convergeparm=.000001 maxiter=50 outest=setpr; parms b0 = 10. b1=0.4 ; model y = b0*(x12**b1); der.bo = (x12**b1); der.b1 = b0*(x12**b1)*log(x12); run; proc print data=setpr; run; data setpr2; set setpr; if(_TYPE_="FINAL"); sser=_SSE_; keep sser; run; data setp2; set setp; if(_TYPE_="FINAL"); sse=_SSE_; keep sse; run; data both; merge setp2 setpr2; run; proc iml; start ff; use both; read all into sse; f = (sse[1,2]-sse[1,1])/(sse[1,1]/15); pvalue = 1-probf(f,1,15); print sse f pvalue; finish; run; run ff;