/*
Multiple Regression Examples/Nonlinear Regression Examples
FILE: MReg-NLinReg-Examples.doc
Examples (with modifications) from
Piegorsch and Bailer (2005)
ANALYZING ENVIRONMENTAL DATA
Wiley
*/
* ======;
* SAS code to fit quadratic regression to Rice Yield Data;
data rice;
input yield mint @@;
mintbar = 23.1975;
x = mint - mintbar;
x2 = x*x;
datalines;
2.8 27.2 3.1 28.1 2.3 29.2 3.0 23.4
2.1 25.7 2.3 26.2
3.6 26.3 2.4 25.9 3.1 26.0
3.8 26.2 2.4 26.4 4.4 23.4
3.5 24.5 3.0 23.7 2.6 23.1
4.2 21.7 3.1 22.5 4.8 23.1
3.1 23.2 3.3 23.6 4.5 23.7
2.9 23.9 3.8 24.4 3.1 24.0
3.2 23.9 3.7 23.5 3.3 22.9
4.7 20.0 3.5 23.7 3.4 22.5
3.2 23.0 3.2 22.4 4.5 21.2
5.0 19.2 6.2 17.4 6.2 19.0
6.0 19.0 6.1 18.8 7.3 18.0 6.6 18.0
proc reg;
model yield = x x2 ;
plot yield*x p.*x/overlay;
run;
* ======;
* SAS code to fit equal-slopes ANCOVA;
* equal-slopes ANCOVA for Faba Bean Genotoxicity Data
data faba;
input etrt $ time percent @@;
y = log(percent / (100 - percent));
timebar = 13.47826;
covar = time - timebar;
datalines;
FokI 0 20.5 FokI 0 13.4 FokI 0 17.2 FokI 0 19.6
FokI 10 65.4 FokI 10 56.0 FokI 10 57.3 FokI 10 59.2
FokI 20 68.4 FokI 20 54.7 FokI 20 62.8 FokI 20 61.8
FokI 30 78.1 FokI 30 72.8 FokI 30 81.9 FokI 30 79.4
EcoRI 0 6.8 EcoRI 0 13.0 EcoRI 10 10.5 EcoRI 10 12.6
EcoRI 20 49.0 EcoRI 20 40.6 EcoRI 30 62.1 EcoRI 30 50.2
MgCl2 0 2.6 MgCl2 0 2.4 MgCl2 0 9.9 MgCl2 0 14.4
MgCl2 0 18.2 MgCl2 0 21.9 MgCl2 10 24.5 MgCl2 10 28.0
MgCl2 20 56.2 MgCl2 20 56.0 MgCl2 20 52.7 MgCl2 20 42.6
MgCl2 20 28.0 MgCl2 20 31.6 MgCl2 30 48.2 MgCl2 30 47.8
MnCl2 0 58.9 MnCl2 0 65.2
MnCl2 10 91.5 MnCl2 10 89.9 MnCl2 20 96.6 MnCl2 20 96.6
proc gplot;
plot y*covar=etrt;
run;
proc glm;
class etrt;
model Y = covar etrt covar*etrt / solution;
run;
* ======;
* SAS code to fit Michaelis-Menten curve;
data puromycn;
input conc vel @@;
datalines;
0.02 76 0.02 47 0.06 97 0.06 107
0.11 123 0.11 139 0.22 159 0.22 152
0.56 191 0.56 201 1.10 207 1.10 200
proc nlin method=marquardt;
parameters b1 = 195.80271
b2 = 0.04841 ;
model vel = b1*conc/(b2 + conc) ;
der.b1 = conc/(b2 + conc);
der.b2 = -b1*conc/((b2 + conc)**2);
output out=mmout predicted=vel_hat residual=vel_res;
run;
proc gplot;
plot vel*conc vel_hat*conc/overlay;
plot vel_res*vel_hat / vref=0;
run;
* ======;
* SAS code to fit segmented regr. model;
* continuous three-segment model to Pond Water Level Data;
data dpond;
input aquifer pond @@;
datalines;
41.54 54.267 43.60 54.227 42.92 54.207
42.45 54.217 42.61 54.257 43.01 54.347
43.88 54.267 43.31 54.21 41.87 54.207
41.52 54.237 41.29 54.217 44.53 54.257
44.22 54.247 44.03 54.227 43.45 54.217
42.85 54.227 42.25 54.197 41.64 54.177
40.63 54.197 42.21 54.257 42.63 54.247
42.37 54.237 41.92 54.227 41.38 54.227
41.09 54.227 40.80 54.237 40.64 54.257
41.54 54.247 41.61 54.237 43.60 54.267
43.99 54.267 43.78 54.247 43.58 54.267
43.13 54.227 43.17 54.247 42.75 54.247
42.37 54.237 41.95 54.207 41.49 54.187
40.90 54.157 40.38 54.087 39.99 54.007
39.93 53.987 39.74 53.957 39.41 53.827
40.20 53.957 40.79 53.987 39.71 53.897
39.28 53.807 38.51 53.627 38.06 53.487
36.58 53.527 36.58 53.487 36.58 53.487
36.56 53.487 36.71 53.487
36.98 53.487 37.18 53.487 38.15 53.607
/*
E(Y) = b0 [aquifer <= tau1 ]
= b0 + b1*(aquifer-tau1) [tau1 < aquifer <= tau2 ]
= b0 + b1*(tau2-tau1) + b3*(aquifer-tau2)
[aquifer > tau2 ]
*/
procnlindata=dpond method=dud;
parameters b0 = 53.5
b1 = 0.1641
b3 = 0.0114
tau1 = 38.0
tau2 = 41.0 ;
if aquifer <= tau1 thendo; *left segment;
model pond = b0;
end;
elseif tau1 < aquifer <= tau2 thendo; * middle segment;
model pond = b0 + b1*(aquifer-tau1) ;
end;
elsedo; *right segment;
model pond = b0 + b1*(tau2-tau1)+b3*(aquifer-tau2);
end;
outputout=pondout predicted=pond_hat residual=pond_res;
run;
procgplotdata=pondout;
plot pond*aquifer pond_hat*aquifer/overlay;
plot pond_res*pond_hat / vref=0;
run;
* ======;
* SAS code to fit monoexponential model;
data pine;
input diam trees @@;
datalines;
1.0 120.905 3.0 72.846 5.0 41.289 7.0 25.103
9.0 15.020 11.0 10.489 13.0 7.630 15.0 7.090
17.0 4.302 19.0 2.563 21.0 2.838 29.0 0.107
proc nlin method=marquardt;
parameters b0 = 0.107 b1 = 110.767782 b2 = 0.19358046 ;
model trees = b0 + b1*exp(-b2*diam) ;
der.b0 = 1;
der.b1 = exp(-b2*diam);
der.b2 = -diam*b1*exp(-b2*diam);
outputout=pineout predicted=trees_hat residual=trees_res;
run;
procgplotdata=pineout;
plot trees*diam trees_hat*diam /overlay;
plot trees_res*trees_hat / vref=0;
run;
* ======;
* SAS code to fit logistic growth curve;
data mold;
input time length @@;
datalines;
6 0.00 7 3.25 8 12.99 9 58.44
10 128.25 11 224.03 12 394.58 14 456.33
proc nlin data=mold method=marquardt;
parameters b0 = 456.33
b1 = -13.9130705492334
b2 = 1.29775227563122 ;
model length = b0/(1+exp(-b1-b2*time)) ;
der.b0 = 1/(1+exp(-b1-b2*time));
der.b1 = b0*exp(-b1-b2*time)/((1+exp(-b1-b2*time))**2);
der.b2 = time*b0*exp(-b1-b2*time)/((1+exp(-b1-b2*time))**2);
outputout=moldout predicted=length_hat residual=length_res;
run;
procgplotdata=moldout;
plot length*time length_hat*time /overlay;
plot length_res*length_hat / vref=0;
run;