" Exercise 18.3

 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
 Version 1, 21/08/2015 "

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

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

" Plot proportion and logit(proportion) against wetness period in adjacent windows "
YAXIS 5; TITLE='Proportion infected'
XAXIS 5; TITLE='Wetness period (h)'
DGRAPH [WINDOW=5; KEY=0] Y=NInf/16; X=Wetness; PEN=Shelf
YAXIS 6; TITLE='Logit(proportion infected)'
XAXIS 6; TITLE='Wetness period (h)'
DGRAPH [WINDOW=6; KEY=0; SCREEN=keep] Y=log((NInf+1)/(17-NInf)); X=Wetness; PEN=Shelf

" Create factor version of Wetness variate "
GROUPS VECTOR=Wetness; FACTOR=fWetness

" Fit logistic GLM to treatment means to check for over-dispersion "
MODEL [DISTRIBUTION=binomial; LINK=logit] Y=NInf; NBINOMIAL=16
FIT [PRINT=mode,summary; FPROBABILITY=yes] Shelf+fWetness

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

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

" Sequential ANOVA table "
FITINDIVIDUALLY [PRINT=accumulated; FPROBABILITY=yes] Shelf+fWetness

" Plot fitted model on linear predictor scale "
RGRAPH [BACKTRANSFORM=none] GROUP=fWetness

" Model fitting linear relationship plus check for lack of fit "
MODEL [DISTRIBUTION=binomial; LINK=logit] Y=NInf; NBINOMIAL=16
FITINDIVIDUALLY [PRINT=accumulated; FPROBABILITY=y; TPROBABILITY=y] Shelf+Wetness+fWetness

" Model fitting linear relationship "
MODEL [DISTRIBUTION=binomial; LINK=logit] Y=NInf; NBINOMIAL=16
FITINDIVIDUALLY [PRINT=#,accumulated; FPROBABILITY=y; TPROBABILITY=y] Shelf+Wetness

" Recheck residuals "
RKEEP RESIDUALS=res; FITTEDVALUES=fval
DRESIDUALS [RESIDUALS=res; FITTEDVALUES=fval] METHOD=fittedvalues,normal,histogram,absresidual

" Make prediction for 36 hr of wetness on logit scale "
PREDICT [PREDICTIONS=Pred; SE=SEpred; BACKTRANSFORM=no] CLASSIFY=Wetness; LEVELS=36

" Calculate 95% CIs for prediction on logit scale "
RKEEP DF=ResDF
CALCULATE Upper=Pred+EDT(0.975;ResDF)*SEpred
CALCULATE Lower=Pred-EDT(0.975;ResDF)*SEpred
PRINT STRUCTURE=Lower,Pred,Upper

" Back-transform to proportion scale "
CALCULATE btPred=ILOGIT(Pred)/100
CALCULATE btUpper=ILOGIT(Upper)/100
CALCULATE btLower=ILOGIT(Lower)/100
PRINT STRUCTURE=btLower,btPred,btUpper

" Plot of fitted model on logit and back-transformed scale "
RGRAPH [BACKTRANSFORM=no]
RGRAPH [BACKTRANSFORM=link]

" End of File "