" Exercise 18.8

 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 J. Clarkson, University of Warwick
 Version 1, 22/08/2015 "

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

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

" 1. Emergence "

" Plot the data "
TRELLIS [PENGROUP=Rep; GROUPS=Fungicide,BCA; YPEN=same; KEY=0] Y=Emerged; X=Variety; METHOD=point

" Fit full log-linear model with no over-dispersion "
MODEL [DISTRIBUTION=poisson] Y=Emerged
FITINDIVIDUALLY [PRINT=model,summary; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety*Fungicide*BCA

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

" Re-fit full log-linear model  & allow over-dispersion "
MODEL [DISTRIBUTION=poisson; DISPERSION=*] Y=Emerged
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety*Fungicide*BCA

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

" Backwards selection - without recalculating residual deviance "
" Drop 3-way interaction "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide.BCA
" Test 2-way interactions "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
" Drop Fungicide.BCA "
DROP [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
" Re-test other 2-way interactions "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
" Drop Variety.BCA "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
" Test Variety.Fungicide and BCA "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] BCA
" Drop Variety.Fungicide "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
" Test main effects "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] BCA
" Cannot drop any main effects "

" Final predictive model "
MODEL [DISTRIBUTION=poisson; DISPERSION=*] Y=Emerged
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety+Fungicide+BCA

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

" Predict on linear predictor and natural scales "
PREDICT [BACKTRANSFORM=none] CLASSIFY=Variety
PREDICT [PRINT=predictions] CLASSIFY=Variety
PREDICT [BACKTRANSFORM=none] CLASSIFY=Fungicide
PREDICT [PRINT=predictions] CLASSIFY=Fungicide
PREDICT [BACKTRANSFORM=none] CLASSIFY=BCA
PREDICT [PRINT=predictions] CLASSIFY=BCA

" 2. Disease incidence "

" Plot the data "
TRELLIS [PENGROUP=Rep; GROUPS=Fungicide,BCA; YPEN=same] Y=Disease/Emerged; X=Variety; METHOD=point

" Fit logistic regression model without over-dispersion "
MODEL [DISTRIBUTION=binomial; LINK=logit] Y=Disease; NBINOMIAL=Emerged
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety*Fungicide*BCA

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

" Re-fit full logistic regression model allowing over-dispersion "
MODEL [DISTRIBUTION=binomial; LINK=logit; DISPERSION=*] Y=Disease; NBINOMIAL=Emerged
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety*Fungicide*BCA

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

" Backwards selection "
" Drop 3-way interaction "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide.BCA
" Test 2-way interactions "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
" Drop Variety.BCA "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
" Re-test other 2-way interactions "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
" Drop Fungicide.Variety "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
" Test Fungicide.BCA and Variety main effect "
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety
" Cannot drop any more terms "

" Final predictive model "
MODEL [DISTRIBUTION=binomial; LINK=logit; DISPERSION=*] Y=Disease; NBINOMIAL=Emerged
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety+Fungicide*BCA

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

" Predict on linear predictor and natural scales "
PREDICT [PRINT=#,lsd; BACKTRANSFORM=none] CLASSIFY=Fungicide,BCA
PREDICT [PRINT=predictions] CLASSIFY=Fungicide,BCA
PREDICT [BACKTRANSFORM=none] CLASSIFY=Variety
PREDICT [PRINT=predictions] CLASSIFY=Variety


" 3. Analyse number healthy plants to assess best treatment "

" Calculate number of healthy plants "
CALCULATE [PRINT=summary] Healthy=Emerged-Disease

" Fit log-linear model with no over-dispersion "
MODEL [DISTRIBUTION=poisson] Y=Healthy
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety*Fungicide*BCA

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

" Re-fit full log-linear model (dispersion estimated) "
MODEL [DISTRIBUTION=poisson; DISPERSION=*] Y=Healthy
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety*Fungicide*BCA

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

" Backwards selection "
" Drop 3-way interaction "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide.BCA
" Test 2-way interactions "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
" Drop Variety.BCA "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.BCA
" Re-test other 2-way interactions "
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
" Drop Variety.Funngicide "
DROP [PRINT=accumulated; FPROBABILITY=yes] Variety.Fungicide
" Try Fungicide.BCA and Variety "
TRY [PRINT=accumulated; FPROBABILITY=yes] Fungicide.BCA
TRY [PRINT=accumulated; FPROBABILITY=yes] Variety
" Cannot drop any further terms "

" Final predictive model "
MODEL [DISTRIBUTION=poisson; DISPERSION=*] Y=Healthy
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Rep+Variety+Fungicide*BCA

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

" Predict on linear predictor and natural scales "
PREDICT [PRINT=#,lsd; BACKTRANSFORM=none] CLASSIFY=Fungicide,BCA
PREDICT [PRINT=predictions] CLASSIFY=Fungicide,BCA
PREDICT [BACKTRANSFORM=none] CLASSIFY=Variety
PREDICT [PRINT=predictions] CLASSIFY=Variety

" End of File "