" Exercise 15.1

 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 R. Harrington & C. Shortall, Rothamsted Research
 Version 1, 26/08/2014 "

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

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

" Plot data by site "
TRELLIS [GROUPS=Site; YTITLE='Wet weight'; XTITLE='Year'] Y=WetWeight; X=Year

" Calculate log10-transform of biomass as used in previous questions "
CALCULATE logWt=LOG10(WetWeight+0.5)

" Regression with groups using log-transformed data "
MODEL Y=logWt
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Year*Site

" Get parameter estimates in abbreviated format "
FIT [PRINT=estimates; TPROBABILITY=yes; CONSTANT=omit] Site/Year

" 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

" Investigate apparent variance heterogeneity further "
DGRAPH [WINDOW=5; KEYWINDOW=0; TITLE='Fitted values plot'] Y=Resid; X=Fitted; PEN=Site
DGRAPH [WINDOW=6; KEYWINDOW=0; SCREEN=keep; TITLE='Absolute residuals plot'] \
  Y=ABS(Resid); X=Fitted; PEN=Site

" For comparison: regression with groups with untransformed data - much worse residual plots "
MODEL Y=WetWeight
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Site*Year
RKEEP [RMETHOD=DEVIANCE] RESIDUALS=Resid; FITTEDVALUES=Fitted
DRESIDUALS [RESIDUALS=Resid; FITTEDVALUES=Fitted] \
  METHOD=fittedvalues,normal,histogram,absresidual

" Extra: taking things further "

" Combine across non-Hereford sites? "
FACTOR [LEV=2; LABELS=!T(Hereford,'Not Hereford')] Hereford
calc Hereford = NEWLEVEL(Site; !(1,2,2,2))
MODEL Y=logWt
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Year*(Hereford/Site)

" Final model: separate intercepts, common slope for non-Hereford sites "
FIT [PRINT=#,accumulated; FPROBABILITY=yes; TPROBABILITY=yes] Year*Hereford+Hereford.Site

" End of File "