title 'Vitamin B2 in turnip greens';
/*
Data from: Draper & Smith (1966), p. 229
See also: Joiner (1981) 'Lurking variables: Some examples';
*/
data turnip;
	input x1 x2 x3 y;	
	obs = _n_;
	label x1='Radiation'
		x2='Soil moisture'
		x3='Air temperature'
		y='Vitamin B2'
		obs = 'Observation order';
datalines;
    1.76          0.070        7.8        110.4
    1.55          0.070        8.9        102.8
    2.73          0.070        8.9        101.0
    2.73          0.070        7.2        108.4
    2.56          0.070        8.4        100.7
    2.80          0.070        8.7        100.3
    2.80          0.070        7.4        102.0
    1.84          0.070        8.7         93.7
    2.16          0.070        8.8         98.9
    1.98          0.020        7.6         96.6
    0.59          0.020        6.5         99.4
    0.80          0.020        6.7         96.2
    0.80          0.020        6.2         99.0
    1.05          0.020        7.0         88.4
    1.80          0.020        7.3         75.3
    1.80          0.020        6.5         92.0
    1.77          0.020        7.6         82.4
    2.30          0.020        8.2         77.1
    2.03          0.474        7.6         74.0
    1.91          0.474        8.3         65.7
    1.91          0.474        8.2         56.8
    1.91          0.474        6.9         62.1
    0.76          0.474        7.4         61.0
    2.13          0.474        7.6         53.2
    2.13          0.474        6.9         59.4
    1.51          0.474        7.5         58.7
    2.05          0.474        7.6         58.0
;
*-- Check for possible two-variable interactions, and quadratic terms;
proc rsreg data=turnip;
	model y = x1 x2 x3 / noopt;

proc rsreg data=turnip;
	model y = x1 x3 x2 / covar=2 noopt;

*-- Quadratic X2 term, with proc reg;
data turnip;
	set turnip;
	x22 = (x2-.188)**2;   *-- Avoid colinearity;
	label 
		x22 = 'X2**2';
	
proc reg data=turnip;
	model y = x1 x2 x3 x22;
	output out=results r=residual p=yhat;

*include goptions;
proc gplot data=results;
	plot residual * yhat = 1
		/ vaxis=axis1 vref=0;
	plot y * obs = 2
	     y * obs = 3
		/ vaxis=axis1 vm=1 overlay;
	
	axis1 label=(a=90 r=0) length=5in;
	symbol1 v=circle h=2 c=black;
	symbol2 v=circle h=2 c=black i=join;
	symbol3 v=none ci=red i=rl w=3;

/*
proc g3d data=results;
	plot x1 *x2 = yhat;
	label yhat='Pred';
	run;
*/