REGRESSION
1. Regression with confidence Intervals
2. Surface Response Model
3. Cate-Nelson (precursor to Linear Plateau)
4. Linear-Plateau
5. Linear-Linear
6. Difference in Slope and Intercept
2013 Josh Pittman (Linear Plateau)
2013 Romulo
Lollato (Oklahoma Wheat Yields,
YR vs. avg Yield, Linear Plateau)
Cate-Nelson
(ppt file) Cate
Nelson SSSAJ Article
Surface Response
Model, + Scatter
(both in same program)
linear and quadratic relationships of x and y with z and a linear
interaction term.
Z = x x2 y y2 xy
data one;
input year yield ndvi5 ndvi7;
cards;
2017 23.86 0.31 0.32
2017 12.68 0.26 0.31
2017
38.18 0.36 0.41
2017 36.44 0.32 0.40
2017
66.91 0.39 0.63
2017 67.78 0.39 0.65
2017
83.73 0.30 0.61
2017 47.13 0.23 0.39
2017
55.68 0.37 0.51
2017 60.80 0.35 0.53
2017
71.69 0.39 0.64
2017 83.77 0.42 0.66
2017
89.72 0.39 0.66
2017 55.96 0.35 0.51
2017
15.23 0.27 0.26
2017 16.70 0.24 0.28
2017
29.50 0.30 0.34
2017 35.25 0.27 0.45
2017
52.98 0.27 0.48
2017 79.72 0.36 0.69
2017
83.07 0.30 0.63
2017 44.10 0.28 0.43
2017
71.35 0.34 0.60
2017 73.39 0.37 0.65
2017
59.98 0.40 0.59
2017 78.23 0.33 0.60
2017
73.08 0.29 0.54
2017 77.88 0.35 0.54
2017 9.36
0.24 0.25
2017 18.25 0.31 0.29
2017 33.32 0.39
0.42
2017 56.18 0.43 0.58
2017 65.86 0.37 0.63
2017 71.17 0.31 0.60
2017 75.54 0.41 0.67
2017
64.45 0.34 0.57
2017 65.30 0.40 0.55
2017
70.46 0.35 0.60
2017 67.78 0.44 0.68
2017
64.41 0.45 0.62
2017 78.33 0.40 0.70
2017
58.04 0.42 0.59
2017 15.29 0.28 0.29
2017
23.62 0.31 0.30
2017 22.75 0.27 0.33
2017
41.38 0.32 0.42
2017 58.58 0.36 0.58
2017
66.71 0.36 0.59
2017 76.53 0.38 0.68
2017
48.02 0.29 0.45
2017 60.78 0.32 0.57
2017
73.81 0.40 0.61
2017 60.00 0.40 0.55
2017
73.87 0.39 0.62
2017 81.09 0.38 0.68
2017
56.18 0.36 0.57
2015 25.35 0.35 0.36
2015
21.25 0.31 0.30
2015 28.09 0.39 0.42
2015
33.89 0.39 0.51
2015 29.10 0.37 0.56
2015
29.29 0.37 0.59
2015 32.68 0.27 0.40
2015
22.95 0.30 0.35
2015 27.49 0.35 0.53
2015
28.26 0.38 0.56
2015 39.29 0.40 0.61
2015
34.51 0.33 0.53
2015 36.66 0.37 0.63
2015
27.67 0.38 0.59
2015 23.55 0.33 0.33
2015
25.86 0.38 0.40
2015 28.67 0.38 0.42
2015
24.03 0.41 0.54
2015 35.80 0.45 0.67
2015
33.04 0.30 0.45
2015 31.03 0.28 0.41
2015
24.16 0.34 0.41
2015 36.50 0.33 0.53
2015
38.62 0.41 0.64
2015 30.77 0.42 0.58
2015
36.85 0.37 0.56
2015 35.05 0.36 0.57
2015
36.08 0.32 0.53
2015 25.23 0.37 0.40
2015
26.43 0.38 0.39
2015 28.13 0.39 0.47
2015
32.86 0.42 0.74
2015 20.37 0.39 0.64
2015
28.08 0.33 0.53
2015 30.97 0.26 0.37
2015
33.11 0.30 0.42
2015 26.03 0.37 0.60
2015
35.99 0.37 0.62
2015 36.62 0.33 0.56
2015
22.61 0.40 0.58
2015 32.21 0.31 0.48
2015
26.31 0.41 0.62
2015 28.32 0.38 0.40
2015
27.15 0.39 0.41
2015 30.41 0.39 0.45
2015
32.24 0.46 0.59
2015 27.12 0.39 0.58
2015
30.72 0.34 0.54
2015 29.22 0.32 0.51
2015
37.88 0.39 0.54
2015 37.72 0.38 0.58
2015
32.61 0.40 0.60
2015 36.16 0.40 0.59
2015
35.92 0.37 0.57
2015 30.10 0.38 0.61
2015
38.76 0.39 0.62
data one; set one; if year = 2017 then
delete;
sndvi=ndvi5 + ndvi7;
proc corr; var sndvi ndvi5 ndvi7 yield;
run;
proc print;
proc rsreg data = one out
= two;
model yield = ndvi5 ndvi7 /predict;
proc g3grid data = two out = three;
grid
ndvi5*ndvi7=yield/spline;
proc g3d
data = three gout=new;
plot ndvi5*ndvi7=yield;
proc g3d data = one
gout=new2;
scatter
ndvi5*ndvi7=yield/shape='pyramid';
scatter
ndvi5*ndvi7=yield/shape='diamond';
run;
Linear-Plateau Program (Romulo data)
data one;
input x y;
/* Avg. wheat yield data over time */
cards;
1950 604.892513
1951 638.4976526
1952 1243.390166
1953 806.5233506
1954 1008.154188
1955 537.6822338
1956 1108.969607
1957 840.1284902
1958 1747.46726
1959 1344.205584
1960 1747.46726
1961 1613.046701
1962 1276.995305
1963 1411.415864
1964 1545.836422
1965 1881.887818
1966 1411.415864
1967 1142.574747
1968 1545.836422
1969 1881.887818
1970 1747.46726
1971 1344.205584
1972 1545.836422
1973 2016.308377
1974 1411.415864
1975 1613.046701
1976 1613.046701
1977 1814.677539
1978 1814.677539
1979 2553.99061
1980 2016.308377
1981 1814.677539
1982 2217.939214
1983 2352.359773
1984 2419.570052
1985 2016.308377
1986 1949.098097
1987 1814.677539
1988 2419.570052
1989 1814.677539
1990 2217.939214
1991 1814.677539
1992 1949.098097
1993 2016.308377
1994 1814.677539
1995 1411.415864
1996 1276.995305
1997 2150.728935
1998 2553.99061
1999 2352.359773
2000 2285.149493
2001 2217.939214
2002 1881.887818
2003 2621.20089
2004 2352.359773
2005 2150.728935
2006 1613.046701
2007 1881.887818
2008 2486.780331
2009 1478.626143
2010 2083.518656
2011 1478.626143
2012 2419.570052
proc print;
proc nlin data = one best = 3;
parms b0=500 to 1000 by 100 b1=10 to 50 by
10 njoint=1970 to 1990 by 5;
if x<njoint
then do;
model y = b0 + b1*x;
der.b0=1;
der.b1=x;
der.njoint=0;
end;
else do;
model
y=b0+b1*njoint;
der.b0=1;
der.b1=njoint;
der.njoint=b1;
end;
file print;
if _obs_ =1 and _model_ =0 then do;
plateau = b0 + b1*njoint;
put
plateau=;
end;
plateau=b0+b1*njoint;
id plateau;
output out = new p = pry
parms=b0 b1 njoint sse=sse;
run;
proc plot;
plot y*x='+'
pry*x='*'/overlay;
run;
proc means noprint;
var y sse b0 b1 njoint
plateau;
output out = new2 n = tdf
mean = y sse b0 b1 njoint plateau
css=csst;
data new3; set new2;
intercpt=b0; slope=b1; joint=njoint;
rsq=(csst-sse)/csst;
edf=tdf-3;
ssr=csst-sse;
msr=ssr/2;
mse=sse/edf;
f=msr/mse;
probf=1-(probf(f,2,edf));
keep intercpt slope joint plateau
rsq f probf;
proc print;
run;
Linear-Plateau Program (r2 peak, GDD, Yield, BRUNO/Jagman)
data one;
input x y;
cards;
50 0.56
66 0.5
69 0.63
75 0.71
80 0.82
86 0.82
100 0.89
105 0.93
112 0.94
119 0.6
126 0.69
PROC PRINT;
data one; set one;
if
exp = 'e502' then delete;
/* LINEAR PLATEAU */
PROC NLIN DATA = ONE
BEST = 3;
PARMS B0=0.3 to 0.5 by 0.01 B1=0.01 to 0.03 by 0.001 NJOINT=80
to 90 by 1;
IF X<NJOINT THEN DO;
MODEL Y = B0 + B1*X;
DER.B0=1;
DER.B1=X;
DER.NJOINT=0;
END;
ELSE DO;
MODEL Y=B0+B1*NJOINT;
DER.B0=1;
DER.B1=NJOINT;
DER.NJOINT=B1;
END;
FILE PRINT;
IF
_obs_ =1 AND _MODEL_ =0 THEN DO;
PLATEAU = B0 + B1*NJOINT;
PUT
PLATEAU=;
END;
PLATEAU=B0+B1*NJOINT;
ID PLATEAU;
OUTPUT OUT =
NEW P = PRY PARMS=B0 B1 NJOINT SSE=SSE;
RUN;
PROC PLOT;
PLOT
Y*X='+' PRY*X='*'/OVERLAY;
RUN;
PROC MEANS NOPRINT;
VAR Y SSE B0 B1
NJOINT PLATEAU;
OUTPUT OUT = NEW2 N = TDF
MEAN = Y SSE B0 B1 NJOINT
PLATEAU
CSS=CSST;
DATA NEW3; SET NEW2;
INTERCPT=B0; SLOPE=B1;
JOINT=NJOINT;
RSQ=(CSST-SSE)/CSST;
EDF=TDF-3;
SSR=CSST-SSE;
MSR=SSR/2;
MSE=SSE/EDF;
F=MSR/MSE;
PROBF=1-(PROBF(F,2,EDF));
KEEP INTERCPT SLOPE JOINT PLATEAU RSQ F PROBF;
PROC PRINT;
RUN;
data one;
input rep trt x y;
cards;
1950 604.892513
1951 638.4976526
1952 1243.390166
1953 806.5233506
1954 1008.154188
1955 537.6822338
1956 1108.969607
1957 840.1284902
1958 1747.46726
1959 1344.205584
1960 1747.46726
1961 1613.046701
1962 1276.995305
1963 1411.415864
1964 1545.836422
1965 1881.887818
1966 1411.415864
1967 1142.574747
1968 1545.836422
1969 1881.887818
1970 1747.46726
1971 1344.205584
1972 1545.836422
1973 2016.308377
1974 1411.415864
1975 1613.046701
1976 1613.046701
1977 1814.677539
1978 1814.677539
1979 2553.99061
1980 2016.308377
1981 1814.677539
1982 2217.939214
1983 2352.359773
1984 2419.570052
1985 2016.308377
1986 1949.098097
1987 1814.677539
1988 2419.570052
1989 1814.677539
1990 2217.939214
1991 1814.677539
1992 1949.098097
1993 2016.308377
1994 1814.677539
1995 1411.415864
1996 1276.995305
1997 2150.728935
1998 2553.99061
1999 2352.359773
2000 2285.149493
2001 2217.939214
2002 1881.887818
2003 2621.20089
2004 2352.359773
2005 2150.728935
2006 1613.046701
2007 1881.887818
2008 2486.780331
2009 1478.626143
2010 2083.518656
2011 1478.626143
2012 2419.570052
Proc print;
proc nlin data = one best = 2;
parms b0=500 to 1000 by 100 b1=10 to 50 by 5 joint=1970 to 1990 by 5
b2 = 200 to 300 by 10;
if x<joint then do;
model y = b0 + b1*x;
der.b0=1;
der.b1=x;
der.joint=0;
der.b2=0;
end;
else do;
model y=b0+(b1-b2)*joint+b2*x;
der.b0=1;
der.b1=joint;
der.joint=b1-b2;
der.b2=x-joint;
end;
file print;
if _obs_ =1
and _model_ =0 then do;
joinlev = b0
+ b1*joint;
put joinlev=;
end;
joinlev=b0+b1*joint;
id joinlev;
output out =
new p = pry parms=b0 b1 joint b2 sse=sse;
run;
proc plot;
plot y*x='+'
pry*x='*'/overlay;
run;
proc means noprint;
var y sse b0 b1 joint b2
joinlev;
output out = new2 n = tdf
mean = y sse b0 b1 joint b2 joinlev
css=csst;
data new3; set new2;
intercpt=b0; slope=b1; joint=joint; slope2=b2; jresp=joinlev;
rsq=(csst-sse)/csst;
edf=tdf-4;
ssr=csst-sse;
msrg=ssr/3;
mse=sse/edf;
f=msrg/mse;
probf=1-(probf(f,2,edf));
keep intercpt slope joint slope2 joinlev rsq msrg mse edf f probf;
proc print;
run;
Linear Plateau (class project) SOIL 5112 ***
Stability Analysis (treatment by environment interactions) click to go to this page
Regression Analysis with Confidence Intervals (point and line)
data one;
input red nir NDVI yield height;
cards;
.5 .6 .7 4000 32
.4 .7 .8 5000 39
.5 .8 .9 5500 41
.4 .5 .7 3500 30
.3 .5 .6 3400 29
.2 .4 .5 2200 25
.24 .45 .55 2400 28
.29 .46 .58 2600 29
proc reg data = one;
model yield=height/p clm;
output out = dog p = pyd l95m=l95m u95m=u95m;
proc gplot; plot pyd*height yield*height='x' l95m*height
u95m*height/overlay;
run;
cli 95% confidence limits for an individual
predicted value
clm 95% confidence limits for the expected value of the dependent variable
(mean)
(if you use cli, get rid of all the "m's" on the end of the l95 and u95
values)
(can match this up with the Magruder example of NPK vs NP) Magruder excel file, and examples below.
input trt $ envmean yld;
if trt ='np' then intc_dif=0;
if trt ='npk' then intc_dif=1;
slop_dif=intc_dif*envmean;
cards;
np 8.7 6.5
np 28.15 28.4
np 25.36666667 28.6
np 22.25 22.9
np 16.06666667 18
np 24.48333333 26.1
np 19.7 20.2
np 29.01666667 30.3
np 10.58333333 11.7
np 24.23333333 24.4
np 28.31666667 30.6
np 6.8 8.7
np 9.566666667 10.9
np 9.983333333 11.9
np 22.51666667 24.1
np 8 6.1
np 16.85 20.9
np 21.2 22.8
np 29.75 34.4
np 17.35 17.4
np 23.75 26.4
np 21.25 21.4
np 16.55 17.1
np 26.41666667 32
np 14.63333333 12.5
np 5.583333333 5.4
np 16.2 15.1
np 16.05 15.8
np 33.41666667 36.9
np 36.91666667 39.5
np 27.7 34
np 24.1 26.1
np 23.95 28.5
np 34.35 41.5
np 16.58333333 20.7
np 31.83333333 30.7
np 35.65 49.3
np 9.35 10.1
np 19.38333333 23.5
np 21.46666667 25.4
np 24.93333333 23.5
np 30.95 36.2
np 30.3 38.9
np 34.78333333 44.1
np 29.78333333 38.8
np 38.48333333 51.4
np 36.21666667 45.6
np 21.86 32.3
np 26.78333333 32.2
np 44.88333333 52.6
np 35.68333333 43
np 31.26666667 38.3
np 36.75 32.2
np 24.38333333 27.9
np 33.15 31.8
np 21.6 22.2
np 17.23333333 13.2
np 12.38333333 11.7
np 25.76666667 25
np 20.11166667 16.3
np 28.53333333 31.9
np 28.06666667 23.7
np 21.90483333 24.0368
np 30.59855 29.1516
np 22.15076667 31.946
np 5.911566667 9.258
np 21.42166667 22.19
np 46.43333333 60.8
np 28.97 33.81
np 38.34555598 52.06795346
np 31.45561744 33.1039001
np 26.78102444 24.87453073
np 32.81616667 41.665
np 41.4 51.6
np 45.7383689 54.67134146
np 32.16666667 31
np 35.33551989 43.63420513
np 3.909166667 6.068
np 42.8735 45.706
np 4.661666667 2.61
np 31.6285833 35.67635604
np 17.65385104 23.14757086
np 34.99166667 44.37
np 39.20333333 52.08
np 28.09405289 32.9500485
npk 8.7 5.5
npk 28.15 32.3
npk 25.36666667 22.7
npk 22.25 25.1
npk 16.06666667 21.9
npk 24.48333333 27
npk 19.7 20.6
npk 29.01666667 32.2
npk 10.58333333 12.4
npk 24.23333333 26.7
npk 28.31666667 33.6
npk 6.8 8.2
npk 9.566666667 9.9
npk 9.983333333 10.9
npk 22.51666667 23.1
npk 8 9.9
npk 16.85 15.1
npk 21.2 24.1
npk 29.75 34.4
npk 17.35 19.7
npk 23.75 21.4
npk 21.25 24.2
npk 16.55 16.7
npk 26.41666667 32.1
npk 14.63333333 15.3
npk 5.583333333 2.5
npk 16.2 15.6
npk 16.05 17
npk 33.41666667 35.7
npk 36.91666667 39.4
npk 27.7 35.2
npk 24.1 27.6
npk 23.95 27
npk 34.35 32.3
npk 16.58333333 22.2
npk 31.83333333 29.9
npk 35.65 34.5
npk 9.35 9.9
npk 19.38333333 23.8
npk 21.46666667 27.1
npk 24.93333333 31
npk 30.95 29.6
npk 30.3 37.1
npk 34.78333333 43.3
npk 29.78333333 30.4
npk 38.48333333 47.8
npk 36.21666667 45.3
npk 21.86 23.8
npk 26.78333333 33.7
npk 44.88333333 50.3
npk 35.68333333 37
npk 31.26666667 32.6
npk 36.75 40.3
npk 24.38333333 25.4
npk 33.15 32.6
npk 21.6 23.4
npk 17.23333333 21.3
npk 12.38333333 12.3
npk 25.76666667 29.7
npk 20.11166667 25.07
npk 28.53333333 32.2
npk 28.06666667 42.1
npk 21.90483333 31.5786
npk 30.59855 36.9428
npk 22.15076667 31.6512
npk 5.911566667 8.5061
npk 21.42166667 24.05
npk 46.43333333 62.6
npk 28.97 37.34
npk 38.34555598 52.84454497
npk 31.45561744 38.11772693
npk 26.78102444 28.03528683
npk 32.81616667 40.716
npk 41.4 59.5
npk 45.7383689 55.55670732
npk 32.16666667 38
npk 35.33551989 45.08231665
npk 3.909166667 4.789
npk 42.8735 45.857
npk 4.661666667 8.27
npk 31.6285833 37.01767818
npk 17.65385104 23.81803027
npk 34.99166667 46.67
npk 39.20333333 49.43
npk 28.09405289 31.94973331
data two; set one;
proc sort; by trt;
proc reg;
model yld = envmean intc_dif slop_dif;
run;
proc reg;
by trt;
model yld=envmean;
run;
data one;
input gdd ndvi yld;
if gdd = 87 then intc_dif=0;
if gdd = 106 then intc_dif=1;
slop_dif=intc_dif*ndvi;
cards;
87 0.363 2142.42
87 0.217 772.01
87 0.439 3285.61
87 0.448
2850.91
87 0.515 4684.25
87 0.510 4641.91
87 0.543 5315.12
87
0.368 2646.27
87 0.481 4358.23
87 0.492 4259.44
87 0.500 5277.02
87 0.515 5247.38
87 0.562 5384.28
87 0.499 4307.42
87 0.334 1589.17
87 0.378 1563.77
87 0.416 2745.06
87 0.527 2431.75
87 0.489 3240.45
87 0.417 5521.18
87 0.471 5483.07
87 0.393 3322.30
87 0.457 5233.26
87 0.495 5202.21
87 0.471 4966.52
87 0.499 5260.08
87 0.604 5629.85
87 0.434 5569.16
87 0.380 2273.68
87 0.444 2721.07
87 0.469 3908.01
87 0.486 4409.04
87 0.551 4811.27
87 0.533 5039.91
87 0.547 5387.10
87 0.416 4486.66
87 0.531 4980.63
87 0.507 4852.20
87 0.530 5346.17
87 0.550 5000.39
87 0.574 5325.00
87 0.529 4898.78
87 0.386 2248.27
87 0.421 2695.67
87 0.427 2852.33
87 0.515 3549.53
87 0.494 4948.17
87 0.562 5272.78
87 0.590 5267.14
87 0.433 4126.77
87 0.458 4698.36
87 0.536 5175.40
87 0.508 4222.74
87 0.489 5381.45
87 0.640 5377.22
87 0.483 5007.45
106 0.390 2142.42
106 0.232 772.01
106 0.509
3285.61
106 0.505 2850.91
106 0.653 4684.25
106 0.697 4641.91
106 0.731 5315.12
106 0.439 2646.27
106 0.613 4358.23
106 0.641
4259.44
106 0.681 5277.02
106 0.721 5247.38
106 0.772 5384.28
106 0.651 4307.42
106 0.341 1589.17
106 0.449 1563.77
106 0.396
2745.06
106 0.523 2431.75
106 0.665 3240.45
106 0.712 5521.18
106 0.779 5483.07
106 0.505 3322.30
106 0.712 5233.26
106 0.703
5202.21
106 0.653 4966.52
106 0.712 5260.08
106 0.777 5629.85
106 0.662 5569.16
106 0.417 2273.68
106 0.447 2721.07
106 0.531
3908.01
106 0.647 4409.04
106 0.729 4811.27
106 0.717 5039.91
106 0.752 5387.10
106 0.557 4486.66
106 0.700 4980.63
106 0.751
4852.20
106 0.751 5346.17
106 0.705 5000.39
106 0.714 5325.00
106 0.652 4898.78
106 0.364 2248.27
106 0.389 2695.67
106 0.407
2852.33
106 0.573 3549.53
106 0.665 4948.17
106 0.765 5272.78
106 0.785 5267.14
106 0.574 4126.77
106 0.648 4698.36
106 0.717
5175.40
106 0.689 4222.74
106 0.723 5381.45
106 0.811 5377.22
106 0.713 5007.45
data two; set one;
proc sort; by gdd;
proc reg;
model yld = ndvi intc_dif slop_dif;
run;
proc reg;
by gdd;
model yld=ndvi;
run;
EXTRA
Keeping a SUBSET of the ENTIRE DATA SET (keeps data for only the varieties entered)
data
awnless;set gpyt;
if cv = 'DELIVER' or cv = 'WM135' or cv = 'TAM401' or cv = 'LONGHORN' ;
run
;proc print; run;