/* Polynomial regression on the Rabbit Jawbone Data (from book by Freund and Wilson) */


data rabbit;
input age length;
cards;
0.01 15.5
0.2 26.1
0.2 26.3
0.21 26.7
0.23 27.5
0.24 27
0.24 27
0.25 26
0.26 28.6
0.34 29.8
0.41 29.7
0.83 37.7
1.09 41.5
1.17 41.9
1.39 48.9
1.53 45.4
1.74 48.3
2.01 50.7
2.12 50.6
2.29 49.2
2.52 49
2.61 45.9
2.64 49.8
2.87 49.4
3.39 51.4
3.41 49.7
3.52 49.8
3.65 49.9
;
run;



/* Initial scatterplot of data: */

proc sgplot data=rabbit;
scatter y=length x=age;
run;

/* Not a linear trend? */

* Defining transformed variables for the polynomial terms:;

data rabbit; set rabbit;
age2 = age**2;
age3 = age**3;
age4 = age**4;
run;

proc reg data=rabbit;
model length = age age2 age3 age4; * quartic regression;
run;


** Looks like the fourth-degree term is not needed.;

proc reg data=rabbit;
model length = age age2 age3; * cubic regression;
run;

* Plotting the fitted curve (this plotting approach only works with linear, quadratic, or cubic);

title 'Cubic Regression Plot';
proc sgplot data=rabbit;
reg y=length x=age / degree=3;
run;
title;

/*****************************************/

/* Try with centered predictors ... */

PROC MEANS DATA = rabbit; 
VAR age;
OUTPUT OUT = agestats MEAN(age) = meanage;
run;

DATA rabbit2;
IF _N_=1 THEN SET agestats;
SET rabbit;
cage = age - meanage;
cage2 = cage**2;
cage3 = cage**3;
cage4 = cage**4;
run;

* Defining transformed variables for the polynomial terms:;


proc reg data=rabbit2;
model length = cage cage2 cage3 cage4; * quartic regression;
run;


** Looks like the fourth-degree term is not needed.;

proc reg data=rabbit2;
model length = cage cage2 cage3; * cubic regression;
run;

* Plotting the fitted curve (this plotting approach only works with linear, quadratic, or cubic);

title 'Cubic Regression Plot, centered age';
proc sgplot data=rabbit2;
reg y=length x=cage / degree=3;
run;
title;