/* Use the NLIN procedure in SAS to do the analysis for problem 3 on assignment 9, Spring 2002 */ data set1; infile 'c:\courses\st511\snew\kinetics.dat'; input y x; run; proc print data=set1; run; /* Get starting values */ data set2; set set1; yinv = 1/y; xinv = 1/x; run; proc reg data=set2 outest=set3; model yinv = xinv ; run; proc print data=set3; run; /* Add an extra case to the data file to get an estimated mean for x=25 and 95% coonfidence interval for the mean response and a 95% prediction interval */ data temp; input y x; cards; . 25 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 = 29.62 b1=13.45; model y = (b0*x)/(b1+x); der.bo = x/(b1+x); der.b1 = -(b0*x)/((b1+x)**2); 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 x 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 "Concentration") 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*x / vaxis=axis3 haxis=axis1; TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Concentration"; 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; ;