/* meatph.sas */
options ls=72 ps=60;
data both;
   infile 'meatph.dat';
   input carcass treatment $ time ph;
   logtime=log(time);
   if time=24 then delete;
   conv = (treatment='conven');
   eshb = (treatment='eshb');
   ctime = conv*logtime;
   etime = eshb*logtime;
   ctime2 = ctime*ctime;
   etime2 = etime*etime;   
proc plot data=both;
   plot ph*time=treatment;
   plot ph*logtime=treatment;
title 'meat pH plots';   
proc print data=both;
   var ph logtime conv ctime ctime2 eshb etime etime2;
title 'meat pH variables';
/* fit two quadratics in log(time) */
proc glm data=both;
   model ph=conv ctime ctime2 eshb etime etime2/ noint ss1 ss2 clparm;
      output out=new1 p=yhat1 r=resid1;
title '2 separate quadratics';      
proc plot data=new1;
   plot resid1*time=treatment / vref=0;
proc univariate data=new1 plot normal;
   var resid1;
proc rank data=new1 normal=blom out=new1a;
   var resid1;
   ranks nscore;
proc plot data=new1a;
   plot nscore*resid1;
/* fit linear in log(time) for conv; quadratic for eshb */
proc glm data=both;
   model ph=conv ctime eshb etime etime2/ noint ss1 ss2 clparm;
      output out=new2 p=yhat2 r=resid2 lcl=clpl
             ucl=clpu lclm=clml uclm=clmu;
title 'linear for conv; quadratic for eshb';      
proc plot data=new2;
   plot resid2*time=treatment / vref=0;
proc univariate data=new2 plot normal;
   var resid2;
proc rank data=new2 normal=blom out=new2a;
   var resid2;
   ranks nscore;
proc plot data=new2a;
   plot nscore*resid2;
proc print data=new2;
   var treatment time yhat2 clml clmu clpl clpu;
/* estimate when pH reaches 6 */
data addon;
   input carcass treatment $ time ph;
   logtime=log(time);
   conv = (treatment='conven');
   eshb = (treatment='eshb');
   ctime = conv*logtime;
   etime = eshb*logtime;
   ctime2 = ctime*ctime;
   etime2 = etime*etime;   
cards;
. conven 2.9 .
. conven 3 .
. conven 3.1 .
. conven 3.2 .
. conven 3.3 .
. conven 3.4 .
. conven 3.5 .
. conven 3.6 .
. conven 3.7 .
. conven 3.8 .
. conven 3.9 .
. conven 4.1 .
. conven 4.2 .
. conven 4.3 .
. conven 4.4 .
. conven 4.5 .
. conven 4.6 .
. conven 4.7 .
. conven 4.8 .
. conven 4.9 .
. conven 5   .
. conven 5.1 .
. eshb 1.9 .
. eshb 2 .
. eshb 2.1 .
. eshb 2.2 .
. eshb 2.3 .
. eshb 2.4 .
. eshb 2.5 .
. eshb 2.6 .
. eshb 2.7 .
. eshb 2.8 .
. eshb 2.9 .
. eshb 3   .
. eshb 3.1 .
. eshb 3.2 .
. eshb 3.3 .
;
data when6;
   set both addon;
proc glm data=when6;
   model ph=conv ctime eshb etime etime2/ noint ss1 ss2 clparm;
      output out=new3 p=yhat2 r=resid2 lcl=clpl
             ucl=clpu lclm=clml uclm=clmu;
title 'linear for conv, quadratic for eshb, find when pH is 6';
proc sort data=new3 out=new4;
  by treatment time;
proc print data=new4;
   var treatment ph time yhat2 clml clmu clpl clpu;
run;