" Exercise 17.2

 Statistical Methods in Biology: Design and Analysis of Experiments and Regression
 by S.J. Welham, S.A. Gezan, S.J. Clark & A. Mead (2014)
 Chapman & Hall/CRC Press, Boca Raton, Florida. ISBN: 978-1-4398-0878-8
 Data from V. Buchanan-Wollaston (PRESTA), University of Warwick
 Version 1, 11/10/2015 "

" Set working directory - change to location of your data file "
\SET [WORKINGDIRECTORY='d:/stats4biol/data/']

" Read data from working directory "
FILEREAD [NAME='SENESCENCE.DAT'; IMETHOD=read] FGROUPS=2(no),yes,3(no)

" Plot data "
DGRAPH Y=CATMA2A31585; X=Day
DGRAPH Y=CATMA1A09000; X=Day

" Exercise 17.2a "

" Calculate powers of centered day variate "
CALCULATE Day2,Day3,Day4=(Day-MEAN(Day))**2,3,4

" Calculate factor version of Day "
GROUPS VECTOR=Day; FACTOR=fDay

" Polynomial regression on CATMA2A31585 " 

" Set up pens for graphics "
PEN NUMBER=11,12,13,14; METHOD=monotonic; SYMBOLS=0; COLOUR='black','red','green','blue'
PEN NUMBER=15; METHOD=point; SYMBOLS=2; COLOUR='purple'; CFILL='purple'

" Fit quadratic polynomial with observation 32 removed "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit2Aquad
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit2Aquad] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='Quadratic model for CATMA2A31585'] Y=CATMA2A31585,fit2Aquad; X=Day; PEN=1,11
" Check for lack of fit "
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+fDay

" Fit cubic with observation 32 removed "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+Day3
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit2Acub
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit2Acub] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='Cubic model for CATMA2A31585'] Y=CATMA2A31585,fit2Acub; X=Day; PEN=1,11
" Check for lack of fit "
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+Day3+fDay

" Fit quartic with observation 32 removed - just to check "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+Day3+Day4
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=Fitted
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=Fitted] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='Quartic model for CATMA2A31585'] Y=CATMA2A31585,Fitted; X=Day; PEN=1,11


" Polynomial regression on CATMA1A09000 " 

" Fit SLR "
MODEL Y=CATMA1A09000
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=Fitted
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=Fitted] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='SLR for CATMA1A09000'] Y=CATMA1A09000,Fitted; X=Day; PEN=1,11
" Check for lack of fit "
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+fDay

" Fit quadratic "
MODEL Y=CATMA1A09000
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=Fitted
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=Fitted] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='Quadratic for CATMA1A09000'] Y=CATMA1A09000,Fitted; X=Day; PEN=1,11
" Check for lack of fit "
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+fDay

" Fit cubic "
MODEL Y=CATMA1A09000
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+Day3
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit1Acub
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit1Acub] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='Cubic for CATMA1A09000'] Y=CATMA1A09000,fit1Acub; X=Day; PEN=1,11
" Check for lack of fit "
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+Day3+fDay

" Fit quartic - just to check "
MODEL Y=CATMA1A09000
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Day+Day2+Day3+Day4
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=Fitted
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=Fitted] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
DGRAPH [TITLE='Quartic for CATMA1A09000'] Y=CATMA1A09000,Fitted; X=Day; PEN=1,11


" Exercise 17.2b "

" Non-linear model for CATMA2A31585 - again exclude observation 32 " 

" Fit single factor (saturated) model (i.e. day means) to get pure error "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FIT [FPROBABILITY=yes; TPROBABILITY=yes] fDay
" Save pure error sum of squares and df "
RKEEP DEVIANCE=PESS; DF=PEDF; FITTED=fit2ALOF
" Calculate and print pure error mean square "
CALCULATE PEMS=PESS/PEDF
PRINT STRUCTURE=PESS,PEDF,PEMS; DECIMALS=*,0,*

" Fit exponential model "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FITCURVE [CURVE=exponential; FPROBABILITY=yes] Day
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit2Aexp
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit2Aexp] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
RGRAPH
" Save residual sum of squares and df "
RKEEP DEVIANCE=ResSS; DF=ResDF
" Calculate lack-of-fit SS and df "
CALCULATE LOFSS=ResSS-PESS
CALCULATE LOFDF=ResDF-PEDF
PRINT STRUCTURE=LOFSS,LOFDF
CALCULATE F=(LOFSS/LOFDF)/PEMS
CALCULATE prF=CUF(F;LOFDF;PEDF)
PRINT STRUCTURE=F,prF,LOFDF,PEDF; DECIMALS=2(*,0)

