/* corn.sas */
options ls=72 ps=60;
data corn;
   infile 'corn.dat';
   input year yield;
   newyear=year - 57.5;
   early = (year < 58);
   late = (year > 57);
   earlyyr = early*newyear;
   lateyr = late*newyear;
proc plot data=corn;
   plot yield*year/href=57.5;
title 'plot of Iowa corn yield data';
proc print data=corn;
   var  yield newyear early late earlyyr lateyr;
title 'Iowa corn yield variables';
proc glm data=corn;
   model yield = earlyyr lateyr / p ss1 ss2 clparm;
      output out=new residual=resid predicted=yhat;
      estimate 'Eslope-Lslope' earlyyr 1 lateyr -1;
title 'piecewise linear model: joinpoint 57.5';
proc print data=new;
   var year newyear yield yhat resid;
proc plot data=new;
   plot resid*year / vref=0 href=57.5;
proc univariate data=new plot normal;
   var resid;
title 'descriptive statistics for residuals';   
proc rank data=new normal=blom out=new1;
   var resid;
   ranks nscore;
proc plot data=new1;
   plot nscore*resid;
title 'normal probability plot for residuals';
data without;
   set corn;
   if year = 47 then delete;
proc glm data=without;
   model yield = earlyyr lateyr / p ss1 ss2 clparm;
      output out=new2 residual=resid predicted=yhat;
      estimate 'Eslope-Lslope' earlyyr 1 lateyr -1;
title 'piecewise linear model: joinpoint 57.5 omitting 1947';
proc print data=new2;
   var year newyear yield yhat resid;
proc plot data=new2;
   plot resid*year / vref=0 href=57.5;
run;