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)Transpiration
Transpiration analysis
- Prepare data
- Visualize
Transpiration 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 = "Trmmol",
groupvars = c("temp", "eco")
)
library(ggplot2)
ggplot(
data = summary_stat,
aes(x = temp, y = Trmmol, colour = eco)
) +
geom_errorbar(aes(ymin = Trmmol - se, ymax = Trmmol + se), # add error bars
width = 0.1 # width of error bars
) +
geom_line(aes(group=eco)) +
geom_point() +
labs(y = "Transpiration")+
scale_color_manual(values=c("deepskyblue3", "darkgoldenrod"))
Boxplots
library(ggplot2)
ggplot(data=yp_data, aes(x=temp,y=Trmmol,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("Transp"))+
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(Trmmol ~ temp + eco, data = yp_data)
interaction <- aov(Trmmol ~ 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 79.43 0.00 0.91 0.91 -33.05
interaction 7 83.99 4.56 0.09 1.00 -31.49-> “No interaction” model is the best fit.
- Run ANOVA model
nointeract <- aov(Trmmol ~ temp + eco, data = yp_data)
summary(nointeract) Df Sum Sq Mean Sq F value Pr(>F)
temp 2 79.46 39.73 36.01 2.35e-07 ***
eco 1 23.31 23.31 21.12 0.000175 ***
Residuals 20 22.07 1.10
---
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(nointeract) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = Trmmol ~ temp + eco, data = yp_data)
$temp
diff lwr upr p adj
27.5°C-20°C 3.8075 2.47874 5.13626 0.0000015
35°C-20°C 3.9100 2.58124 5.23876 0.0000010
35°C-27.5°C 0.1025 -1.22626 1.43126 0.9792432
$eco
diff lwr upr p adj
LR-HR 1.970833 1.076313 2.865354 0.000175- The difference between 20 and 27.5, and between 20 and 35deg, is stat significant. The effect of ecotype is confirmed.
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(Trmmol ~ treatment, yp_data)
anova(yp_model)Analysis of Variance Table
Response: Trmmol
Df Sum Sq Mean Sq F value Pr(>F)
treatment 5 105.438 21.0876 19.577 1.036e-06 ***
Residuals 18 19.389 1.0772
---
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(Trmmol ~ 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 105.44 21.09 19.577 1.04e-06 ***
treatment: Eco 1 23.31 23.31 21.635 0.000199 ***
treatment: temp lin 1 61.15 61.15 56.771 5.69e-07 ***
treatment: temp quad 1 18.30 18.30 16.991 0.000640 ***
treatment: e*t lin 1 0.96 0.96 0.892 0.357550
treatment: e*t quad 1 1.72 1.72 1.595 0.222789
Residuals 18 19.39 1.08
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant response of transp to ecotype.
The response of transp to temperature may have a quadratic component.