QK Box 6.5

Medley and Clements (1998) sampled stations on six streams in a region of the Rocky Mountains of Colorado, USA known to be polluted by heavy metals. They recorded zinc concentration, and species richness and species diversity of the diatom community and proportion of diatom cells that were the early-successional species, Achnanthes minutissima. We’ll focus on the species diversity.

Diatom, Achnanthes longipes. Photo by Nic Craig.

These data were used in the first edition and the file is here.

Medley, C. N. & Clements, W. H. (1998). Responses of diatom communities to heavy metals in streams: The influence of longitudinal variation. Ecological Applications, 8, 631-44.

Preliminaries

First, make sure the core set of packages is loaded, then load any that are specific to this Box.

Load the required packages (car, sjstats, lmPerm)

Import medley data file (medley.csv)

medley <- read.csv("../data/medley.csv")
head(medley,10)

#make zinc is treated as a factor and arrange in sensible order
medley$zinc<-factor(medley$zinc)
medley$zinc<-fct_relevel(medley$zinc,c("BACK", "LOW", "MED", "HIGH"))

Fit model for diversity

medley.aov <- aov(diversity~zinc, data=medley)

Check diagnostics

plot(medley.aov)

Examine model results

glance(medley.aov)

tidy(medley.aov)

#summary(medley.aov)

Get effect sizes

effectsize(medley.aov)

For one-way between subjects designs, partial eta squared is equivalent to eta squared. Returning eta squared.

# Effect Size for ANOVA

Parameter | Eta2 |       95% CI
-------------------------------
zinc      | 0.28 | [0.04, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

effectsize(medley.aov, type = "omega")

For one-way between subjects designs, partial omega squared is equivalent to omega squared. Returning omega squared.

# Effect Size for ANOVA

Parameter | Omega2 |       95% CI
---------------------------------
zinc      |   0.21 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

effectsize(medley.aov, type = "f")

For one-way between subjects designs, partial eta squared is equivalent to eta squared. Returning eta squared.

# Effect Size for ANOVA

Parameter | Cohen's f |      95% CI
-----------------------------------
zinc      |      0.63 | [0.19, Inf]

- One-sided CIs: upper bound fixed at [Inf].

emmeans(medley.aov, ~zinc)

 zinc emmean    SE df lower.CL upper.CL
 BACK   1.80 0.165 30    1.461     2.13
 LOW    2.03 0.165 30    1.696     2.37
 MED    1.72 0.155 30    1.401     2.04
 HIGH   1.28 0.155 30    0.961     1.60

Confidence level used: 0.95 

Get default reference coding contrasts

summary.lm(medley.aov)

Call:
aov(formula = diversity ~ zinc, data = medley)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.03750 -0.22896  0.07986  0.33222  0.79750 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.79750    0.16478  10.909 5.81e-12 ***
zincLOW      0.23500    0.23303   1.008   0.3213    
zincMED     -0.07972    0.22647  -0.352   0.7273    
zincHIGH    -0.51972    0.22647  -2.295   0.0289 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4661 on 30 degrees of freedom
Multiple R-squared:  0.2826,    Adjusted R-squared:  0.2108 
F-statistic: 3.939 on 3 and 30 DF,  p-value: 0.01756

Check to make sure reference is Background

contrasts(medley$zinc) <- contr.treatment(4, base = 1)
medley2.aov<-aov(data=medley, diversity~zinc)
summary.lm(medley2.aov)

Call:
aov(formula = diversity ~ zinc, data = medley)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.03750 -0.22896  0.07986  0.33222  0.79750 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.79750    0.16478  10.909 5.81e-12 ***
zinc2        0.23500    0.23303   1.008   0.3213    
zinc3       -0.07972    0.22647  -0.352   0.7273    
zinc4       -0.51972    0.22647  -2.295   0.0289 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4661 on 30 degrees of freedom
Multiple R-squared:  0.2826,    Adjusted R-squared:  0.2108 
F-statistic: 3.939 on 3 and 30 DF,  p-value: 0.01756

Unplanned comparisons

Get Tukey pairwise comparisons and Dunnett test

med_emm <-emmeans(medley.aov, ~zinc)
contrast(med_emm, "tukey")

 contrast    estimate    SE df t.ratio p.value
 BACK - LOW   -0.2350 0.233 30  -1.008  0.7457
 BACK - MED    0.0797 0.226 30   0.352  0.9847
 BACK - HIGH   0.5197 0.226 30   2.295  0.1219
 LOW - MED     0.3147 0.226 30   1.390  0.5153
 LOW - HIGH    0.7547 0.226 30   3.333  0.0117
 MED - HIGH    0.4400 0.220 30   2.003  0.2096

P value adjustment: tukey method for comparing a family of 4 estimates 

contrast(med_emm, "dunnett")   #Note already have BACK as default reference group

 contrast    estimate    SE df t.ratio p.value
 LOW - BACK    0.2350 0.233 30   1.008  0.6140
 MED - BACK   -0.0797 0.226 30  -0.352  0.9493
 HIGH - BACK  -0.5197 0.226 30  -2.295  0.0761

P value adjustment: dunnettx method for 3 tests 

Run REGW test for comparison

This test uses the agricolae package

rw<-REGW.test(medley.aov, "zinc", alpha=0.05, console=FALSE, group=FALSE)
rw

$statistics
    MSerror Df     Mean      CV
  0.2172137 30 1.694118 27.5106

$parameters
  test name.t ntr alpha
  REGW   zinc   4  0.05

$regw
NULL

$means
     diversity       std r        se  Min  Max    Q25   Q50    Q75
BACK  1.797500 0.4852613 8 0.1647778 0.76 2.27 1.6575 1.935 2.0875
HIGH  1.277778 0.4268717 9 0.1553540 0.63 1.90 1.0400 1.250 1.4300
LOW   2.032500 0.4449960 8 0.1647778 1.40 2.83 1.7875 1.990 2.2300
MED   1.717778 0.5030104 9 0.1553540 0.80 2.19 1.6200 1.940 2.0600

$comparison
             difference pvalue signif.         LCL        UCL
BACK - HIGH  0.51972222 0.0570       . -0.01336042  1.0528049
BACK - LOW  -0.23500000 0.5594         -0.79329860  0.3232986
BACK - MED   0.07972222 0.9847         -0.53606192  0.6955064
HIGH - LOW  -0.75472222 0.0046      ** -1.28780487 -0.2216396
HIGH - MED  -0.44000000 0.1443         -0.99829860  0.1182986
LOW - MED    0.31472222 0.3191         -0.21836042  0.8478049

$groups
NULL

attr(,"class")
[1] "group"

Non-parametric approach

Do Kruskal-Wallis test

kruskal.test(diversity~zinc, data=medley)

    Kruskal-Wallis rank sum test

data:  diversity by zinc
Kruskal-Wallis chi-squared = 9, df = 3, p-value = 0.03

Randomization test (using lmPerm)

medley.aovp <- aovp(diversity~zinc, perm="Exact", data=medley)

[1] "Settings:  unique SS "

summary(medley.aovp)

Component 1 :
            Df R Sum Sq R Mean Sq Iter Pr(Prob)  
zinc2        3     2.57     0.856 4418    0.022 *
Residuals   30     6.52     0.217                
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Levene’s test

leveneTest(medley$diversity, medley$zinc)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3    0.02      1
      30