" Fit inverse linear model (ldl = linear-divided-by-linear) "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FITCURVE [CURVE=ldl; FPROBABILITY=yes] Day
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit2Aldl
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit2Aldl] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
RGRAPH
" Save residual sum of squares and df "
RKEEP DEVIANCE=ResSS; DF=ResDF
" Calculate lack-of-fit SS and df "
CALCULATE LOFSS=ResSS-PESS
CALCULATE LOFDF=ResDF-PEDF
PRINT STRUCTURE=LOFSS,LOFDF
CALCULATE F=(LOFSS/LOFDF)/PEMS
CALCULATE prF=CUF(F;LOFDF;PEDF)
PRINT STRUCTURE=F,prF,LOFDF,PEDF; DECIMALS=2(*,0)

" Fit quadratic-divided-by-linear (qdl) model "
MODEL Y=MVINSERT(CATMA2A31585;ID.EQ.32)
FITCURVE [CURVE=qdl; FPROBABILITY=yes] Day
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit2Aqdl
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit2Aqdl] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
RGRAPH
" Save residual sum of squares and df "
RKEEP DEVIANCE=ResSS; DF=ResDF
" Calculate lack-of-fit SS and df "
CALCULATE LOFSS=ResSS-PESS
CALCULATE LOFDF=ResDF-PEDF
PRINT STRUCTURE=LOFSS,LOFDF
CALCULATE F=(LOFSS/LOFDF)/PEMS
CALCULATE prF=CUF(F;LOFDF;PEDF)
PRINT STRUCTURE=F,prF,LOFDF,PEDF; DECIMALS=2(*,0)

" Plot the fitted exponential, inverse linear (ldl) and qdl curves with the cubic model and the data "
DGRAPH Y=CATMA2A31585,fit2Acub,fit2Aexp,fit2Aldl,fit2Aqdl,fit2ALOF; X=Day; PEN=1,11,12,13,14,15

" Non-linear model for CATMA1A09000 " 

" Fit single factor (saturated) model (i.e. day means) - get pure error "
MODEL Y=CATMA1A09000
FIT [FPROBABILITY=yes; TPROBABILITY=yes] fDay
" Save pure error sum of squares and df "
RKEEP DEVIANCE=PESS; DF=PEDF; FITTED=fit1ALOF
" Calculate and print pure error mean square "
CALCULATE PEMS=PESS/PEDF
PRINT STRUCTURE=PESS,PEDF,PEMS; DECIMALS=*,0,*

" Fit logistic model "
MODEL Y=CATMA1A09000
FITCURVE [CURVE=logistic; FPROBABILITY=yes] Day
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit1Alog
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit1Alog] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
RGRAPH
" Save residual sum of squares and df "
RKEEP DEVIANCE=ResSS; DF=ResDF
" Calculate lack-of-fit SS and df "
CALCULATE LOFSS=ResSS-PESS
CALCULATE LOFDF=ResDF-PEDF
PRINT STRUCTURE=LOFSS,LOFDF
CALCULATE F=(LOFSS/LOFDF)/PEMS
CALCULATE prF=CUF(F;LOFDF;PEDF)
PRINT STRUCTURE=F,prF,LOFDF,PEDF; DECIMALS=2(*,0)

" Fit Gompertz model "
MODEL Y=CATMA1A09000
FITCURVE [CURVE=gompertz; FPROBABILITY=yes] Day
" Save standardized residuals and fitted values and produce residual plots "
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=fit1Agom
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=fit1Agom] METHOD=fittedvalues,normal,histogram,absresidual
" Plot fitted model "
RGRAPH
" Save residual sum of squares and df "
RKEEP DEVIANCE=ResSS; DF=ResDF
" Calculate lack-of-fit SS and df "
CALCULATE LOFSS=ResSS-PESS
CALCULATE LOFDF=ResDF-PEDF
PRINT STRUCTURE=LOFSS,LOFDF
CALCULATE F=(LOFSS/LOFDF)/PEMS
CALCULATE prF=CUF(F;LOFDF;PEDF)
PRINT STRUCTURE=F,prF,LOFDF,PEDF; DECIMALS=2(*,0)

" Plot the fitted exponential and Gompertz curves with the cubic model and the data "
DGRAPH Y=CATMA1A09000,fit1Acub,fit1Alog,fit1Agom,fit1ALOF; X=Day; PEN=1,11,12,13,15

" End of File "