yp_data <- read.csv("SDlluteus_rev.csv")
yp_data$eco <- as.factor(yp_data$eco)
yp_data$temp <-as.factor(yp_data$temp)
levels(yp_data$temp) <- c("20°C","27.5°C","35°C")
yp_data$treatment <- as.factor(yp_data$treatment)Photosynthesis rate
Photosynthesis Rate analysis
- Prepare data
- Visualize
Photosynthesis rate mean +/- SE
library(Rmisc)Loading required package: latticeLoading required package: plyrlibrary(plyr)
library(lattice)
# compute mean and standard error of the mean by subgroup
summary_stat <- summarySE(yp_data,
measurevar = "photoS",
groupvars = c("temp", "eco")
)
library(ggplot2)
ggplot(
data = summary_stat,
aes(x = temp, y = photoS, colour = eco)
) +
geom_errorbar(aes(ymin = photoS - se, ymax = photoS + se), # add error bars
width = 0.1 # width of error bars
) +
geom_line(aes(group=eco)) +
geom_point() +
labs(y = "photoS")+
scale_color_manual(values=c("deepskyblue3", "darkgoldenrod"))
Boxplots
library(ggplot2)
ggplot(data=yp_data, aes(x=temp,y=photoS,fill=eco))+
geom_boxplot(outlier.shape = NA)+
geom_point(position=position_jitterdodge(),size=3,alpha=0.5)+
theme(
plot.tag = element_text(size = 16),
panel.border = element_rect(colour = "black", fill=NA, size=1),
panel.background = element_rect(fill = "white"),
axis.title.x = element_blank(),
axis.text.x = element_text(colour = "black", size=12),
axis.title.y = element_text(colour = "black", size=14),
axis.text.y = element_text(colour = "black", size=12)
)+
labs(y=expression("photoS"))+
scale_fill_manual(values=c("deepskyblue3", "darkgoldenrod"))Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.
“Interaction” v. “No interaction” model
no interaction: photoS ~ temp + eco
interaction: photoS ~ temp * eco
-> Identify best model using AIC
nointeract <- aov(photoS ~ temp + eco, data = yp_data)
interaction <- aov(photoS ~ temp * eco, data = yp_data)
#find best model (best listed first via aictab)
library(AICcmodavg)
model.set <- list(nointeract, interaction)
model.names <- c("no interaction", "interaction")
aictab(model.set, modnames = model.names)
Model selection based on AICc:
K AICc Delta_AICc AICcWt Cum.Wt LL
no interaction 5 142.34 0.00 0.97 0.97 -64.50
interaction 7 149.17 6.84 0.03 1.00 -64.09-> “no Interaction” model is the best fit.
- Run ANOVA model
nointeract <- aov(photoS ~ temp + eco, data = yp_data)
summary(nointeract) Df Sum Sq Mean Sq F value Pr(>F)
temp 2 103.9 52.0 3.424 0.0526 .
eco 1 452.5 452.5 29.824 2.4e-05 ***
Residuals 20 303.5 15.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The ecotype has a statistically significant effect on photosynthesis rate (p=2.4e-05).
Because the p-value for temperature is close to 0.05, we decided to run a A Tukey post-hoc test.
-> Run a Tukey post-hoc test.
TukeyHSD(nointeract) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = photoS ~ temp + eco, data = yp_data)
$temp
diff lwr upr p adj
27.5°C-20°C 2.124989 -2.802511 7.0524897 0.5304162
35°C-20°C -2.949742 -7.877243 1.9777585 0.3055437
35°C-27.5°C -5.074731 -10.002232 -0.1472306 0.0428359
$eco
diff lwr upr p adj
LR-HR 8.684521 5.367332 12.00171 2.4e-05- The test confirmed significant pairwise differences between the 2 ecotypes (+ 8.68 under 28830), and a statistically significant difference between 27.5 and 35°C (-5.07 under 35°C).
To explore this further, we combined ecotype treatment (High Rainfall + Low Rainfall) and temperature Treatment (20, 27.5, 35) to create 6 treatments: HR20, HR27.5, HR35, L20, L27.5, L35. We used orthogonal polynomial contrasts to investigate differences between levels of treatment.
- Apply orthogonal polynomial contrasts
-> Run ANOVA:
yp_data$treatment <- as.factor(yp_data$treatment)
str(yp_data$treatment) Factor w/ 6 levels "H20","H275","H35",..: 4 4 4 4 1 1 1 1 5 5 ...# Run ANOVA for outcome
yp_model <- aov(photoS ~ treatment, yp_data)
anova(yp_model)Analysis of Variance Table
Response: photoS
Df Sum Sq Mean Sq F value Pr(>F)
treatment 5 566.77 113.355 6.9606 0.0008857 ***
Residuals 18 293.14 16.285
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1->Set up orthogonal contrasts:
With 6 treatments (H20, H275, H35, L20, L275, L35), 5 orthogonal comparisons (i.e., statistically independent) are possible:
-> Assign contrast coef based on research question:
| 1: Contrast ‘ecotype’ (1,1,1,-1,-1,-1) -> H0 = Mean photoS is the same for both ecotypes |
| 2: Contrast ‘temp linear’ (1,0,-1,1,0,-1) -> H0 = Mean photoS for 20deg = mean photoS for 35deg |
| 3: Contrast ‘temp quadratic’ (1,-2,1,1,-2,1) ->H0 = Mean photoS at 27.5deg = Mean photoS (20deg,35deg) |
| 4: Contrast ‘ecotype * temp linear’ (1,0,-1,-1,0,1) -> H0 = Linear component of photoS ~ temp is the same for both ecotypes |
| 5: Contrast ‘ecotype * temp quadratic’ (1,-2,1,-1,2,-1) -> H0 = Quadratic component of photoS ~ temp is the same for both ecotypes |
-> Create contrast matrix based on coef:
contrastmatrix <- cbind( c(1,1,1,-1,-1,-1),
c(1,0,-1,1,0,-1),
c(1,-2,1,1,-2,1),
c(1,0,-1,-1,0,1),
c(1,-2,1,-1,2,-1))-> Apply matrix to treatments:
contrasts(yp_data$treatment) <- contrastmatrix
yp_data$treatment [1] L20 L20 L20 L20 H20 H20 H20 H20 L275 L275 L275 L275 H275 H275 H275
[16] H275 L35 L35 L35 L35 H35 H35 H35 H35
attr(,"contrasts")
[,1] [,2] [,3] [,4] [,5]
H20 1 1 1 1 1
H275 1 0 -2 0 -2
H35 1 -1 1 -1 1
L20 -1 1 1 -1 -1
L275 -1 0 -2 0 2
L35 -1 -1 1 1 -1
Levels: H20 H275 H35 L20 L275 L35-> Run contrast analysis:
yp_contrast_aov <- aov(photoS ~ treatment, yp_data)
summary(yp_contrast_aov, split = list(treatment = list("Eco" = 1,
"temp lin"= 2, "temp quad" = 3, "e*t lin" = 4, "e*t quad" = 5))) Df Sum Sq Mean Sq F value Pr(>F)
treatment 5 566.8 113.4 6.961 0.000886 ***
treatment: Eco 1 452.5 452.5 27.787 5.18e-05 ***
treatment: temp lin 1 34.8 34.8 2.137 0.161010
treatment: temp quad 1 69.1 69.1 4.244 0.054148 .
treatment: e*t lin 1 9.0 9.0 0.550 0.467927
treatment: e*t quad 1 1.4 1.4 0.084 0.774736
Residuals 18 293.1 16.3
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- There is a significant response of photosynthesis rate to ecotype.