"
Exercise 11.4
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 P Lutman, Rothamsted Research.
Version 1, 3/05/2015
Set working directory - use setwd() function or from Session menu in RStudio
e.g. SET [WORKINGDIRECTORY=d:/stats4biol/data']
"
" Read data from working directory "
FILEREAD [NAME='SCAB.DAT'; IMETHOD=read] FGROUPS=no,5(yes),no
" Create new factor Type "
FACTOR [LEVELS=2; LABELS=!T(Control,Sulphur)] IDENTIFIER=Type
CALCULATE Type = (Treatment.ne.1) + 1
" ANOVA (logit scale) with factorial + added control treatment structure "
TREATMENTSTRUCTURE Type/(Sulphur*Timing)
ANOVA [FPROBABILITY=yes; PSE=means,lsd,differences] Y=LOGIT(Scab)
" Save standardized residuals "
AKEEP [RESIDUALS=Stres; RMETHOD=standardized]
" Plot residuals against row number "
TRELLIS [PEN=Col; TITLE='Residuals vs row number'; YPEN=same; KEY=0] Stres; Row; METHOD=point,line
" Plot residuals against column number "
TRELLIS [PEN=Row; TITLE='Residuals vs column number'; YPEN=same; KEY=0] Stres; Col; METHOD=point,line
" Get treatment allocation to rows and columns "
TABULATE [CLASS=Row,Col] !(#Treatment); MEAN=Tmean
PRINT Tmean; DECIMALS=0; FIELD=6
" Add column into model as extraneous factor - not balanced so use GenStat regression "
MODEL LOGIT(Scab)
FIT [PRINT=model,summary,accumulated; FPROB=yes] Col + Type/(Sulphur*Timing)
RCHECK
" Use type deviance to get standardized residuals "
RKEEP [RMETHOD=deviance] RESIDUALS=Stres2
" Plot residuals against row number "
TRELLIS [PEN=Col; TITLE='Residuals vs row number'; YPEN=same; KEY=0] Stres2; Row; METHOD=point,line
" Plot residuals against column number "
TRELLIS [PEN=Row; TITLE='Residuals vs column number'; YPEN=same; KEY=0] Stres2; Col; METHOD=point,line
" Get predictions "
PREDICT [COMBINATIONS=present; PREDICT=Tpred; SE=Tse] Type,Sulphur,Timing
" Remove missing combinations "
TABULATE [CLASS=Type,Sulphur,Timing; COUNT=Tc]
CALC Tpred,Tse = mvinsert(Tpred,Tse; Tc.eq.0)
PRINT Tpred,Tse
" End of file "