" Exercise 18.11

 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 D. Gray, Horticulture Research International
 Version 1, 31/08/2015 "

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

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

" Plot data "
DGRAPH Count; Days; PEN=Condition

"  Fit logistic model assuming no over-dispersion "
MODEL [DISTRIBUTION=binomial; LINK=logit] Y=Count; NBINOMIAL=100
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Days*Condition

" Save residual deviance and df "
RKEEP DEVIANCE=ResDev; DF=ResDF
" Test for over-dispersion "
PRINT STRUCTURE=CUCHI(ResDev;ResDF); DECIMALS=4

"  Re-Fit logistic model with dispersion estimated "
MODEL [DISTRIBUTION=binomial; LINK=logit; DISPERSION=*] Y=Count; NBINOMIAL=100
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Days*Condition

" Residual plots based on standardized deviance residuals "
RKEEP RESIDUALS=res; FITTEDVALUES=fval
DRESIDUALS [RESIDUALS=res; FITTEDVALUES=fval] METHOD=fittedvalues,normal,histogram,absresidual

" Re-plot residuals vs Days grouped by Condition "
TRELLIS [GROUP=Condition] res; Days

" Examine fitted model with data "
RGRAPH

" Add quadratic term in Days to model "
CALC DaySqrd = Days*Days
MODEL [DISTRIBUTION=binomial; LINK=logit; DISPERSION=*] Y=Count; NBINOMIAL=100
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Condition*(Days+DaySqrd)

" Still over-dispersion present? Yes "
RKEEP DEVIANCE=ResDev; DF=ResDF
PRINT STRUCTURE=CUCHI(ResDev;ResDF); DECIMALS=4

" Residual plots based on standardized deviance residuals "
RKEEP RESIDUALS=res; FITTEDVALUES=fval
DRESIDUALS [RESIDUALS=res; FITTEDVALUES=fval] METHOD=fittedvalues,normal,histogram,absresidual

" Re-plot residuals vs Days grouped by Condition "
TRELLIS [GROUP=Condition] res; Days

" Plot fitted with data "
TRELLIS [GROUP=Condition; KEY=0] Count,fval; Days; METHOD=point,line

TRELLIS [PEN=Condition; YPEN=same; KEY=0] Count,fval; Days; METHOD=point,line

" Check for lack of fit "
GROUP Days; fDays
MODEL [DISTRIBUTION=binomial; LINK=logit; DISPERSION=*] Y=Count; NBINOMIAL=100
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Condition*(Days+DaySqrd+fDays)

" Compare fit with quadratic model "
RKEEP FITTEDVALUES=fvalfac
TRELLIS [GROUP=Condition; KEY=0] fval,fvalfac; Days; METHOD=line,point

" End of File "