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)Leaf Water Potential
Leaf Water Potential analysis
- Prepare data
- Visualize
Leaf Wat Pot 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 = "lwatpot",
groupvars = c("temp", "eco")
)
library(ggplot2)
ggplot(
data = summary_stat,
aes(x = temp, y = lwatpot, colour = eco)
) +
geom_errorbar(aes(ymin = lwatpot - se, ymax = lwatpot + se), # add error bars
width = 0.1 # width of error bars
) +
geom_line(aes(group=eco)) +
geom_point() +
labs(y = "lwatpot")+
scale_color_manual(values=c("deepskyblue3", "darkgoldenrod"))
Boxplots
library(ggplot2)
ggplot(data=yp_data, aes(x=temp,y=lwatpot,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("lwatpot"))+
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: lwatpot ~ temp + eco
interaction: lwatpot ~ temp * eco
-> Identify best model using AIC
nointeract <- aov(lwatpot ~ temp + eco, data = yp_data)
interaction <- aov(lwatpot ~ 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
interaction 7 -54.19 0.00 1 1 37.60
no interaction 5 -32.55 21.64 0 1 22.94-> “Interaction” model is the best fit.
- Run ANOVA model
interaction <- aov(lwatpot ~ temp * eco, data = yp_data)
summary(interaction) Df Sum Sq Mean Sq F value Pr(>F)
temp 2 1.2140 0.6070 178.38 1.36e-12 ***
eco 1 0.2817 0.2817 82.78 3.75e-08 ***
temp:eco 2 0.1465 0.0732 21.52 1.69e-05 ***
Residuals 18 0.0612 0.0034
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Response to both ecotype and temperature is significant.
-> Run a Tukey post-hoc test.
TukeyHSD(interaction) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = lwatpot ~ temp * eco, data = yp_data)
$temp
diff lwr upr p adj
27.5°C-20°C -0.51250 -0.586938092 -0.4380619 0.0000000
35°C-20°C -0.43125 -0.505688092 -0.3568119 0.0000000
35°C-27.5°C 0.08125 0.006811908 0.1556881 0.0312056
$eco
diff lwr upr p adj
LR-HR 0.2166667 0.1666343 0.266699 0
$`temp:eco`
diff lwr upr p adj
27.5°C:HR-20°C:HR -0.3500 -0.48108725 -0.21891275 0.0000014
35°C:HR-20°C:HR -0.2625 -0.39358725 -0.13141275 0.0000686
20°C:LR-20°C:HR 0.4375 0.30641275 0.56858725 0.0000000
27.5°C:LR-20°C:HR -0.2375 -0.36858725 -0.10641275 0.0002314
35°C:LR-20°C:HR -0.1625 -0.29358725 -0.03141275 0.0104777
35°C:HR-27.5°C:HR 0.0875 -0.04358725 0.21858725 0.3202302
20°C:LR-27.5°C:HR 0.7875 0.65641275 0.91858725 0.0000000
27.5°C:LR-27.5°C:HR 0.1125 -0.01858725 0.24358725 0.1178189
35°C:LR-27.5°C:HR 0.1875 0.05641275 0.31858725 0.0029060
20°C:LR-35°C:HR 0.7000 0.56891275 0.83108725 0.0000000
27.5°C:LR-35°C:HR 0.0250 -0.10608725 0.15608725 0.9892050
35°C:LR-35°C:HR 0.1000 -0.03108725 0.23108725 0.1996196
27.5°C:LR-20°C:LR -0.6750 -0.80608725 -0.54391275 0.0000000
35°C:LR-20°C:LR -0.6000 -0.73108725 -0.46891275 0.0000000
35°C:LR-27.5°C:LR 0.0750 -0.05608725 0.20608725 0.4792908- The difference between 27.5 and 35deg is stat significant. Furthermore, the difference for high rainfall ecotype is significant between 20 and 27.5deg, and between 20 deg and 35deg.
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(lwatpot ~ treatment, yp_data)
anova(yp_model)Analysis of Variance Table
Response: lwatpot
Df Sum Sq Mean Sq F value Pr(>F)
treatment 5 1.64208 0.32842 96.514 2.36e-12 ***
Residuals 18 0.06125 0.00340
---
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 lwatpot for 20deg = mean lwatpot for 35deg |
| 3: Contrast ‘temp quadratic’ (1,-2,1,1,-2,1) ->H0 = Mean lwatpot at 27.5deg = Mean lwatpot (20deg,35deg) |
| 4: Contrast ‘ecotype * temp linear’ (1,0,-1,-1,0,1) -> H0 = Linear component of lwatpot ~ temp is the same for both ecotypes |
| 5: Contrast ‘ecotype * temp quadratic’ (1,-2,1,-1,2,-1) -> H0 = Quadratic component of lwatpot ~ 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(lwatpot ~ 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 1.6421 0.3284 96.514 2.36e-12 ***
treatment: Eco 1 0.2817 0.2817 82.776 3.75e-08 ***
treatment: temp lin 1 0.7439 0.7439 218.617 1.64e-11 ***
treatment: temp quad 1 0.4701 0.4701 138.138 7.05e-10 ***
treatment: e*t lin 1 0.1139 0.1139 33.474 1.75e-05 ***
treatment: e*t quad 1 0.0326 0.0326 9.566 0.00628 **
Residuals 18 0.0612 0.0034
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant response of lwatpot to ecotype.
The response of lwatpot to temperature may have a quadratic component that depends on ecotype.