" Exercise 6.6

 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. Webster, Rothamsted Research & previously Ecole Polytechnique Fédérale de Lausanne
 Version 1, 20/08/2014 "

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

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

" Check allocation of samples to land use categories "
TABULATE [PRINT=count; CLASS=LandUse]


" Analysis of Cadmium "

" One-way ANOVA of variate Cd with residual plots"
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=Cd
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] Cd

" Formal test for homogeneity of variances "
VHOMOGENEITY [PRINT=test,variance; GROUPS=LandUse] Cd


" One-way ANOVA of log10-transformed Cd with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=LOG10(Cd)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse


" One-way ANOVA of square-root-transformed Cd with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=SQRT(Cd)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] sqrt(Cd)

" Formal test for homogeneity of variances "
VHOMOGENEITY [PRINT=test,variance; GROUPS=LandUse] sqrt(Cd)


" One-way ANOVA of cube-root-transformed Cd with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=(Cd)**(1/3)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] (Cd)**(1/3)

" Formal test for homogeneity of variances "
VHOMOGENEITY [PRINT=test,variance; GROUPS=LandUse] (Cd)**(1/3)

" Save residuals from this model for use later "
AKEEP [RESIDUALS=ResidcubeCd; RMETHOD=standardized]
" Save means and back-transform "
AKEEP TERMS=LandUse; MEANS=cubeCdMeans; SEDMEANS=SED
PRINT cubeCdMeans,SED; DECIMALS=3,4
PRINT STRUCTURE=cubeCdMeans**3


" Analysis of Chromium "

" One-way ANOVA of variate Cr with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=Cr
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] Cr

" Formal test for homogeneity of variances "
VHOMOGENEITY [PRINT=test,variance; GROUPS=LandUse] Cr

" Save residuals from this model for use later "
AKEEP [RESIDUALS=ResidCr; RMETHOD=standardized]
" Save and print means "
AKEEP TERMS=LandUse; MEANS=CrMeans; SEDMEANS=SED
PRINT CrMeans,SED; DECIMALS=2,3


" Analysis of copper concentrations "

" Remove missing values from analysis "
RESTRICT Cu; Cu.NE.c('*')

" One-way ANOVA of variate Cu with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=Cu
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] Cu


" One-way ANOVA of log10(Cu) with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=LOG10(Cu)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] log10(Cu)

" Formal test for homogeneity of variances "
VHOMOGENEITY [PRINT=test,variance; GROUPS=LandUse] log10(Cu)

" Save residuals from this model for use later "
AKEEP [RESIDUALS=ResidlogCu; RMETHOD=standardized]
" Save means and back-transform "
AKEEP TERMS=LandUse; MEANS=log10CuMeans; SEDMEANS=SED
PRINT log10CuMeans,SED; DEC=3,4
PRINT STRUCTURE=10**log10CuMeans; DECIMALS=3

" Other transformations "

" One-way ANOVA of sqrt(Cu) with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=sqrt(Cu)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" One-way ANOVA of 1/Cu (reciprocal transformation) with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=(1/Cu)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse


" Analysis of Zinc concentrations " 

" Remove missing values from analysis "
RESTRICT Zn; Zn.NE.c('*')

" One-way ANOVA of variate Zn with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=Zn
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] Zn


" One-way ANOVA of log10(Zn) with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=LOG10(Zn)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse

" Check within groups means and variances "
DESCRIBE [SELECT=nval,mean,var; GROUP=LandUse] log10(Zn)

" Formal test for homogeneity of variances "
VHOMOGENEITY [PRINT=test,variance; GROUPS=LandUse] log10(Zn)

" Save residuals from this model for use later "
AKEEP [RESIDUALS=ResidlogZn; RMETHOD=standardized]
" Save means and back-transform "
AKEEP TERMS=LandUse; MEANS=log10ZnMeans; SEDMEANS=SED
PRINT log10ZnMeans,SED; DEC=3,4
PRINT STRUCTURE=10**log10ZnMeans; DECIMALS=2


" Other transformations "

