SEAGRASS
Alan
December 2, 2014
First read in SEAGRASS.csv provided courtesy of Maite from the study Elso, Garc-Jim’enez, and Robaina (2012).
SEAGRASS <- repmis::source_data("https://raw.github.com/alanarnholt/Data/master/SEAGRASS.csv")Downloading data from: https://raw.github.com/alanarnholt/Data/master/SEAGRASS.csv
SHA-1 hash of the downloaded data file is:
58d3636d52e812be37ef6d7b14f4bd25baf79ccchead(SEAGRASS) fluorescence spermidine salinity
1 0.6541250 NO 05PSU
2 0.5437917 NO 05PSU
3 0.5454167 NO 05PSU
4 0.5878333 NO 05PSU
5 0.7058750 YES 05PSU
6 0.6670000 YES 05PSUstr(SEAGRASS)'data.frame': 32 obs. of 3 variables:
$ fluorescence: num 0.654 0.544 0.545 0.588 0.706 ...
$ spermidine : Factor w/ 2 levels "NO","YES": 1 1 1 1 2 2 2 2 1 1 ...
$ salinity : Factor w/ 4 levels "05PSU","11PSU",..: 1 1 1 1 1 1 1 1 2 2 ...Note: this is a two-factorial design…not a one-way ANOVA…if it was analyzed as a one-way ANOVA originally, then it slipped through the review process. There are two factors: spermidine with two levels and salinity with four levels. That is, it is a completely randomized, fixed effects, 2 \(\times\) 4 factorial experiment.
Note to self: spermidine regulates plant growth and the presence of spermidine is known to help plants tolerate salinity so one might argue that it is a block…although one does not test for the effects of blocks…the lower the fluorescence value the more stress a plant is under.
model <- aov(fluorescence ~ spermidine*salinity, data = SEAGRASS)
summary(model) Df Sum Sq Mean Sq F value Pr(>F)
spermidine 1 0.01123 0.01123 8.270 0.00832 **
salinity 3 0.14379 0.04793 35.285 5.83e-09 ***
spermidine:salinity 3 0.00680 0.00227 1.668 0.20044
Residuals 24 0.03260 0.00136
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(TukeyHSD(model, which = c("spermidine", "salinity")))
with(data = SEAGRASS,
tapply(fluorescence, list(spermidine, salinity), mean)
) 05PSU 11PSU 18PSU 36PSU
NO 0.5827917 0.7074583 0.7736958 0.7818125
YES 0.6700104 0.7239271 0.8026875 0.7990209Checking model assumptions- normality of errors
library(car)
qqPlot(model, envelope = FALSE)
shapiro.test(rstandard(model))
Shapiro-Wilk normality test
data: rstandard(model)
W = 0.9603, p-value = 0.2791Checking constant variance
leveneTest(model)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 7 0.8973 0.5243
24 Is there any interaction?
library(ggplot2)
ggplot(data = SEAGRASS, aes(x = salinity, y = fluorescence, color = spermidine,
group = spermidine, linetype = spermidine)) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line") +
theme_bw()
ggplot(data = SEAGRASS, aes(x = spermidine, y = fluorescence, color = salinity,
group = salinity, linetype = salinity)) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line") +
theme_bw()
library(PASWR2)
with(data = SEAGRASS,
twoway.plots(fluorescence, spermidine, salinity)
)
Model effects and means
model.tables(model, type = "effects", se = TRUE,
cterms = c("spermidine", "salinity"))Tables of effects
spermidine
spermidine
NO YES
-0.018736 0.018736
salinity
salinity
05PSU 11PSU 18PSU 36PSU
-0.10377 -0.01448 0.05802 0.06024
Standard errors of effects
spermidine salinity
0.009214 0.013031
replic. 16 8model.tables(model, type = "means", se = TRUE,
cterms = c("spermidine", "salinity"))Tables of means
Grand mean
0.7301755
spermidine
spermidine
NO YES
0.7114 0.7489
salinity
salinity
05PSU 11PSU 18PSU 36PSU
0.6264 0.7157 0.7882 0.7904
Standard errors for differences of means
spermidine salinity
0.01303 0.01843
replic. 16 8Determine power assuming the true means for the eight treatments are:
MEANS <- c(0.6, 0.70, 0.72, 0.73, 0.77, 0.8, 0.79, 0.8)
M <- matrix(MEANS, byrow = TRUE, nrow = 2)
dimnames(M) <- list(c("NO", "YES"), c("05PSU", "11PSU", "18PSU", "36PSU"))
M 05PSU 11PSU 18PSU 36PSU
NO 0.60 0.7 0.72 0.73
YES 0.77 0.8 0.79 0.80if there are n = 6 measurements per treatment assuming \(\sigma = 0.1\).
First we must find \(\lambda\) for each hypothesis.
alpha <- 0.05
Sigma <- 0.1
n <- 6
a <- 2
b <- 4
Y <- rep(MEANS, each = n)
Tsper <- factor(rep(rep(c("NO", "YES"), each = n), b))
Tsper [1] NO NO NO NO NO NO YES YES YES YES YES YES NO NO NO NO NO
[18] NO YES YES YES YES YES YES NO NO NO NO NO NO YES YES YES YES
[35] YES YES NO NO NO NO NO NO YES YES YES YES YES YES
Levels: NO YESTsali <- factor(rep(c("05PSU","11PSU","18PSU","36PSU"), each = n*a))
NDF <- data.frame(Y, Tsper, Tsali)
str(NDF)'data.frame': 48 obs. of 3 variables:
$ Y : num 0.6 0.6 0.6 0.6 0.6 0.6 0.7 0.7 0.7 0.7 ...
$ Tsper: Factor w/ 2 levels "NO","YES": 1 1 1 1 1 1 2 2 2 2 ...
$ Tsali: Factor w/ 4 levels "05PSU","11PSU",..: 1 1 1 1 1 1 1 1 1 1 ...summary(aov(Y ~ Tsper*Tsali, data = NDF)) Df Sum Sq Mean Sq F value Pr(>F)
Tsper 1 0.01687 0.01687 3.975e+30 <2e-16 ***
Tsali 3 0.16043 0.05348 1.260e+31 <2e-16 ***
Tsper:Tsali 3 0.01642 0.00547 1.290e+30 <2e-16 ***
Residuals 40 0.00000 0.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1SShypTsper <- summary(aov(Y ~ Tsper*Tsali, data = NDF))[[1]][1, 2]
SShypTsper[1] 0.016875lambdaTsper <- SShypTsper/Sigma^2
lambdaTsper[1] 1.6875critF1 <- qf(1 - alpha, a-1, a*b*(n-1))
critF1[1] 4.084746TheoPowerH1 <- pf(critF1, a-1, a*b*(n-1), lambdaTsper, lower = FALSE)
TheoPowerH1[1] 0.2450903SShypTsali <- summary(aov(Y ~ Tsper*Tsali, data = NDF))[[1]][2, 2]
SShypTsali[1] 0.160425lambdaTsali <- SShypTsali/Sigma^2
lambdaTsali[1] 16.0425critF2 <- qf(1 - alpha, b-1, a*b*(n-1))
critF2[1] 2.838745TheoPowerH2 <- pf(critF2, b-1, a*b*(n-1), lambdaTsali, lower = FALSE)
TheoPowerH2[1] 0.9078258Note that with 6 observations per treatment the power to detect differences in PSU levels is 0.9078258.
References
Elso, M Zarranz, P Garc-Jim’enez, and RR Robaina. 2012. “Endogenous Polyamine Content and Photosynthetic Performance Under Hypo-Osmotic Conditions Reveal Cymodocea Nodosa as an Obligate Halophyte.” Aquat. Biol. 17 (1). Inter-Research Science Center: 7–17. doi:10.3354/ab00454. http://dx.doi.org/10.3354/ab00454.