setwd("C:/Users/theco/Desktop")
library(plotrix)
library(car)
## Warning: package 'car' was built under R version 3.3.3
Data <- read.csv("oxygenlightF.csv", header = TRUE)
Data
## Treatment chamber O2.flux cockles ulva
## 1 black Dark 1.1500000
## 2 black Light 1.1400000
## 3 blue Dark 1.1300000 cockles
## 4 blue Light 1.1800000 cockles
## 5 green Dark 1.5100000 cockles ulva
## 6 green Light 1.1000000 cockles ulva
## 7 yellow Dark -1.5700000 ulva
## 8 yellow Light 0.1200000 ulva
## 9 black Dark -3.6400000
## 10 black Light 1.4600000
## 11 blue Dark -0.8100000 cockles
## 12 blue Light -1.0300000 cockles
## 13 green Dark 1.4285714 cockles ulva
## 14 green Light 1.5259740 cockles ulva
## 15 blue Dark 1.3900000 cockles
## 16 blue Light 0.1200000 cockles
## 17 green Dark -3.2600000 cockles ulva
## 18 green Light -1.0800000 cockles ulva
## 19 yellow Dark -3.6800000 ulva
## 20 yellow Light -0.4200000 ulva
## 21 black Dark -1.9100000
## 22 black Light 0.0200000
## 23 blue Dark 0.4600000 cockles
## 24 blue Light 0.9700000 cockles
## 25 green Dark -3.7700000 cockles ulva
## 26 green Light -2.6700000 cockles ulva
## 27 yellow Dark -3.6600000 ulva
## 28 yellow Light 0.8200000 ulva
## 29 black Dark -2.6900000
## 30 black Light 0.5000000
## 31 blue Dark -0.2300000 cockles
## 32 blue Light -0.2800000 cockles
## 33 green Dark -0.3900000 cockles ulva
## 34 green Light 0.2400000 cockles ulva
## 35 yellow Dark -3.2100000 ulva
## 36 yellow Light -1.1000000 ulva
## 37 black Dark -1.2000000
## 38 black Light 1.0700000
## 39 blue Dark -0.6400000 cockles
## 40 blue Light 0.6700000 cockles
## 41 green Dark 0.6000000 cockles ulva
## 42 green Light -1.1500000 cockles ulva
## 43 yellow Dark -1.4100000 ulva
## 44 yellow Light 1.6200000 ulva
## 45 black Dark -2.6900000
## 46 black Light 2.4400000
## 47 blue Dark -2.7800000 cockles
## 48 blue Light -1.0600000 cockles
## 49 green Dark -0.5600000 cockles ulva
## 50 green Light -0.4100000 cockles ulva
## 51 yellow Dark -2.0616883 ulva
## 52 yellow Light 0.8847403 ulva
summary(Data)
## Treatment chamber O2.flux cockles ulva
## black :12 Dark :26 Min. :-3.7700 :24 :26
## blue :14 Light:26 1st Qu.:-1.4500 cockles:28 ulva:26
## green :14 Median :-0.3350
## yellow:12 Mean :-0.4964
## 3rd Qu.: 0.9950
## Max. : 2.4400
head(Data)
## Treatment chamber O2.flux cockles ulva
## 1 black Dark 1.15
## 2 black Light 1.14
## 3 blue Dark 1.13 cockles
## 4 blue Light 1.18 cockles
## 5 green Dark 1.51 cockles ulva
## 6 green Light 1.10 cockles ulva
#detach(Data)
attach(Data)
hist(O2.flux)
# black = nothing; blue = +c-u, yellow = -c+u, green = +c+u
### Cockles, Ulva and dO
# Interaction effect insignificant; therefore -u and +u data can be combined. See "MARI301 R 2017.R"
## ALL DATA (cockles vs no cockles?)
all.yc <- O2.flux[Data$cockles == "cockles"]
all.nc <- O2.flux[Data$cockles != "cockles"]
all.means = c(mean(all.nc), mean(all.yc))
mp <- barplot(all.means, axes=FALSE, axisnames=FALSE, ylim=c(-1.5,.5), main="All chambers, no ulva", xlab="# cockles", ylab="Oxygen flux (mM/L)")
axis(1, labels=c("0","300"), at=mp)
axis(2, at=seq(-.5,1.5,by=.5))
box()
std.error(all.yc)
## [1] 0.2782361
std.error(all.nc)
## [1] 0.3873289
SE <- matrix(c(std.error(all.nc), std.error(all.yc)), 2)
segments(mp, all.means-SE, mp, all.means+SE, lwd=2)
segments(mp-0.1, all.means-SE, mp+0.1, all.means-SE, lwd=2)
segments(mp-0.1, all.means+SE, mp+0.1, all.means+SE, lwd=2)
abline(h=0, lty=2)
# stats
leveneTest(O2.flux, Treatment)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.3792 0.2604
## 48
# group 3 0.9292 0.4432
shapiro.test(O2.flux)
##
## Shapiro-Wilk normality test
##
## data: O2.flux
## W = 0.9387, p-value = 0.009922
# W = 0.96253, p-value = 0.4438
model1 <- lm(O2.flux ~ Treatment)
summary(model1)
##
## Call:
## lm(formula = O2.flux ~ Treatment)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2782 -0.9333 0.0918 1.3023 2.8025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.3625 0.4867 -0.745 0.460
## Treatmentblue 0.2975 0.6633 0.449 0.656
## Treatmentgreen -0.1293 0.6633 -0.195 0.846
## Treatmentyellow -0.7764 0.6883 -1.128 0.265
##
## Residual standard error: 1.686 on 48 degrees of freedom
## Multiple R-squared: 0.05391, Adjusted R-squared: -0.005222
## F-statistic: 0.9117 on 3 and 48 DF, p-value: 0.4424
# Residual standard error: 1.084 on 22 degrees of freedom
cockles <- Data$cockles == "cockles"
# cockles2 <- Data$Treatment == "green" | Data$Treatment == "blue"
ulva <- Data$ulva == "ulva"
model2 <- lm(O2.flux~cockles + ulva + cockles:ulva)
Anova(model2, type="III")
## Anova Table (Type III tests)
##
## Response: O2.flux
## Sum Sq Df F value Pr(>F)
## (Intercept) 1.577 1 0.5547 0.4600
## cockles 0.572 1 0.2012 0.6558
## ulva 3.617 1 1.2724 0.2649
## cockles:ulva 0.395 1 0.1389 0.7110
## Residuals 136.447 48
model3 <- lm(O2.flux~cockles+ulva)
Anova(model3, type="II")
## Anova Table (Type II tests)
##
## Response: O2.flux
## Sum Sq Df F value Pr(>F)
## cockles 2.883 1 1.0322 0.3146
## ulva 4.497 1 1.6104 0.2104
## Residuals 136.842 49
plot(model3)
# Conclusion: addition of cockles is insignificant in the absense of ulva when dark and light are taken together
# Next step: is dark significantly different than light? Should they be separately analyzed?
## ALL DATA Dark vs Light
all.d <- O2.flux[Data$chamber == "Dark"]
all.l <- O2.flux[Data$chamber == "Light"]
all.means.dvl = c(mean(all.d), mean(all.l))
mp.dvl <- barplot(all.means.dvl, axes=FALSE, axisnames=FALSE, ylim=c(-2,1), main="Dark versus light, all chambers", xlab="Chamber type", ylab="Oxygen flux (mM/L)", col=c("white", "gray"))
axis(1, labels=c("Dark", "Light"), at=mp)
axis(2, at=seq(-2,1))
box()
SE.dvl <- matrix(c(std.error(all.d), std.error(all.l)), 2)
segments(mp.dvl, all.means.dvl-SE.dvl, mp.dvl, all.means.dvl+SE.dvl, lwd=2)
segments(mp.dvl-0.1, all.means.dvl-SE.dvl, mp.dvl+0.1, all.means.dvl-SE.dvl, lwd=2)
segments(mp.dvl-0.1, all.means.dvl+SE.dvl, mp.dvl+0.1, all.means.dvl+SE.dvl, lwd=2)
abline(h=0, lty=2)
cockles <- Data$cockles == "cockles"
light <- Data$chamber == "Light"
model2 <- lm(O2.flux~cockles + light + cockles:light)
Anova(model2, type="III")
## Anova Table (Type III tests)
##
## Response: O2.flux
## Sum Sq Df F value Pr(>F)
## (Intercept) 58.838 1 31.611 9.413e-07 ***
## cockles 20.735 1 11.140 0.001639 **
## light 51.411 1 27.621 3.358e-06 ***
## cockles:light 22.485 1 12.080 0.001092 **
## Residuals 89.343 48
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Thus we have established dark is significantly different than light!
#DARK
all.yc.d <- O2.flux[Data$cockles == "cockles" & Data$chamber == "Dark"]
all.nc.d <- O2.flux[Data$cockles != "cockles" & Data$chamber == "Dark"]
all.means.d = c(mean(all.nc.d), mean(all.yc.d))
mp.d <- barplot(all.means.d, axes=FALSE, axisnames=FALSE, ylim=c(-3,1.5), xlab="Presence of cockles", ylab="Oxygen flux (mM/L)")
axis(1, labels=c("No (-c)","Yes (+c)"), at=mp)
axis(2, at=seq(-3,1.5,by=0.5))
box()
SE.d <- matrix(c(std.error(all.nc.d), std.error(all.yc.d)), 2)
segments(mp.d, all.means.d-SE.d, mp.d, all.means.d+SE.d, lwd=2)
segments(mp.d-0.1, all.means.d-SE.d, mp.d+0.1, all.means.d-SE.d, lwd=2)
segments(mp.d-0.1, all.means.d+SE.d, mp.d+0.1, all.means.d+SE.d, lwd=2)
abline(h=0, lty=2)
all.model.d <- lm(O2.flux[Data$chamber=="Dark"] ~ cockles[Data$chamber=="Dark"])
summary(all.model.d) # p = .00871
##
## Call:
## lm(formula = O2.flux[Data$chamber == "Dark"] ~ cockles[Data$chamber ==
## "Dark"])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3470 -0.8657 0.0928 0.9815 3.3643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.2143 0.4602 -4.812 6.7e-05
## cockles[Data$chamber == "Dark"]TRUE 1.7913 0.6271 2.856 0.00871
##
## (Intercept) ***
## cockles[Data$chamber == "Dark"]TRUE **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.594 on 24 degrees of freedom
## Multiple R-squared: 0.2537, Adjusted R-squared: 0.2226
## F-statistic: 8.159 on 1 and 24 DF, p-value: 0.008706
#LIGHT
all.yc.l <- O2.flux[Data$cockles == "cockles" & Data$chamber == "Light"]
all.nc.l <- O2.flux[Data$cockles != "cockles" & Data$chamber == "Light"]
all.means.l = c(mean(all.nc.l), mean(all.yc.l))
mp.l <- barplot(all.means.l, axes=FALSE, axisnames=FALSE, ylim=c(-3,1.5), xlab="Presence of cockles", ylab="Oxygen flux (mM/L)")
axis(1, labels=c("No (-c)","Yes (+c)"), at=mp)
axis(2, at=seq(-3,1.5,by=0.5))
box()
SE.l <- matrix(c(std.error(all.nc.l), std.error(all.yc.l)), 2)
segments(mp.l, all.means.l-SE.l, mp.l, all.means.l+SE.l, lwd=2)
segments(mp.l-0.1, all.means.l-SE.l, mp.l+0.1, all.means.l-SE.l, lwd=2)
segments(mp.l-0.1, all.means.l+SE.l, mp.l+0.1, all.means.l+SE.l, lwd=2)
abline(h=0, lty=2)
all.model.l <- lm(O2.flux[Data$chamber=="Light"] ~ cockles[Data$chamber=="Light"])
summary(all.model.l)
##
## Call:
## lm(formula = O2.flux[Data$chamber == "Light"] ~ cockles[Data$chamber ==
## "Light"])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.5361 -0.8453 0.1395 0.7897 1.7271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.7129 0.3138 2.272 0.0323
## cockles[Data$chamber == "Light"]TRUE -0.8468 0.4276 -1.980 0.0592
##
## (Intercept) *
## cockles[Data$chamber == "Light"]TRUE .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.087 on 24 degrees of freedom
## Multiple R-squared: 0.1405, Adjusted R-squared: 0.1046
## F-statistic: 3.922 on 1 and 24 DF, p-value: 0.05924
# Adding cockles was significant with respect to dO for both light chambers and dark chambers regardless of ulva present
# +c = nutrient flux +/-?
### Cockles, Ulva and chambers
## Interaction effects, 3way anova
# 3-way anova code taken from www.quantide.com/wp-content/uploads/2017/02/Three-way-Anova-with-R.pdf
O2.m1 <- aov(O2.flux~cockles*ulva*chamber, data=Data)
summary(O2.m1)
## Df Sum Sq Mean Sq F value Pr(>F)
## cockles 1 2.88 2.883 1.502 0.226895
## ulva 1 4.50 4.497 2.343 0.132994
## chamber 1 29.51 29.511 15.376 0.000305 ***
## cockles:ulva 1 0.39 0.395 0.206 0.652369
## cockles:chamber 1 22.48 22.485 11.715 0.001351 **
## ulva:chamber 1 0.00 0.000 0.000 0.988393
## cockles:ulva:chamber 1 0.00 0.000 0.000 0.995841
## Residuals 44 84.45 1.919
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
O2.m2 <- update(O2.m1, . ~ . - cockles:ulva:chamber)
Anova(O2.m1, O2.m2, type="III")
## Anova Table (Type III tests)
##
## Response: O2.flux
## Sum Sq Df F value Pr(>F)
## (Intercept) 20.093 1 10.4690 0.0023076 **
## cockles 8.464 1 4.4098 0.0415002 *
## ulva 1.772 1 0.9234 0.3418351
## chamber 25.843 1 13.4644 0.0006534 ***
## cockles:ulva 0.193 1 0.1005 0.7527321
## cockles:chamber 11.277 1 5.8754 0.0195319 *
## ulva:chamber 0.000 1 0.0002 0.9890624
## cockles:ulva:chamber 0.000 1 0.0000 0.9958409
## Residuals 84.450 45
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p-value: .0415 *, 0.3418, 0.0007 ***, 0.7527, 0.0195 *, 0.9891, 0.9958
plot(O2.m2) # residuals