" One-way ANOVA of 1/Zn (reciprocal transformation) with residual plots "
TREATMENTSTRUCTURE LandUse
ANOVA [FPROBABILITY=yes] Y=(1/Zn)
APLOT [RMETHOD=standardized] METHOD=fittedvalues,normal,histogram,absresidual

" Boxplot of distribution of residuals "
AKEEP [RESIDUALS=Res; RMETHOD=standardized]
BOXPLOT Res; GROUPS=LandUse


" Investigating spatial correlation "

" Plot the co-ordinates of the spatial locations "
DGRAPH [TITLE='Co-ordinates of sampling locations'] Y=Y; X=X

" Form factor versions of X and Y taking into account that sample points are not on a perfect grid "
VARIATE IDENTIFIER=limits; VALUES=!(1,3,5,7,9,11...37)
CALCULATE limits=limits/8
GROUPS VECTOR=Y; FACTOR=fY; LIMITS=limits
GROUPS VECTOR=X; FACTOR=fX; LIMITS=limits

" Map data onto full grid (x ordered within y) "
TABULATE [CLASSIFICATION=fY,fX] ResidcubeCd; mean=Tmean
PRINT Tmean; DECIMALS=2; FIELD=6
" Get shade plot of residuals "
DSHADE [KEY=0; YORIENTATION=normal] Tmean; PEN=!(1)

" Plot residuals against previous residual within row or column "

" Set up window & pens for plotting "
CALCULATE Max = max(abs(ResidcubeCd))*1.1
YAXIS WINDOW=1; TITLE='Standardized residual'; LOWER=-Max; UPPER=Max
XAXIS WINDOW=1; TITLE='Previous standardized residual'; LOWER=-Max; UPPER=Max
PEN NUMBER=1; METHOD=point; LINESTYLE=0; SYMBOL=2; COLOUR='black'; CFILL='black'
PEN NUMBER=20; METHOD=line; LINESTYLE=1; SYMBOL=0; COLOUR='black'

" Plot residual vs residual from adjacent within-row sample (x-neighbours with same y-coordinate) "

" Map data onto full grid (x ordered within y) "
TABULATE [CLASSIFICATION=fY,fX] ResidcubeCd; mean=Tmean
DSHADE [YORIENTATION=normal] Tmean; PEN=!(4,2)
" Extract to variate with matching classifying factors "
VTABLE Tmean; VARIATE=xresCd; CLASS=!P(xfY,xfX)
" Get previous value "
CALCULATE prevYresCd = shift(xresCd)
" Get rid of any differences that go over two rows of grid "
CALCULATE prevYresCd = mvinsert(prevYresCd; xfY.eq.min(xfY))
" Plot points "
DGRAPH [WINDOW=1; KEYWINDOW=0; TITLE='Comparisons within rows'] Y=xresCd; prevYresCd; pen=1
" Add grid "
DGRAPH [WINDOW=1; KEYWINDOW=0; SCREEN=keep] Y=Max*!(0,0),!(-1,1); Max*!(-1,1),!(0,0); pen=20,20

" Plot residual vs residual from adjacent within-column sample (y-neighbours with same x-coordinate) "

" Map data onto full grid (y ordered within x) "
TABULATE [CLASSIFICATION=fX,fY] ResidcubeCd; mean=Tmean
" Extract to variate with matching classifying factors "
VTABLE Tmean; VARIATE=yresCd; CLASS=!P(yfX,yfY)
" Get previous value "
CALCULATE prevXresCd = shift(yresCd)
" Get rid of any differences that go over two columns of grid "
CALCULATE prevXresCd = mvinsert(prevXresCd; yfX.eq.min(yfX))
" Plot points "
DGRAPH [WINDOW=1; KEYWINDOW=0; TITLE='Comparisons within columns'] Y=yresCd; prevXresCd; pen=1
" Add grid "
DGRAPH [WINDOW=1; KEYWINDOW=0; SCREEN=keep] Y=Max*!(0,0),!(-1,1); Max*!(-1,1),!(0,0); pen=20,20

" Reset graphics to default set up "
YAXIS [RESET=yes] 1: XAXIS [RESET=yes] 1: PEN [RESET=yes] 1...20

" End of File "