/* Filename: splnfit.sas Purpose: fit four linear regression models (splines) on Nepal data*/

data a; infile "c:\nepal.dat";
 input id age gender wt ht arm bf day month year mage lt;
run;

/*delete the irrlevant variables and the missing values*/
data b; set a (keep=age wt);
  if wt>80 then delete;
run;

/*find the 10th, 33rd, 67th, and 90th percentile*/
proc univariate noprint;
 var age;
 output out=pctout pctlpts=10 33 67 90 pctlpre=p pctlname=p10 p33 p67 p90;
quit;

/*look at the percentiles*/
proc print data=pctout; quit;

/*create the variables needed for fitting the data using cubic polynomial curve and splines*/
data c; set b;

/*create variables for fitting various splines wiith knots at 10th, 33rd 67th and 90th
 percentiles. The values 14, 27, 48, and 63 are obtained in the PROC UNIVARIATE */

 if age<=14 then age14=0; else age14=age-14;
 if age<=27 then age27=0; else age27=age-27;
 if age<=48 then age48=0; else age48=age-48;
 if age<=63 then age63=0; else age63=age-63;
 ages=age*age;
 ageq=age*age*age;
 ageq14=age14*age14*age14;
 ageq27=age27*age27*age27;
 ageq48=age48*age48*age48;
 ageq63=age63*age63*age63;
 term2=ageq14-ageq48*((63-14)/(63-48))+ageq63*((48-14)/(63-48));
 term3=ageq27-ageq48*((63-27)/(63-48))+ageq63*((48-27)/(63-48));
run;

/*Fitting cubic polynomial*/
proc glm;
 model wt=age ages ageq;
 output out=out1 p=pred1;
run;

/*Fitting linear spline with two knots*/
proc glm;
 model wt=age age27 age48;
 output out=out2 p=pred2;
run;

/*Fitting cubic spline with two knots*/
proc glm;
 model wt=age ages ageq ageq27 ageq48;
 output out=out3 p=pred3;
run;

/*Fitting cubic spline with four knots*/
proc glm;
 model wt=age term2 term3;
 output out=out4 p=pred4;
run;

/*create the data set for plotting*/
data plotdat; merge out1 out2 out3 out4; keep age wt pred1 pred2 pred3 pred4; run;

/*sort the data by age for the purpose of drawing fitted line*/
proc sort data=plotdat; by age; run;

goptions reset=all device=win target=hplj3si vpos=40 hpos=100 hsize=8 in vsize=6 in
   rotate=landscape vorigin=1.5 in horigin=1.5 in;
 symbol1 interpol=none value=point;
 axis1 minor=none order=(0 to 20 by 2) label=(a=90 'weight');
 axis2 minor=none label=('age in weeks');
title1 'Scatter Plot of Weight Versus Age';
proc gplot data=plotdat;
 plot wt*age/nolegend frame haxis=axis2 vaxis=axis1;
run;
title1 'Predicted Weight Versus Age (Cubic Polynomial)';
 symbol1 interpol=join;
 plot pred1*age/nolegend frame haxis=axis2 vaxis=axis1;
run;
title1 'Predicted Weight Versus Age (Linear Spline With Two Knots)';
 plot pred2*age/nolegend frame haxis=axis2 vaxis=axis1; run;
title1 'Predicted Weight Versus Age (Cubic Spline With Three Knots)';
 plot pred3*age/nolegend frame haxis=axis2 vaxis=axis1; run;
title1 'Predicted Weight Versus Age (Cubic Spline With Four Knots)';
 plot pred4*age/nolegend frame haxis=axis2 vaxis=axis1; run;
quit;