moss_anova.txt

				JFM  2/6/2012



ANOVA using moss cover data on 4 kinds of surfaces.  

Data collected by ESCI 340 students on 1/23/2012

In R.



ESCI 340 ANOVA & Tukey tests, using moss cover data

Calculations Friday 3 February and Monday 6 February 2012



> # Enter data

> moss.df <- c(10, 30, 0, 40, 5, 0, 25, 50, 0, 12)

> length(moss.df)   # Check that sample size = 10

[1] 10

> moss.bm <- c(60, 100, 15, 10, 80, 100, 100, 60, 70, 0)

> length(moss.bm)

[1] 10

> moss.bm

 [1]  60 100  15  10  80 100 100  60  70   0

> moss.rock <- c(30, 50, 100, 70, 100, 100, 7, 90,0, 60)

> moss.rock

 [1]  30  50 100  70 100 100   7  90   0  60

> length(moss.rock)

[1] 10

> moss.log <- c(30, 90, 25, 60, 0, 100, 45, 100, 10, 85)

> moss.log

 [1]  30  90  25  60   0 100  45 100  10  85

> length(moss.log)

[1] 10

> moss <- c(moss.df, moss.bm, moss.rock, moss.log)

> length(moss)

[1] 40

> moss

 [1]  10  30   0  40   5   0  25  50   0  12  60 100  15  10  80 100 100  60  70

[20]   0  30  50 100  70 100 100   7  90   0  60  30  90  25  60   0 100  45 100

[39]  10  85

> moss.surf <- rep(c("Dfir", "Maple", "Rock", "Log"), c(10, 10, 10, 10))

> moss.surf <- factor(moss.surf)

> plot(moss.surf, moss)



> # Use ANOVA function in R: 

>

> moss.aov <- aov(moss ~ moss.surf)

> anova(moss.aov)

Analysis of Variance Table



Response: moss

          Df Sum Sq Mean Sq F value  Pr(>F)  

moss.surf  3  12844  4281.4  3.6322 0.02182 *

Residuals 36  42435  1178.7                  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

>

> # ***********************************************************

> # Perform ANOVA calculations individually:

 

> moss.mean <- c(mean(moss.df), mean(moss.bm), mean(moss.rock), mean(moss.log))

> moss.mean

[1] 17.2 59.5 60.7 54.5

> moss.gss <- 3*var(moss.mean)

> moss.gss

[1] 1284.428

> moss.gss <- 3* 10 * var(moss.mean)

> moss.gss

[1] 12844.28

> moss.ess <- 9*(var(moss.df) + var(moss.bm) + var(moss.rock) + var(moss.log))

> moss.ess

[1] 42434.7

> moss.gms <- moss.gss / 3

> moss.gms

[1] 4281.425

> moss.ems <- moss.ess / (40-4)

> moss.ems

[1] 1178.742

> moss.f <- moss.gms / moss.ems

> moss.f

[1] 3.632200

> # Compare with F-statistic w/ 3,36 degrees of freedom

> # F_0.025(1),3,35 = 3.52     F_0.01(1),3,35 = 4.40

> # 0.025 > P > 0.01

> 1 - pf(moss.f, 3, 36)    # precise P-value

[1] 0.02181648



> # Tukey Multiple comparison test, to determine which means differ significantly

> moss.se <- sqrt(moss.ems / 10)

> moss.se

[1] 10.85699

> # q_0.001,30,4 = 5.156    so P < 0.001

> (mean(moss.df) - mean(moss.rock)) / moss.se

[1] -4.006636

> # q_0.05,30,4 = 3.845    so P < 0.05

> (mean(moss.df) - mean(moss.bm)) / moss.se

[1] -3.896109

> # P < 0.05

> (mean(moss.df) - mean(moss.log)) / moss.se

[1] -3.435576

>

*************************************************************************

Kruskal-Wallis NP ANOVA



> rank.log <- c(3.5, 10, 14.5, 17, 20, 24.5, 30, 31.5, 36.5, 36.5)

> sum(rank.df)

[1] 111.5

> sum(rank.bm)

[1] 241.5

> sum(rank.rock)

[1] 243

> sum(rank.log)

[1] 224

> r.sum <- sum(c(111.5^2, 241.5^2, 243^2, 224^2))/10

> r.sum

[1] 17997.95

> (12/(40*41))*r.sum - 3*41      # Kruskal-Wallis H-statistic

[1] 8.692317

> ti <- c(6, 3, 2, 3, 2, 4, 2, 2, 8)

> tsum <- sum(ti^3 - ti)

> tsum

[1] 846

> cfactor <- 1 - (tsum/(40^3 - 40))

> cfactor

[1] 0.986773

> hc <- ((12/(40*41))*r.sum - 3*41) / cfactor

> hc

[1] 8.808832

> # Compare with Chi-squared statistic w/ k-1 degrees of freedom

> # Chi-sq_0.05,3 = 7.815    Chi-sq_0.025,3 = 9.345

> # 0.05 > P > 0.025

> 1 - pchisq(hc, df = 3)       # precise P-value

[1] 0.03194357