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)Plant hydraulic conductance
Plant hydraulic conductance analysis
- Prepare data
- Visualize
Plant Hydraulic Conductance 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 = "phc",
groupvars = c("temp", "eco")
)
library(ggplot2)
ggplot(
data = summary_stat,
aes(x = temp, y = phc, colour = eco)
) +
geom_errorbar(aes(ymin = phc - se, ymax = phc + se), # add error bars
width = 0.1 # width of error bars
) +
geom_line(aes(group=eco)) +
geom_point() +
labs(y = "phc")+
scale_color_manual(values=c("deepskyblue3", "darkgoldenrod"))
Boxplots
library(ggplot2)
ggplot(data=yp_data, aes(x=temp,y=phc,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("phc"))+
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: phc ~ temp + eco
interaction: phc ~ temp * eco
-> Identify best model using AIC
nointeract <- aov(phc ~ temp + eco, data = yp_data)
interaction <- aov(phc ~ 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 89.09 0.0 1 1 -34.04
no interaction 5 123.19 34.1 0 1 -54.93-> “Interaction” model is the best fit.
- Run ANOVA model
interaction <- aov(phc ~ temp * eco, data = yp_data)
summary(interaction) Df Sum Sq Mean Sq F value Pr(>F)
temp 2 28.84 14.42 10.82 0.000819 ***
eco 1 261.86 261.86 196.57 3.97e-11 ***
temp:eco 2 112.70 56.35 42.30 1.57e-07 ***
Residuals 18 23.98 1.33
---
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 = phc ~ temp * eco, data = yp_data)
$temp
diff lwr upr p adj
27.5°C-20°C -2.61775 -4.0905838 -1.1449162 0.0007126
35°C-20°C -1.82650 -3.2993338 -0.3536662 0.0141079
35°C-27.5°C 0.79125 -0.6815838 2.2640838 0.3762875
$eco
diff lwr upr p adj
LR-HR 6.606333 5.616392 7.596275 0
$`temp:eco`
diff lwr upr p adj
27.5°C:HR-20°C:HR 1.89575 -0.6979455 4.489445 0.2354583
35°C:HR-20°C:HR 2.84925 0.2555545 5.442945 0.0266135
20°C:LR-20°C:HR 12.73250 10.1388045 15.326195 0.0000000
27.5°C:LR-20°C:HR 5.60125 3.0075545 8.194945 0.0000261
35°C:LR-20°C:HR 6.23025 3.6365545 8.823945 0.0000062
35°C:HR-27.5°C:HR 0.95350 -1.6401955 3.547195 0.8458063
20°C:LR-27.5°C:HR 10.83675 8.2430545 13.430445 0.0000000
27.5°C:LR-27.5°C:HR 3.70550 1.1118045 6.299195 0.0029392
35°C:LR-27.5°C:HR 4.33450 1.7408045 6.928195 0.0005808
20°C:LR-35°C:HR 9.88325 7.2895545 12.476945 0.0000000
27.5°C:LR-35°C:HR 2.75200 0.1583045 5.345695 0.0339182
35°C:LR-35°C:HR 3.38100 0.7873045 5.974695 0.0068260
27.5°C:LR-20°C:LR -7.13125 -9.7249455 -4.537555 0.0000009
35°C:LR-20°C:LR -6.50225 -9.0959455 -3.908555 0.0000034
35°C:LR-27.5°C:LR 0.62900 -1.9646955 3.222695 0.9690908- The difference between 20 and 27.5 (both temp only and temp*eco) and between 20 and 35deg (temp only and temp*eco) 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(phc ~ treatment, yp_data)
anova(yp_model)Analysis of Variance Table
Response: phc
Df Sum Sq Mean Sq F value Pr(>F)
treatment 5 403.40 80.679 60.564 1.257e-10 ***
Residuals 18 23.98 1.332
---
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 phc is the same for both ecotypes |
| 2: Contrast ‘temp linear’ (1,0,-1,1,0,-1) -> H0 = Mean phc for 20deg = mean phc for 35deg |
| 3: Contrast ‘temp quadratic’ (1,-2,1,1,-2,1) ->H0 = Mean phc at 27.5deg = Mean phc (20deg,35deg) |
| 4: Contrast ‘ecotype * temp linear’ (1,0,-1,-1,0,1) -> H0 = Linear component of phc ~ temp is the same for both ecotypes |
| 5: Contrast ‘ecotype * temp quadratic’ (1,-2,1,-1,2,-1) -> H0 = Quadratic component of phc ~ 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(phc ~ 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 403.4 80.68 60.56 1.26e-10 ***
treatment: Eco 1 261.9 261.86 196.57 3.97e-11 ***
treatment: temp lin 1 13.3 13.34 10.02 0.005359 **
treatment: temp quad 1 15.5 15.50 11.63 0.003118 **
treatment: e*t lin 1 87.5 87.45 65.65 2.04e-07 ***
treatment: e*t quad 1 25.2 25.24 18.95 0.000383 ***
Residuals 18 24.0 1.33
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant response of phc to ecotype.
The response of phc to temperature may have a quadratic component that depends on ecotype.