Growth rates in A. cervicornis under elevated N and N+P and subsequent heat stress
Ana Palacio
Feb 13, 2020
General project set-up
# Load all libraries and sources required to run the script
library(tidyverse)
library(ggthemes)
library(lmerTest)
library(emmeans)
library(multcomp)
library(gridExtra)
library(rstatix)
# Default ggplot settings
Fill.colour<-scale_colour_manual(values = c("black", "gray70", "gray35"))
ggthe_bw<-theme(plot.background=element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)+
theme_bw()Import and clean growth data
# Import data:
## Growth long
BW.Tall<-read.csv("Growth_corrected.csv", header = TRUE)
## Genotype information
MOTE<-read.csv("Genotype_Info.csv", header = TRUE)
BW.Tall<-merge(BW.Tall, MOTE, by="Gen", all.x=TRUE)
## Data type
BW.Tall$Initial.Date<-as.Date(BW.Tall$Initial.Date, "%m/%d/%y")
BW.Tall$Final.Date<-as.Date(BW.Tall$Final.Date, "%m/%d/%y")
# BW.Tall$Initial.Date<-as.Date(BW.Tall$Initial.Date, "%Y-%m-%d")
# BW.Tall$Final.Date<-as.Date(BW.Tall$Final.Date, "%Y-%m-%d")
BW.Tall$Days<-(as.numeric(as.Date(BW.Tall$Final.Date))-17485)
BW.Tall$DayF<-(as.factor(BW.Tall$Days))
# BW.Tall$TimePoint[BW.Tall$Day=="-28"]<-"3"
# BW.Tall$TimePoint[BW.Tall$Day=="-7"]<-"4"
# BW.Tall$TimePoint[BW.Tall$Day=="28"]<-"9"
# BW.Tall$TimePoint[BW.Tall$Day=="62"]<-"11"
# BW.Tall$TimePoint[BW.Tall$Day=="75"]<-"13"
# BW.Tall$TimePoint[BW.Tall$Day=="91"]<-"16"
# BW.Tall$TimePoint[BW.Tall$Day=="100"]<-"17"
# BW.Tall$TimePoint[BW.Tall$Day=="113"]<-"22"
BW.Tall$TimePoint[BW.Tall$DayF=="-28"]<-"1"
BW.Tall$TimePoint[BW.Tall$DayF=="-7"]<-"NA"
BW.Tall$TimePoint[BW.Tall$DayF=="28"]<-"2"
BW.Tall$TimePoint[BW.Tall$DayF=="62"]<-"3"
BW.Tall$TimePoint[BW.Tall$DayF=="75"]<-"4"
BW.Tall$TimePoint[BW.Tall$DayF=="91"]<-"5"
BW.Tall$TimePoint[BW.Tall$DayF=="100"]<-"6"
BW.Tall$TimePoint[BW.Tall$DayF=="113"]<-"7"
BW.Tall$TimePoint<-as.numeric(as.character(BW.Tall$TimePoint))
BW.Tall$Treatment <- as.character(BW.Tall$Treatment)
BW.Tall$Treatment[BW.Tall$Treatment == "Control"] <- "A"
BW.Tall$Treatment[BW.Tall$Treatment == "NP"] <- "N+P"
BW.Tall$Treatment <- as.factor(BW.Tall$Treatment)
BW.Tall$Nutrients <- "Nutrients"
BW.Tall$Nutrients[BW.Tall$Treatment == "A"] <- "Ambient"
BW.Tall$Colony <- as.factor(BW.Tall$Gen)
BW.Tall$MoteGen<-factor(BW.Tall$MoteGen,
levels=c("G_48", "G_62","G_31", "G_08","G_07", "G_50"))
BW.Tall$SampleID <- paste("Ac", BW.Tall$Fra, BW.Tall$TimePoint, sep = "_")
## Data subsets
BW.Initial<-subset(BW.Tall, Initial.Date=="2017-08-30") # Only baseline
BW.Initial<-droplevels(BW.Initial)
BW.Control<-subset(BW.Tall, Treatment=="A") # Only Ambient
BW.Control<-subset(BW.Control, Days<80)
BW.Control<-subset(BW.Control, Initial.Date!="2017-10-18")
BW.Control<-droplevels(BW.Control)
BW.FinalAC<-subset(BW.Control, DayF=="75")
BW.FinalAC<-droplevels(BW.FinalAC)
#BW.nutrients<-subset(BW.Tall, Days>-10) # Remove baseline
#BW.nutrients<-subset(BW.nutrients, Days<76)
BW.nutrients<-subset(BW.Tall, Days<76) # with baseline
BW.nutrients<-droplevels(BW.nutrients)
BW.100<-subset(BW.Tall, DayF=="100")
BW.91<-subset(BW.Tall, DayF=="91")
## Remove recovery data points
BW.Tall<-subset(BW.Tall, Initial.Date!="2017-10-18")
BW.Tall<-subset(BW.Tall, Initial.Date!="2018-02-23") # After heat stress
BW.selected<-subset(BW.Tall, Days!="28")
BW.selected<-subset(BW.selected, Days!="62")
N.fragments<-BW.Tall %>%
group_by(MoteGen) %>% count(DayF)
N.fragmentsGrowth over time (during Control and Heat stress) by nutrient treatment
Figure 3a
# Treatment effect over time
BW_Treat<- ggplot(BW.Tall, aes (Days, dAW.gd, shape=factor(Treatment),
colour=factor(Treatment))) +
ggthe_bw + Fill.colour + ggtitle("a")+
stat_summary(fun.data = "mean_cl_boot", geom = "errorbar",
position=position_dodge(width=5))+
stat_summary(fun.y=mean, geom="line",
position=position_dodge(width=5), linetype=2) +
stat_summary(fun.y=mean, geom="point", size =2,
position=position_dodge(width=5)) +
geom_hline(yintercept = 0, linetype=3)+
scale_x_continuous(name="Days in the experiment",
breaks = seq(0, 115, by=15)) +
scale_y_continuous(name="Growth (mg / g d)",
breaks = seq(-1, 4, by=1)) +
theme(legend.position=c(0.2, 0.5),
legend.title = element_blank(),
strip.background =element_rect(fill=NA)) +
annotate("segment", x = -30, xend = -2, y = -1.5, yend = -1.5,
colour = "gray90")+
annotate("text", x = -16, y = -1.2, label = "Baseline", size=3)+
annotate("segment", x = 2, xend = 91, y = -1.5, yend = -1.5,
colour = "gray90")+
annotate("text", x = 45, y = -1.2, label = "Nutrients", size=3)+
annotate("segment", x = 79, xend = 91, y = -1.4, yend = -1,
colour = "gray90")+
annotate("segment", x = 91, xend = 110, y = -1, yend = -1,
colour = "gray90")+
annotate("text", x = 99, y = -1.2, label = "Heat", size=3)
BW_Treat
# ggsave(file="Outputs/Figure3_BW_corrected.svg", plot=BW_Treat, dpi = 300, units="in", width=3.5, height=3.0)Figure S2a: All genotypes by nutrient treatment over time
Figure S2a: All genotypes by nutrient treatment over time (differrence)
# Claculate relative changes (1 - Control)
BW.Tall$Nutrients2<-BW.Tall$Nutrients
BW.Tall$Nutrients2[BW.Tall$DayF=="-28"]<-"Ambient"
BW.Summary <-BW.Tall %>%
group_by(MoteGen, Nutrients2, DayF) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
BW.Summary BW.Summary<-subset(BW.Summary, Nutrients2=="Ambient" )
BW.Summary<-BW.Summary[, -3]
BW.Summary<-BW.Summary[, -3]
BW.Summary<-BW.Summary[, -3]
BW.Tall<-merge(BW.Tall, BW.Summary, by=c("MoteGen", "DayF"), all.x = T)
BW.Tall$Difference<-BW.Tall$dAW.gd-BW.Tall$mean
Colour.colour<-scale_colour_manual(values = c("black", "gray70"))
Fill.colour<-scale_fill_manual(values = c("black", "gray70"))
BW_Diff_Col <- ggplot(BW.Tall, aes(Days, Difference, fill=Nutrients2,
shape=factor(Nutrients2))) +
#annotate("segment", x = -30, xend = -2, y = -1.5, yend = -1.5,
# colour = "gray90")+
# annotate("text", x = -16, y = -1.2, label = "BL", size=3)+
annotate("segment", x = 2, xend = 91, y = -3.5, yend = -3.5,
colour = "gray90")+
annotate("text", x = 45, y = -3.9, label = "Nutrients", size=4, colour="gray")+
annotate("segment", x = 79, xend = 91, y = -3.5, yend = 1,
colour = "gray90")+
annotate("segment", x = 91, xend = 110, y =1, yend =1,
colour = "gray90")+
annotate("text", x = 99, y = 0.8, label = "H", size=4, colour="gray")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position=position_dodge(width=5), width = 10)+
stat_summary(fun.y=mean, geom="point", size =2,
position=position_dodge(width=5), alpha=0.8) +
stat_summary(fun.y=mean, geom="line",
position=position_dodge(width=5), linetype=1) +
scale_shape_manual(values=c(21, 14),
labels=c("A", "N and N+P"))+
ggthe_bw + Colour.colour + Fill.colour +
scale_x_continuous(name="Days in the experiment",
limits = c(-35,113),
breaks = seq(0, 113, 30),
expand = c(0, 0))+
scale_y_continuous(name="Difference in growth (mg/ g d) respecto to Ambient",
breaks = seq(-4, 1.5, by=1)) +
geom_hline(yintercept = 0, linetype=3)+
theme(text=element_text(size=14),
legend.position=c(0.1, 0.2),
legend.title = element_blank(),
strip.background =element_rect(fill=NA)) +
facet_wrap(~MoteGen)
BW_Diff_Col
#ggsave(file="Outputs/S_3_BW_Diff_Col.svg", plot=BW_Diff_Col, width=5.0, height=5)- Percantages
# Claculate percentajes changes (1 - Control)
BW.Tall$percentage<-(BW.Tall$Difference/BW.Tall$mean)*100
Colour.colour<-scale_colour_manual(values = c("black", "gray70"))
Fill.colour<-scale_fill_manual(values = c("black", "gray70"))
BW_perent_Col <- ggplot(BW.Tall, aes(Days, percentage, fill=Nutrients2,
shape=factor(Nutrients2))) +
#annotate("segment", x = -30, xend = -2, y = -1.5, yend = -1.5,
# colour = "gray90")+
# annotate("text", x = -16, y = -1.2, label = "BL", size=3)+
annotate("segment", x = 2, xend = 91, y = -140, yend = -140,
colour = "gray90")+
annotate("text", x = 45, y = -135, label = "Nutrients", size=4, colour="gray")+
annotate("segment", x = 79, xend = 91, y = -139, yend = 28,
colour = "gray90")+
annotate("segment", x = 91, xend = 110, y =28, yend =28,
colour = "gray90")+
annotate("text", x = 99, y = 33, label = "H", size=4, colour="gray")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
position=position_dodge(width=5), width = 10)+
stat_summary(fun.y=mean, geom="point", size =2,
position=position_dodge(width=5), alpha=0.8) +
stat_summary(fun.y=mean, geom="line",
position=position_dodge(width=5), linetype=1) +
scale_shape_manual(values=c(21, 14),
labels=c("A", "N and N+P"))+
ggthe_bw + Colour.colour + Fill.colour +
scale_x_continuous(name="Days in the experiment",
limits = c(-35,113),
breaks = seq(0, 113, 30),
expand = c(0, 0))+
scale_y_continuous(name="Percentage of growth change respecto to Ambient",
breaks = seq(-140, 100, by=20)) +
geom_hline(yintercept = 0, linetype=3)+
theme(text=element_text(size=14),
legend.position=c(0.1, 0.2),
legend.title = element_blank(),
strip.background =element_rect(fill=NA)) +
facet_wrap(~MoteGen)
BW_perent_Col
#ggsave(file="Outputs/S_3_BW_Percentaje_Col.svg", plot=BW_perent_Col, width=5.0, height=5)Figure 3b: By genotype, but N and N+P pooled
# sample means
BW.nutrients1 <- aggregate(dAW.gd~Nutrients+Days+MoteGen, data=BW.Tall, mean, na.rm=TRUE)
names(BW.nutrients1)[4] <- "mean"
# sample SDs
BW.nutrientssd1 <- aggregate(dAW.gd~Nutrients+Days+MoteGen, data=BW.Tall, sd, na.rm=TRUE)
names(BW.nutrientssd1)[4] <- "sd"
# Collect the summary statistics in a dataframe
BW.nutrients1<-merge(BW.nutrients1, BW.nutrientssd1)
BW.nutrients1<-BW.nutrients1[order(BW.nutrients1$Days, BW.nutrients1$MoteGen),]
BW.nutrients1All genotypes (colour) by nutrient treatment over time
Other interesting figures
Compare individual genotypes inside each nutrient treatment * Growth declined in G50, G07 and G08 in A, but increased in G48, G62 and G31
BW_ColonyT<- ggplot(BW.Tall, aes (Days, dAW.gd, colour=MoteGen)) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.8)+
stat_summary(fun.y=mean, geom="point", size =2, alpha=0.5) +
ylab("Growth (mg/ g d)") +
annotate("segment", x = 2, xend = 91, y = -1, yend = -1,
colour = "gray90", linetype=1)+
annotate("segment", x = 79, xend = 91, y = -1, yend = -0.5,
colour = "gray90", linetype=1)+
annotate("segment", x = 91, xend = 110, y = -0.5, yend = -0.5,
colour = "gray90", linetype=1)+
annotate("text", x = 45, y = -1.5, label = "Nutrients", size=3)+
annotate("text", x = 99, y = -1.5, label = "H", size=3)+
stat_summary(fun.y=mean, geom="line") + ggthe_bw + facet_grid (~Treatment)
BW_ColonyT
Compare individual fragments inside each nutrient treatment
BW_ColonyTF<- ggplot(BW.Tall, aes (Days, dAW.gd, colour=factor(Fra), group=Fra)) + stat_summary(fun.y=mean, geom="line") + ggthe_bw + facet_grid (MoteGen~Treatment)+ geom_point()+ theme(legend.position = "none")
BW_ColonyTF
Initial performance vs ambient per genotype: some genotypes declined in growth even without nutrients
BW_Init<- ggplot(BW.Initial, aes (MoteGen, dAW.gd, colour=MoteGen)) +
ggtitle("(a) Baseline")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2) +
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) + ggthe_bw +
scale_y_continuous(name=("Growth (mg/ g d) "),
limits = c(0,6),
breaks = seq(0, 6, by=0.5))+
theme(legend.position=c(0.3, 0.3),
legend.title = element_blank(),
strip.background =element_rect(fill=NA))
BW_Final<- ggplot(BW.FinalAC, aes (MoteGen, dAW.gd, colour=MoteGen)) +
ggtitle("(a) Ambient (day 75)")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2) +
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) + ggthe_bw +
scale_y_continuous(name=("Growth (mg/ g d) "),
limits = c(0,6),
breaks = seq(0, 6, by=0.5))+
ylab("Growth (mg/ g d) ")+
theme(legend.position="none",
legend.title = element_blank(),
strip.background =element_rect(fill=NA))
BW_Colony_Control<- ggplot(BW.Control, aes (Days, dAW.gd, colour=MoteGen)) +
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar",
width = 10, position=position_dodge(width=5)) +
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5,
position=position_dodge(width=5)) +
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=5)) +
ggthe_bw + # geom_jitter()+
ggtitle("(b) Ambient - Control") +
scale_y_continuous(name=("Growth (mg/ g d) "),
limits = c(0,6),
breaks = seq(0, 6, by=0.5))+
theme(legend.position="none")
Captivityeffect<-grid.arrange(BW_Init,BW_Final, nrow=1)
Captivityeffect2<-grid.arrange(BW_Init,BW_Colony_Control, nrow=1)
#ggsave(file="Outputs/S_Genotypes.svg", plot=Captivityeffect, width=6, height=3.5)BW_Colony_Control<-ggplot(BW.Control, aes (Days, dAW.gd, colour=factor(Fra), group=Fra)) +
#stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.8)+
#stat_summary(fun.y=mean, geom="point", size =2, alpha=0.5) +
stat_summary(fun.y=mean, geom="line") + ggthe_bw +
facet_wrap (~MoteGen)+ geom_point()+ theme(legend.position = "none")
BW_Colony_Control
LMER model selection and pairwise comparisons
Treatment
LMER4
# Complex to simple models
LM_1 <- lmer(dAW.gd ~ Treatment*DayF +
(1 + Treatment|Genotype) + (1|Replicate), REML=TRUE, data= BW.Tall)
LM_2 <- lmer(dAW.gd ~ Treatment * DayF +
(Treatment|Genotype), REML=TRUE, data= BW.Tall)
LM_3 <- lmer(dAW.gd ~ Treatment * DayF +
(1|Genotype), REML=TRUE, data= BW.Tall)
# Select model
AIC(LM_1, LM_2, LM_3) # 1=2, but 3 has lower AIC anova(LM_1, LM_2, refit=FALSE) anova(LM_2, LM_3, refit=FALSE) anova(LM_1, LM_3, refit=FALSE)# Final model
LM_Treatment_Days<-lmer(dAW.gd ~ Treatment * DayF +
(1|Genotype) + (1|Fra), data= BW.Tall)
# Removed (1|Replicate) since it was not significant (p=0.1642)
isSingular(LM_Treatment_Days)## [1] FALSE anova(LM_Treatment_Days) ranova(LM_Treatment_Days) #summary(LM_Treatment_Days)
par(mfrow=c(1,2))
plot(fitted(LM_Treatment_Days), residuals(LM_Treatment_Days, type="pearson"))
abline(h=0)
qqnorm(residuals(LM_Treatment_Days, type="pearson"))
abline(0,1)
# Pairwise comparisons
BW.emmc<-emmeans(LM_Treatment_Days, ~Treatment*DayF)
#contrast(BW.emmc, "tukey")
BW_groups<-cld(BW.emmc)
#BW.emmc<-emmeans(LM_Treatment_Days, ~Treatment|DayF)
#contrast(BW.emmc, "tukey")
#BW_groups<-cld(BW.emmc)
#BW_groups<-as.data.frame(BW_groups[complete.cases(BW_groups),])
BW_groups<-BW_groups[order(BW_groups$DayF,BW_groups$Treatment),]
BW_groups#write.csv(BW_groups, "Outputs/BW_groups_corrected.csv", row.names = F)nlme
#library(nlme)
#
# fit1 <- gls(dAW.gd ~ DayF *Treatment, data=BW.Tall,
# correlation=corSymm(form=~TimePoint|Fra),
# weights=varIdent(form=~1|TimePoint),
# na.action=na.exclude,
# control=glsControl(opt="optim"))Only Ambient summary model
# summary(fit1)
#
# intervals(fit1)
# anova(fit1, type="marginal")
# getVarCov(fit1)
# residuals(fit1, type="pearson")
#
# par(mfrow=c(1,2))
# plot(fitted(fit1), residuals(fit1, type="pearson"))
# abline(h=0)
# qqnorm(residuals(fit1, type="pearson"))
# abline(0,1)Predictions
# library(rms)
# clda_Gls_1 <- Gls(dAW.gd ~ DayF * Treatment,
# correlation=corSymm (form = ~TimePoint|Fra),
# weights=varIdent(form=~1|TimePoint),
# na.action=na.exclude,
# control=glsControl(opt="optim"),
# data=BW.Tall)
#
# data1<-expand.grid(DayF=unique(BW.Tall$DayF),
# Treatment=unique(BW.Tall$Treatment))
#
# predictions_1 <- cbind(data1,
# predict (clda_Gls_1, data1, conf.int=0.95))
#
# predictions_1$Days<-as.numeric(as.character(predictions_1$DayF))Model prediction by treatment
# pd <- position_dodge(2)
# limits <- aes(ymax = upper , ymin=lower, shape=Treatment)
#
# pCI1 <- ggplot(predictions_1, aes(y=linear.predictors, x=Days,
# fill=Treatment)) +
# geom_errorbar(limits, width= 2 , position=pd) +
# geom_line(aes(group=Treatment, y=linear.predictors), position=pd, size=0.1) +
# geom_point(aes(fill=Treatment), shape=21, position=pd, size=4, alpha=0.8) +
#
# ggthe_bw+
#
# geom_hline(yintercept = 0, linetype=3)+
#
# scale_x_continuous(name="Days in the experiment",
# breaks = seq(0, 115, by=15)) +
# scale_y_continuous(name="Estimated growth (mg / g d) ",
# breaks = seq(-1, 4, by=1)) +
#
# theme(text=element_text(size=14),
# legend.position="right",
# legend.title = element_blank(),
# strip.background =element_rect(fill=NA)) +
#
# annotate("segment", x = -30, xend = -2, y = -1.5, yend = -1.5,
# colour = "gray90")+
# annotate("text", x = -16, y = -1.2, label = "Baseline", size=3)+
#
# annotate("segment", x = 2, xend = 91, y = -1.5, yend = -1.5,
# colour = "gray90")+
# annotate("text", x = 45, y = -1.2, label = "Nutrients", size=3)+
#
# annotate("segment", x = 79, xend = 91, y = -1.4, yend = -1,
# colour = "gray90")+
#
# annotate("segment", x = 91, xend = 110, y = -1, yend = -1,
# colour = "gray90")+
# annotate("text", x = 99, y = -1.2, label = "Heat", size=3)
#
# pCI1Genets
No nutrients
lmer4
- Baseline differences
# ANOVA
LM_bl <- lmer(dAW.gd ~ MoteGen + (1|Replicate), data= BW.Initial)
anova(LM_bl) ranova(LM_bl) LM_bl <- lm(dAW.gd ~ MoteGen, data= BW.Initial)
anova(LM_bl) summary(LM_bl)##
## Call:
## lm(formula = dAW.gd ~ MoteGen, data = BW.Initial)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.15347 -0.45620 0.03423 0.49835 2.66945
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.12658 0.16193 19.309 < 2e-16 ***
## MoteGenG_62 0.03993 0.22702 0.176 0.860704
## MoteGenG_31 1.16789 0.26853 4.349 3.03e-05 ***
## MoteGenG_08 0.94450 0.38546 2.450 0.015821 *
## MoteGenG_07 0.79777 0.23336 3.419 0.000879 ***
## MoteGenG_50 1.10776 0.28757 3.852 0.000196 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8568 on 112 degrees of freedom
## Multiple R-squared: 0.2577, Adjusted R-squared: 0.2246
## F-statistic: 7.776 on 5 and 112 DF, p-value: 2.599e-06# Pairwise comparisons
BW_bl.emmc<-emmeans(LM_bl, ~MoteGen)
#contrast(BW_bl.emmc, "tukey")
#Tukey groups
BW_bl_groups<-cld(BW_bl.emmc)
#write.csv(BW_bl_groups, "Outputs/BW_baseline_groups.csv", row.names = F)
# Summary stats
# library(rstatix)
Stats_BL<-BW.Initial %>%
group_by(MoteGen) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
Stats_BL- Ambient and Control evolution
summary(BW.Control)## Gen Fra Treatment Replicate Time
## G-O :35 Min. :102.0 A:154 R1:83 1.Sep-Oct:39
## Green :40 1st Qu.:149.0 R2:71 3.Nov-Dec:38
## Orange:20 Median :182.0 4.Dec-Jan:38
## R-O :35 Mean :194.1 5.Jan-Feb:39
## R-Y :16 3rd Qu.:257.8
## Yellow: 8 Max. :290.0
## Period Initial.Date Final.Date dBW.g
## 1.Sep-Oct:39 Min. :2017-08-30 Min. :2017-10-18 Min. :1.242
## 4.Dec-Jan:38 1st Qu.:2017-09-18 1st Qu.:2017-11-01 1st Qu.:2.898
## Jan-Feb :39 Median :2017-11-30 Median :2017-12-30 Median :3.668
## Nov-Dec :38 Mean :2017-11-18 Mean :2017-12-19 Mean :3.568
## 3rd Qu.:2018-01-08 3rd Qu.:2018-01-25 3rd Qu.:4.197
## Max. :2018-01-16 Max. :2018-01-29 Max. :6.337
## dAW.gd Original.Collection.Date MoteGen Collection.Habitat
## Min. :1.252 02/23/08:43 G_48:40 Midchannel:94
## 1st Qu.:2.935 03/18/10:20 G_62:35 Offshore :60
## Median :3.675 05/21/13:35 G_31:20
## Mean :3.578 05/28/10:40 G_08: 8
## 3rd Qu.:4.213 07/30/10:16 G_07:35
## Max. :6.375 G_50:16
## Lat Lon Genotype WB
## Min. :24.52 Min. :-81.64 Green :40 R: 20
## 1st Qu.:24.52 1st Qu.:-81.58 Green and Orange:35 S:134
## Median :24.55 Median :-81.53 Orange :20
## Mean :24.54 Mean :-81.55 Red and Orange :35
## 3rd Qu.:24.56 3rd Qu.:-81.50 Red and Yellow :16
## Max. :24.58 Max. :-81.43 Yellow : 8
## Days DayF TimePoint Nutrients Colony
## Min. :-28.00 -28:39 Min. :1.00 Length:154 G-O :35
## 1st Qu.:-14.00 28 :38 1st Qu.:1.25 Class :character Green :40
## Median : 45.00 62 :38 Median :2.50 Mode :character Orange:20
## Mean : 34.11 75 :39 Mean :2.50 R-O :35
## 3rd Qu.: 71.75 3rd Qu.:3.75 R-Y :16
## Max. : 75.00 Max. :4.00 Yellow: 8
## SampleID
## Length:154
## Class :character
## Mode :character
##
##
## # Model
LM_A_C.0 <- lmer(dAW.gd ~ MoteGen * DayF + (1|Replicate), data= BW.Control)
anova(LM_A_C.0) # Day alone is not significant, but Day:Genet is ranova(LM_A_C.0) # Replicate is not significant LM_A_C <-lmer(dAW.gd ~ MoteGen * DayF + (1|Fra), data= BW.Control)
anova(LM_A_C) # Day alone is not significant, but Day:Genet is ranova(LM_A_C)# Pairwise comparisons
BW_A_C.emmc<-emmeans(LM_A_C, ~MoteGen | DayF)
#contrast(BW_A_C.emmc, "tukey")
#BW_A_C.emmc<-emmeans(LM_A_C, ~DayF|MoteGen)
#contrast(BW_A_C.emmc, "tukey")
#Tukey groups
BW_A_C_groups<-cld(BW_A_C.emmc)
BW_A_C_groups<-BW_A_C_groups[order(BW_A_C_groups$MoteGen,BW_A_C_groups$Day),]
BW_A_C_groups #write.csv(BW_A_C_groups, "Outputs/BW_A_C_groups.csv", row.names = F)
# Summary stats
# library(rstatix)
Stats_Control<-BW.Control %>%
group_by(MoteGen) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
Stats_Control- Ambient final day
summary(BW.FinalAC)## Gen Fra Treatment Replicate Time
## G-O : 9 Min. :102.0 A:39 R1:21 5.Jan-Feb:39
## Green :10 1st Qu.:150.5 R2:18
## Orange: 5 Median :182.0
## R-O : 9 Mean :194.7
## R-Y : 4 3rd Qu.:254.5
## Yellow: 2 Max. :290.0
## Period Initial.Date Final.Date dBW.g
## Jan-Feb:39 Min. :2018-01-16 Min. :2018-01-29 Min. :1.601
## 1st Qu.:2018-01-16 1st Qu.:2018-01-29 1st Qu.:2.476
## Median :2018-01-16 Median :2018-01-29 Median :3.655
## Mean :2018-01-16 Mean :2018-01-29 Mean :3.469
## 3rd Qu.:2018-01-16 3rd Qu.:2018-01-29 3rd Qu.:4.241
## Max. :2018-01-16 Max. :2018-01-29 Max. :5.071
## dAW.gd Original.Collection.Date MoteGen Collection.Habitat
## Min. :1.632 02/23/08:11 G_48:10 Midchannel:24
## 1st Qu.:2.501 03/18/10: 5 G_62: 9 Offshore :15
## Median :3.580 05/21/13: 9 G_31: 5
## Mean :3.436 05/28/10:10 G_08: 2
## 3rd Qu.:4.220 07/30/10: 4 G_07: 9
## Max. :4.996 G_50: 4
## Lat Lon Genotype WB
## Min. :24.52 Min. :-81.64 Green :10 R: 5
## 1st Qu.:24.53 1st Qu.:-81.58 Green and Orange: 9 S:34
## Median :24.55 Median :-81.53 Orange : 5
## Mean :24.54 Mean :-81.55 Red and Orange : 9
## 3rd Qu.:24.56 3rd Qu.:-81.50 Red and Yellow : 4
## Max. :24.58 Max. :-81.43 Yellow : 2
## Days DayF TimePoint Nutrients Colony
## Min. :75 75:39 Min. :4 Length:39 G-O : 9
## 1st Qu.:75 1st Qu.:4 Class :character Green :10
## Median :75 Median :4 Mode :character Orange: 5
## Mean :75 Mean :4 R-O : 9
## 3rd Qu.:75 3rd Qu.:4 R-Y : 4
## Max. :75 Max. :4 Yellow: 2
## SampleID
## Length:39
## Class :character
## Mode :character
##
##
## # Model
LM_A_C.75 <- lmer(dAW.gd ~ MoteGen + (1|Replicate), data= BW.FinalAC)
isSingular(LM_A_C.75)## [1] TRUE anova(LM_A_C.75) # Day alone is not significant, but Day:Genet is ranova(LM_A_C.75) # Replicate is not significant LM_A_C.75 <-lm(dAW.gd ~ MoteGen , data=BW.FinalAC)
anova(LM_A_C.75) # Day alone is not significant, but Day:Genet is summary(LM_A_C.75)##
## Call:
## lm(formula = dAW.gd ~ MoteGen, data = BW.FinalAC)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14549 -0.25813 0.04602 0.39836 0.86661
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.13512 0.21260 19.450 < 2e-16 ***
## MoteGenG_62 -0.27343 0.30890 -0.885 0.382465
## MoteGenG_31 -0.03574 0.36823 -0.097 0.923267
## MoteGenG_08 -2.07257 0.52076 -3.980 0.000356 ***
## MoteGenG_07 -1.64485 0.30890 -5.325 7.07e-06 ***
## MoteGenG_50 -1.41577 0.39774 -3.560 0.001152 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6723 on 33 degrees of freedom
## Multiple R-squared: 0.6024, Adjusted R-squared: 0.5421
## F-statistic: 9.998 on 5 and 33 DF, p-value: 6.847e-06# Pairwise comparisons
BW_A_C75.emmc<-emmeans(LM_A_C.75, ~MoteGen)
#contrast(BW_A_C.emmc, "tukey")
#Tukey groups
BW_A_C75_groups<-cld(BW_A_C75.emmc)
BW_A_C75_groups #write.csv(BW_A_C75_groups, "Outputs/BW_A_C75_groups.csv", row.names = F)
Stats_75_a<-BW.FinalAC %>%
group_by(MoteGen) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
Stats_75_a- Baseline + Ambient-Control overtime
BW.AC<-subset(BW.Tall, Days<80)
BW.AC$Treatment[BW.AC$DayF=="-28"]<-"A"
BW.AC<-subset(BW.AC, Treatment=="A") # Only Ambient
BW.AC<-droplevels(BW.AC)
Stats_Con_BL<-BW.AC %>%
group_by(MoteGen, DayF) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
Stats_Con_BL summary(BW.AC)## MoteGen DayF Gen Fra Treatment Replicate
## G_48:58 -28:118 G-O :55 Min. :101.0 A:233 R1:123
## G_62:55 28 : 38 Green :58 1st Qu.:152.0 R2:110
## G_31:31 62 : 38 Orange:31 Median :186.0
## G_08:12 75 : 39 R-O :52 Mean :196.6
## G_07:52 R-Y :25 3rd Qu.:261.0
## G_50:25 Yellow:12 Max. :300.0
## Time Period Initial.Date Final.Date
## 1.Sep-Oct:118 1.Sep-Oct:118 Min. :2017-08-30 Min. :2017-10-18
## 3.Nov-Dec: 38 4.Dec-Jan: 38 1st Qu.:2017-08-30 1st Qu.:2017-10-18
## 4.Dec-Jan: 38 Jan-Feb : 39 Median :2017-08-30 Median :2017-10-18
## 5.Jan-Feb: 39 Nov-Dec : 38 Mean :2017-10-22 Mean :2017-11-28
## 3rd Qu.:2017-12-15 3rd Qu.:2018-01-16
## Max. :2018-01-16 Max. :2018-01-29
## dBW.g dAW.gd Original.Collection.Date
## Min. :1.242 Min. :1.252 02/23/08:64
## 1st Qu.:2.891 1st Qu.:2.926 03/18/10:31
## Median :3.648 Median :3.645 05/21/13:55
## Mean :3.602 Mean :3.608 05/28/10:58
## 3rd Qu.:4.232 3rd Qu.:4.242 07/30/10:25
## Max. :6.593 Max. :6.594
## Collection.Habitat Lat Lon
## Midchannel:144 Min. :24.52 Min. :-81.64
## Offshore : 89 1st Qu.:24.53 1st Qu.:-81.58
## Median :24.55 Median :-81.53
## Mean :24.54 Mean :-81.55
## 3rd Qu.:24.56 3rd Qu.:-81.50
## Max. :24.58 Max. :-81.43
## Genotype WB Days TimePoint
## Green :58 R: 31 Min. :-28.00 Min. :1.000
## Green and Orange:55 S:202 1st Qu.:-28.00 1st Qu.:1.000
## Orange :31 Median :-28.00 Median :1.000
## Red and Orange :52 Mean : 13.05 Mean :1.991
## Red and Yellow :25 3rd Qu.: 62.00 3rd Qu.:3.000
## Yellow :12 Max. : 75.00 Max. :4.000
## Nutrients Colony SampleID TreatmentG
## Length:233 G-O :55 Length:233 A:233
## Class :character Green :58 Class :character
## Mode :character Orange:31 Mode :character
## R-O :52
## R-Y :25
## Yellow:12
## Nutrients2 mean sd Difference
## Length:233 Min. :1.975 Min. :0.1070 Min. :-2.5847000
## Class :character 1st Qu.:3.127 1st Qu.:0.6070 1st Qu.:-0.3124000
## Mode :character Median :3.837 Median :0.6330 Median : 0.0464000
## Mean :3.608 Mean :0.7274 Mean : 0.0000013
## 3rd Qu.:4.135 3rd Qu.:1.0220 3rd Qu.: 0.4053000
## Max. :4.793 Max. :1.1730 Max. : 2.6698000
## percentage
## Min. :-67.36252
## 1st Qu.: -9.81814
## Median : 1.17685
## Mean : -0.00074
## 3rd Qu.: 11.05532
## Max. : 68.03772# Model
LM_AC_BL <- lmer(dAW.gd ~ MoteGen * DayF + (1|Replicate), data= BW.AC)
isSingular(LM_AC_BL)## [1] TRUE anova(LM_AC_BL) # Day alone is not significant, but Day:Genet is ranova(LM_AC_BL) # Replicate is not significant LM_AC_BL <-lmer(dAW.gd ~ MoteGen * DayF + (1|Fra), data= BW.AC)
anova(LM_AC_BL) # Day alone is not significant, but Day:Genet is ranova(LM_AC_BL)# Pairwise comparisons
BW_AC_BL.emmc<-emmeans(LM_AC_BL, ~MoteGen | DayF)
#contrast(BW_AC_BL.emmc, "tukey")
BW_AC_BL.emmc<-emmeans(LM_AC_BL, ~DayF|MoteGen)
#contrast(BW_AC_BL.emmc, "tukey")
#Tukey groups
BW_AC_BL.emmc<-emmeans(LM_AC_BL, ~MoteGen * DayF)
BW_AC_BL_groups<-cld(BW_AC_BL.emmc)
BW_AC_BL_groups<-BW_AC_BL_groups[order(BW_AC_BL_groups$DayF,BW_AC_BL_groups$MoteGen),]
BW_AC_BL_groups #write.csv(BW_AC_BL_groups, "Outputs/BW_AC_BL_groups", row.names = F)- Ambient heat stress (day 100)
summary(BW.100)## Gen Fra Treatment Replicate Time
## G-O :21 Min. :102.0 A :32 R1:35 7.Feb-Feb:63
## Green :22 1st Qu.:151.0 N :15 R2:28 1.Sep-Oct: 0
## Orange: 7 Median :187.0 N+P:16 2.Oct-Nov: 0
## R-O : 7 Mean :197.5 3.Nov-Dec: 0
## R-Y : 4 3rd Qu.:249.0 4.Dec-Jan: 0
## Yellow: 2 Max. :289.0 5.Jan-Feb: 0
## (Other) : 0
## Period Initial.Date Final.Date
## 1.Sep-Oct : 0 Min. :2018-02-14 Min. :2018-02-23
## 4.Dec-Jan : 0 1st Qu.:2018-02-14 1st Qu.:2018-02-23
## Bleaching : 0 Median :2018-02-14 Median :2018-02-23
## Jan-Feb : 0 Mean :2018-02-14 Mean :2018-02-23
## Nov-Dec : 0 3rd Qu.:2018-02-14 3rd Qu.:2018-02-23
## Oct-Nov : 0 Max. :2018-02-14 Max. :2018-02-23
## Ramping :63
## dBW.g dAW.gd Original.Collection.Date MoteGen
## Min. :-1.6664 Min. :-1.4820 02/23/08: 9 G_48:22
## 1st Qu.:-0.5611 1st Qu.:-0.3182 03/18/10: 7 G_62:21
## Median : 0.5736 Median : 0.8028 03/26/10: 0 G_31: 7
## Mean : 0.4222 Mean : 0.6424 05/10/10: 0 G_08: 2
## 3rd Qu.: 1.4625 3rd Qu.: 1.6950 05/21/13:21 G_07: 7
## Max. : 2.9316 Max. : 3.1725 05/28/10:22 G_50: 4
## 07/30/10: 4
## Collection.Habitat Lat Lon
## Midchannel :34 Min. :24.52 Min. :-81.64
## Offshore :29 1st Qu.:24.52 1st Qu.:-81.64
## Reef Margin: 0 Median :24.53 Median :-81.58
## Mean :24.54 Mean :-81.57
## 3rd Qu.:24.55 3rd Qu.:-81.52
## Max. :24.58 Max. :-81.43
##
## Genotype WB Days DayF TimePoint
## Green :22 R: 7 Min. :100 100 :63 Min. :6
## Green and Orange:21 S:56 1st Qu.:100 -28 : 0 1st Qu.:6
## Orange : 7 Median :100 -7 : 0 Median :6
## Red and Orange : 7 Mean :100 28 : 0 Mean :6
## Red and Yellow : 4 3rd Qu.:100 62 : 0 3rd Qu.:6
## Yellow : 2 Max. :100 75 : 0 Max. :6
## (Other) : 0 (Other): 0
## Nutrients Colony SampleID
## Length:63 G-O :21 Length:63
## Class :character Green :22 Class :character
## Mode :character Orange: 7 Mode :character
## R-O : 7
## R-Y : 4
## Yellow: 2
## BW.100.C<-subset(BW.100, Treatment=="A")
summary(BW.100.C)## Gen Fra Treatment Replicate Time
## G-O :7 Min. :102.0 A :32 R1:17 7.Feb-Feb:32
## Green :8 1st Qu.:149.8 N : 0 R2:15 1.Sep-Oct: 0
## Orange:4 Median :177.5 N+P: 0 2.Oct-Nov: 0
## R-O :7 Mean :189.8 3.Nov-Dec: 0
## R-Y :4 3rd Qu.:235.2 4.Dec-Jan: 0
## Yellow:2 Max. :287.0 5.Jan-Feb: 0
## (Other) : 0
## Period Initial.Date Final.Date
## 1.Sep-Oct : 0 Min. :2018-02-14 Min. :2018-02-23
## 4.Dec-Jan : 0 1st Qu.:2018-02-14 1st Qu.:2018-02-23
## Bleaching : 0 Median :2018-02-14 Median :2018-02-23
## Jan-Feb : 0 Mean :2018-02-14 Mean :2018-02-23
## Nov-Dec : 0 3rd Qu.:2018-02-14 3rd Qu.:2018-02-23
## Oct-Nov : 0 Max. :2018-02-14 Max. :2018-02-23
## Ramping :32
## dBW.g dAW.gd Original.Collection.Date MoteGen
## Min. :0.5736 Min. :0.7557 02/23/08:9 G_48:8
## 1st Qu.:1.1549 1st Qu.:1.3866 03/18/10:4 G_62:7
## Median :1.4625 Median :1.6950 03/26/10:0 G_31:4
## Mean :1.4597 Mean :1.6721 05/10/10:0 G_08:2
## 3rd Qu.:1.7315 3rd Qu.:1.9698 05/21/13:7 G_07:7
## Max. :2.9316 Max. :3.1725 05/28/10:8 G_50:4
## 07/30/10:4
## Collection.Habitat Lat Lon
## Midchannel :20 Min. :24.52 Min. :-81.64
## Offshore :12 1st Qu.:24.53 1st Qu.:-81.58
## Reef Margin: 0 Median :24.55 Median :-81.53
## Mean :24.54 Mean :-81.54
## 3rd Qu.:24.56 3rd Qu.:-81.50
## Max. :24.58 Max. :-81.43
##
## Genotype WB Days DayF TimePoint
## Green :8 R: 4 Min. :100 100 :32 Min. :6
## Green and Orange:7 S:28 1st Qu.:100 -28 : 0 1st Qu.:6
## Red and Orange :7 Median :100 -7 : 0 Median :6
## Orange :4 Mean :100 28 : 0 Mean :6
## Red and Yellow :4 3rd Qu.:100 62 : 0 3rd Qu.:6
## Yellow :2 Max. :100 75 : 0 Max. :6
## (Other) :0 (Other): 0
## Nutrients Colony SampleID
## Length:32 G-O :7 Length:32
## Class :character Green :8 Class :character
## Mode :character Orange:4 Mode :character
## R-O :7
## R-Y :4
## Yellow:2
## # Model
LM_A_H.100 <- lmer(dAW.gd ~ MoteGen + (1|Replicate), data= BW.100.C)
anova(LM_A_H.100) ranova(LM_A_H.100) # Replicate is not significant# Pairwise comparisons
BW_A_H100.emmc<-emmeans(LM_A_H.100, ~MoteGen)
#contrast(BW_A_H100.emmc, "tukey")
#Tukey groups
BW_A_H100_groups<-cld(BW_A_H100.emmc)
BW_A_H100_groups #write.csv(BW_A_H100_groups, "Outputs/BW_A_H100_groups.csv", row.names = F)
Stats_100_a<-BW.100.C %>%
group_by(MoteGen) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
Stats_100_aBW_91<- ggplot(subset(BW.91, Treatment=="A"), aes (MoteGen, dAW.gd, colour=MoteGen)) +
ggtitle("(c) Heat (day 91)")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2) +
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) + ggthe_bw +
scale_y_continuous(name=("Growth (mg/ g d) "),
limits = c(0,6),
breaks = seq(0, 6, by=0.5))+
theme(legend.position="none",
legend.title = element_blank(),
strip.background =element_rect(fill=NA)) #+
facet_grid(Treatment~.)## <ggproto object: Class FacetGrid, Facet, gg>
## compute_layout: function
## draw_back: function
## draw_front: function
## draw_labels: function
## draw_panels: function
## finish_data: function
## init_scales: function
## map_data: function
## params: list
## setup_data: function
## setup_params: function
## shrink: TRUE
## train_scales: function
## vars: function
## super: <ggproto object: Class FacetGrid, Facet, gg>BW_100<- ggplot(subset(BW.100, Treatment=="A"), aes (MoteGen, dAW.gd, colour=MoteGen)) +
ggtitle("(d) Heat (day 100)")+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2) +
stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) + ggthe_bw +
scale_y_continuous(name=("Growth (mg/ g d) "),
limits = c(0,6),
breaks = seq(0, 6, by=0.5))+
theme(legend.position="none",
legend.title = element_blank(),
strip.background =element_rect(fill=NA)) #+facet_grid(Treatment~.)
HeatEffect<-grid.arrange(BW_91, BW_100, nrow=1)
Figure S2
All_ambient<-grid.arrange(BW_Init,BW_Final, BW_91, BW_100, nrow=1)
#ggsave(file="Outputs/S_Genotypes_Ambient.svg", plot=All_ambient, width=8, height=5.5)- Selected timepoints = BL + Ambient 75 (Control) and 100 (Heat)
# Subset_data
BW.A<-BW.Tall
BW.A$Treatment[BW.A$DayF=="-28"]<-"A"
BW.A<-subset(BW.A, Treatment=="A")
BW.A<-subset(BW.A, DayF!="28")
BW.A<-subset(BW.A, DayF!="62")
# BW.A<-subset(BW.A, DayF!="91")
BW.A<-droplevels(BW.A)
# Summary
summary(BW.A)## MoteGen DayF Gen Fra Treatment Replicate
## G_48:54 -28:118 G-O :52 Min. :101.0 A:221 R1:116
## G_62:52 75 : 39 Green :54 1st Qu.:152.0 R2:105
## G_31:29 91 : 32 Orange:29 Median :183.0
## G_08:12 100: 32 R-O :49 Mean :195.7
## G_07:49 R-Y :25 3rd Qu.:249.0
## G_50:25 Yellow:12 Max. :300.0
## Time Period Initial.Date Final.Date
## 1.Sep-Oct:118 1.Sep-Oct:118 Min. :2017-08-30 Min. :2017-10-18
## 5.Jan-Feb: 39 Jan-Feb : 39 1st Qu.:2017-08-30 1st Qu.:2017-10-18
## 6.Feb-Feb: 32 Ramping : 64 Median :2017-08-30 Median :2017-10-18
## 7.Feb-Feb: 32 Mean :2017-11-08 Mean :2017-12-10
## 3rd Qu.:2018-01-29 3rd Qu.:2018-02-14
## Max. :2018-02-14 Max. :2018-02-23
## dBW.g dAW.gd Original.Collection.Date
## Min. :0.5736 Min. :0.7557 02/23/08:61
## 1st Qu.:2.2765 1st Qu.:2.2537 03/18/10:29
## Median :3.1131 Median :3.1206 05/21/13:52
## Mean :3.1393 Mean :3.1412 05/28/10:54
## 3rd Qu.:3.9952 3rd Qu.:3.9831 07/30/10:25
## Max. :6.5930 Max. :6.5938
## Collection.Habitat Lat Lon
## Midchannel:138 Min. :24.52 Min. :-81.64
## Offshore : 83 1st Qu.:24.53 1st Qu.:-81.58
## Median :24.55 Median :-81.53
## Mean :24.54 Mean :-81.55
## 3rd Qu.:24.56 3rd Qu.:-81.50
## Max. :24.58 Max. :-81.43
## Genotype WB Days TimePoint
## Green :54 R: 29 Min. :-28.00 Min. :1.000
## Green and Orange:52 S:192 1st Qu.:-28.00 1st Qu.:1.000
## Orange :29 Median :-28.00 Median :1.000
## Red and Orange :49 Mean : 25.94 Mean :2.833
## Red and Yellow :25 3rd Qu.: 91.00 3rd Qu.:5.000
## Yellow :12 Max. :100.00 Max. :6.000
## Nutrients Colony SampleID TreatmentG
## Length:221 G-O :52 Length:221 A:221
## Class :character Green :54 Class :character
## Mode :character Orange:29 Mode :character
## R-O :49
## R-Y :25
## Yellow:12
## Nutrients2 mean sd Difference
## Length:221 Min. :1.208 Min. :0.2020 Min. :-2.1530000
## Class :character 1st Qu.:2.490 1st Qu.:0.4430 1st Qu.:-0.3124000
## Mode :character Median :3.167 Median :0.6330 Median : 0.0368000
## Mean :3.141 Mean :0.6763 Mean : 0.0000072
## 3rd Qu.:3.924 3rd Qu.:0.9130 3rd Qu.: 0.3779000
## Max. :4.294 Max. :1.0870 Max. : 2.6698000
## percentage
## Min. :-55.56189
## 1st Qu.:-10.76386
## Median : 1.13198
## Mean : 0.00009
## 3rd Qu.: 13.70665
## Max. : 68.03772 Stats_Ambient<-BW.A %>%
group_by(MoteGen, DayF) %>%
get_summary_stats(dAW.gd, type = "mean_sd")
Stats_Ambient# Model
LM_Ambinet <- lmer(dAW.gd ~ MoteGen * DayF + (1|Replicate), data= BW.A)
isSingular(LM_Ambinet)## [1] TRUE anova(LM_Ambinet) ranova(LM_Ambinet) # Replicate is not significant LM_Ambinet <-lmer(dAW.gd ~ MoteGen * DayF + (1|Fra), data= BW.A)
anova(LM_Ambinet) # Day alone is not significant, but Day:Genet is ranova(LM_Ambinet)# Pairwise comparisons
BW_Ambient.emmc<-emmeans(LM_Ambinet, ~MoteGen | DayF)
#contrast(BW_Ambient.emmc, "tukey")
BW_Ambient.emmc<-emmeans(LM_Ambinet, ~DayF|MoteGen)
#contrast(BW_Ambient.emmc, "tukey")
#Tukey groups
BW_AC_BL.emmc<-emmeans(LM_Ambinet, ~MoteGen * DayF)
BW_AC_BL_groups<-cld(BW_AC_BL.emmc)
BW_AC_BL_groups<-BW_AC_BL_groups[
order(BW_AC_BL_groups$DayF,BW_AC_BL_groups$MoteGen),]
BW_AC_BL_groups #write.csv(BW_AC_BL_groups, "Outputs/BW_AC_BL_groups", row.names = F)Treatment and genets
Control all treatments
# 1. Model
LM_4 <- lmer(dAW.gd ~ Treatment*DayF*MoteGen +
(1|Fra), REML=TRUE, data=BW.nutrients)
# Select model
anova(LM_4) ranova(LM_4)# 2. Predict values:
pred_nutrients <- predict(LM_4,re.form = NA)
#3. Bootstrap CI:
nutrients.boot1 <- bootMer(LM_4, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(nutrients.boot1$t, 2, sd)
CI.lo_1 <- pred_nutrients - std.err*1.96
CI.hi_1 <- pred_nutrients + std.err*1.96
#Plot
Model_nutrinets_plot<- ggplot(
BW.nutrients, aes(x = Days, y = dAW.gd, colour = MoteGen)) +
geom_line(aes(y = pred_nutrients),size=2) +
#geom_point(aes(fill=factor(Treatment)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Growth")),
limits = c(0, 6),
breaks = seq(0, 6, by=1), expand = c(0,0))+
scale_x_continuous("Days in the experiment", limits = c(-30, 78),
breaks = seq(-30, 76, by=15), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=0.5) + ggthe_bw
Model_nutrinets_plot + facet_grid(~Treatment)
# Pairwise comparisons
BW.emmc<-emmeans(LM_4, ~Treatment|DayF|MoteGen)
contrast(BW.emmc, "tukey")## DayF = -28, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N -0.15828 0.316 258 -0.500 0.8712
## A - N+P 0.21157 0.316 258 0.669 0.7819
## N - N+P 0.36984 0.325 258 1.139 0.4907
##
## DayF = -7, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N -0.87039 0.316 258 -2.751 0.0174
## A - N+P -0.37396 0.316 258 -1.182 0.4650
## N - N+P 0.49642 0.325 258 1.529 0.2789
##
## DayF = 28, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N 0.63171 0.316 258 1.997 0.1151
## A - N+P 0.25054 0.316 258 0.792 0.7083
## N - N+P -0.38117 0.325 258 -1.174 0.4696
##
## DayF = 62, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N 1.50914 0.316 258 4.770 <.0001
## A - N+P 1.58785 0.316 258 5.019 <.0001
## N - N+P 0.07871 0.325 258 0.242 0.9681
##
## DayF = 75, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N 2.47114 0.316 258 7.811 <.0001
## A - N+P 2.20803 0.316 258 6.979 <.0001
## N - N+P -0.26311 0.325 258 -0.811 0.6968
##
## DayF = -28, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 0.58742 0.316 258 1.857 0.1535
## A - N+P 0.44863 0.316 258 1.418 0.3331
## N - N+P -0.13879 0.308 258 -0.451 0.8941
##
## DayF = -7, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N -0.15313 0.316 258 -0.484 0.8789
## A - N+P -0.37154 0.328 283 -1.134 0.4941
## N - N+P -0.21840 0.320 285 -0.683 0.7733
##
## DayF = 28, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 1.08356 0.323 272 3.358 0.0026
## A - N+P 0.57995 0.323 272 1.797 0.1725
## N - N+P -0.50360 0.308 258 -1.635 0.2327
##
## DayF = 62, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 1.55948 0.316 258 4.929 <.0001
## A - N+P 1.43657 0.316 258 4.541 <.0001
## N - N+P -0.12291 0.308 258 -0.399 0.9160
##
## DayF = 75, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 2.37987 0.316 258 7.522 <.0001
## A - N+P 2.07784 0.316 258 6.568 <.0001
## N - N+P -0.30203 0.308 258 -0.981 0.5897
##
## DayF = -28, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N -0.03156 0.417 258 -0.076 0.9968
## A - N+P -0.76568 0.435 258 -1.758 0.1859
## N - N+P -0.73412 0.417 258 -1.761 0.1850
##
## DayF = -7, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N -0.44643 0.417 258 -1.071 0.5331
## A - N+P -0.60062 0.435 258 -1.379 0.3533
## N - N+P -0.15419 0.417 258 -0.370 0.9274
##
## DayF = 28, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N 1.70606 0.417 258 4.092 0.0002
## A - N+P 0.29158 0.435 258 0.670 0.7814
## N - N+P -1.41448 0.417 258 -3.392 0.0023
##
## DayF = 62, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N 2.97343 0.417 258 7.131 <.0001
## A - N+P 3.25682 0.435 258 7.478 <.0001
## N - N+P 0.28339 0.417 258 0.680 0.7755
##
## DayF = 75, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N 3.32735 0.417 258 7.980 <.0001
## A - N+P 3.04824 0.435 258 6.999 <.0001
## N - N+P -0.27911 0.417 258 -0.669 0.7815
##
## DayF = -28, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N -1.21151 0.666 300 -1.819 0.1652
## A - N+P -0.75803 0.666 300 -1.138 0.4914
## N - N+P 0.45348 0.643 353 0.705 0.7604
##
## DayF = -7, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N -1.21930 0.629 258 -1.940 0.1297
## A - N+P -0.91300 0.629 258 -1.452 0.3157
## N - N+P 0.30630 0.562 258 0.545 0.8492
##
## DayF = 28, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N -0.62010 0.629 258 -0.987 0.5861
## A - N+P -0.45567 0.629 258 -0.725 0.7490
## N - N+P 0.16443 0.562 258 0.292 0.9540
##
## DayF = 62, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N 1.32383 0.629 258 2.106 0.0906
## A - N+P 0.81403 0.629 258 1.295 0.3993
## N - N+P -0.50980 0.562 258 -0.907 0.6365
##
## DayF = 75, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N 1.62885 0.629 258 2.591 0.0272
## A - N+P 1.49258 0.629 258 2.375 0.0479
## N - N+P -0.13627 0.562 258 -0.242 0.9681
##
## DayF = -28, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N -0.21880 0.325 258 -0.674 0.7788
## A - N+P -0.37763 0.335 258 -1.129 0.4973
## N - N+P -0.15883 0.335 258 -0.475 0.8833
##
## DayF = -7, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N -0.65207 0.325 258 -2.009 0.1121
## A - N+P -0.68462 0.335 258 -2.046 0.1034
## N - N+P -0.03255 0.335 258 -0.097 0.9948
##
## DayF = 28, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N 1.30548 0.325 258 4.022 0.0002
## A - N+P 0.88833 0.335 258 2.655 0.0228
## N - N+P -0.41715 0.335 258 -1.247 0.4268
##
## DayF = 62, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N 1.56877 0.331 272 4.743 <.0001
## A - N+P 1.75589 0.341 271 5.156 <.0001
## N - N+P 0.18712 0.335 258 0.559 0.8418
##
## DayF = 75, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N 2.04447 0.331 272 6.181 <.0001
## A - N+P 1.33779 0.335 258 3.998 0.0002
## N - N+P -0.70668 0.341 271 -2.075 0.0970
##
## DayF = -28, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N -1.02727 0.487 258 -2.110 0.0898
## A - N+P 0.45245 0.462 258 0.980 0.5905
## N - N+P 1.47973 0.462 258 3.203 0.0043
##
## DayF = -7, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N -0.50118 0.487 258 -1.029 0.5590
## A - N+P 0.00498 0.462 258 0.011 0.9999
## N - N+P 0.50616 0.462 258 1.096 0.5175
##
## DayF = 28, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N 1.08945 0.487 258 2.238 0.0669
## A - N+P 1.50390 0.462 258 3.256 0.0037
## N - N+P 0.41445 0.462 258 0.897 0.6425
##
## DayF = 62, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N 2.55840 0.487 258 5.254 <.0001
## A - N+P 2.24356 0.478 282 4.693 <.0001
## N - N+P -0.31484 0.478 282 -0.659 0.7876
##
## DayF = 75, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N 3.07633 0.675 455 4.558 <.0001
## A - N+P 2.81281 0.478 282 5.884 <.0001
## N - N+P -0.26351 0.669 468 -0.394 0.9180
##
## P value adjustment: tukey method for comparing a family of 3 estimates BW_groups<-cld(BW.emmc)
#BW_groups<-as.data.frame(BW_groups[complete.cases(BW_groups),])
BW_groups<-BW_groups[order(BW_groups$Day,BW_groups$Treatment),]
BW_groups #write.csv(BW_groups, "Outputs/BW_groups_corrected.csv", row.names = F)Pooled N and N=P Control
BW.Tall$Nutrients <- "Nutrients"
BW.Tall$Nutrients[BW.Tall$Treatment == "A"] <- "Ambient"
BW.Nutrients<-subset(BW.Tall, Days<91)
# 1. Model
LM_5 <- lmer(dAW.gd ~ Nutrients*DayF*MoteGen +
(1|Fra), REML=TRUE, data=BW.Nutrients)
# Select model
anova(LM_5) ranova(LM_5)# 2. Predict values:
pred_nutrients2 <- predict(LM_5, re.form = NA)
#3. Bootstrap CI:
nutrients.boot2 <- bootMer(LM_5, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
std.err <- apply(nutrients.boot2$t, 2, sd)
CI.lo_2 <- pred_nutrients2 - std.err*1.96
CI.hi_2 <- pred_nutrients2 + std.err*1.96
#Plot
Model_nutrinets_plot2<- ggplot(BW.Nutrients,
aes(x = Days, y = dAW.gd,
colour = MoteGen, shape=Nutrients)) +
geom_line(aes(y = pred_nutrients2),size=2) +
#geom_point(aes(fill=factor(Nutrients)),
# shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
geom_ribbon(aes(ymin = CI.lo_2, ymax = CI.hi_2),
size=2, alpha = 0.1, linetype = 0) +
#scale_color_manual(values=my_colours) +
#scale_fill_manual(values=my_colours) +
scale_y_continuous(name=expression(~italic("Growth")),
limits = c(-2, 6),
breaks = seq(-2, 6, by=1), expand = c(0,0))+
scale_x_continuous("Days in the experiment", limits = c(-30, 110),
breaks = seq(-30, 110, by=30), expand = c(0,0))+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
position = position_dodge(1) )+
stat_summary(fun.y=mean, geom="line", position = position_dodge(1),
linetype=1, alpha=0.5) + ggthe_bw
Model_nutrinets_plot2 + facet_grid(~Nutrients)
# Pairwise comparisons
contrast(BW.emmc, "tukey")## DayF = -28, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N -0.15828 0.316 258 -0.500 0.8712
## A - N+P 0.21157 0.316 258 0.669 0.7819
## N - N+P 0.36984 0.325 258 1.139 0.4907
##
## DayF = -7, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N -0.87039 0.316 258 -2.751 0.0174
## A - N+P -0.37396 0.316 258 -1.182 0.4650
## N - N+P 0.49642 0.325 258 1.529 0.2789
##
## DayF = 28, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N 0.63171 0.316 258 1.997 0.1151
## A - N+P 0.25054 0.316 258 0.792 0.7083
## N - N+P -0.38117 0.325 258 -1.174 0.4696
##
## DayF = 62, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N 1.50914 0.316 258 4.770 <.0001
## A - N+P 1.58785 0.316 258 5.019 <.0001
## N - N+P 0.07871 0.325 258 0.242 0.9681
##
## DayF = 75, MoteGen = G_48:
## contrast estimate SE df t.ratio p.value
## A - N 2.47114 0.316 258 7.811 <.0001
## A - N+P 2.20803 0.316 258 6.979 <.0001
## N - N+P -0.26311 0.325 258 -0.811 0.6968
##
## DayF = -28, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 0.58742 0.316 258 1.857 0.1535
## A - N+P 0.44863 0.316 258 1.418 0.3331
## N - N+P -0.13879 0.308 258 -0.451 0.8941
##
## DayF = -7, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N -0.15313 0.316 258 -0.484 0.8789
## A - N+P -0.37154 0.328 283 -1.134 0.4941
## N - N+P -0.21840 0.320 285 -0.683 0.7733
##
## DayF = 28, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 1.08356 0.323 272 3.358 0.0026
## A - N+P 0.57995 0.323 272 1.797 0.1725
## N - N+P -0.50360 0.308 258 -1.635 0.2327
##
## DayF = 62, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 1.55948 0.316 258 4.929 <.0001
## A - N+P 1.43657 0.316 258 4.541 <.0001
## N - N+P -0.12291 0.308 258 -0.399 0.9160
##
## DayF = 75, MoteGen = G_62:
## contrast estimate SE df t.ratio p.value
## A - N 2.37987 0.316 258 7.522 <.0001
## A - N+P 2.07784 0.316 258 6.568 <.0001
## N - N+P -0.30203 0.308 258 -0.981 0.5897
##
## DayF = -28, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N -0.03156 0.417 258 -0.076 0.9968
## A - N+P -0.76568 0.435 258 -1.758 0.1859
## N - N+P -0.73412 0.417 258 -1.761 0.1850
##
## DayF = -7, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N -0.44643 0.417 258 -1.071 0.5331
## A - N+P -0.60062 0.435 258 -1.379 0.3533
## N - N+P -0.15419 0.417 258 -0.370 0.9274
##
## DayF = 28, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N 1.70606 0.417 258 4.092 0.0002
## A - N+P 0.29158 0.435 258 0.670 0.7814
## N - N+P -1.41448 0.417 258 -3.392 0.0023
##
## DayF = 62, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N 2.97343 0.417 258 7.131 <.0001
## A - N+P 3.25682 0.435 258 7.478 <.0001
## N - N+P 0.28339 0.417 258 0.680 0.7755
##
## DayF = 75, MoteGen = G_31:
## contrast estimate SE df t.ratio p.value
## A - N 3.32735 0.417 258 7.980 <.0001
## A - N+P 3.04824 0.435 258 6.999 <.0001
## N - N+P -0.27911 0.417 258 -0.669 0.7815
##
## DayF = -28, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N -1.21151 0.666 300 -1.819 0.1652
## A - N+P -0.75803 0.666 300 -1.138 0.4914
## N - N+P 0.45348 0.643 353 0.705 0.7604
##
## DayF = -7, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N -1.21930 0.629 258 -1.940 0.1297
## A - N+P -0.91300 0.629 258 -1.452 0.3157
## N - N+P 0.30630 0.562 258 0.545 0.8492
##
## DayF = 28, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N -0.62010 0.629 258 -0.987 0.5861
## A - N+P -0.45567 0.629 258 -0.725 0.7490
## N - N+P 0.16443 0.562 258 0.292 0.9540
##
## DayF = 62, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N 1.32383 0.629 258 2.106 0.0906
## A - N+P 0.81403 0.629 258 1.295 0.3993
## N - N+P -0.50980 0.562 258 -0.907 0.6365
##
## DayF = 75, MoteGen = G_08:
## contrast estimate SE df t.ratio p.value
## A - N 1.62885 0.629 258 2.591 0.0272
## A - N+P 1.49258 0.629 258 2.375 0.0479
## N - N+P -0.13627 0.562 258 -0.242 0.9681
##
## DayF = -28, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N -0.21880 0.325 258 -0.674 0.7788
## A - N+P -0.37763 0.335 258 -1.129 0.4973
## N - N+P -0.15883 0.335 258 -0.475 0.8833
##
## DayF = -7, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N -0.65207 0.325 258 -2.009 0.1121
## A - N+P -0.68462 0.335 258 -2.046 0.1034
## N - N+P -0.03255 0.335 258 -0.097 0.9948
##
## DayF = 28, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N 1.30548 0.325 258 4.022 0.0002
## A - N+P 0.88833 0.335 258 2.655 0.0228
## N - N+P -0.41715 0.335 258 -1.247 0.4268
##
## DayF = 62, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N 1.56877 0.331 272 4.743 <.0001
## A - N+P 1.75589 0.341 271 5.156 <.0001
## N - N+P 0.18712 0.335 258 0.559 0.8418
##
## DayF = 75, MoteGen = G_07:
## contrast estimate SE df t.ratio p.value
## A - N 2.04447 0.331 272 6.181 <.0001
## A - N+P 1.33779 0.335 258 3.998 0.0002
## N - N+P -0.70668 0.341 271 -2.075 0.0970
##
## DayF = -28, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N -1.02727 0.487 258 -2.110 0.0898
## A - N+P 0.45245 0.462 258 0.980 0.5905
## N - N+P 1.47973 0.462 258 3.203 0.0043
##
## DayF = -7, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N -0.50118 0.487 258 -1.029 0.5590
## A - N+P 0.00498 0.462 258 0.011 0.9999
## N - N+P 0.50616 0.462 258 1.096 0.5175
##
## DayF = 28, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N 1.08945 0.487 258 2.238 0.0669
## A - N+P 1.50390 0.462 258 3.256 0.0037
## N - N+P 0.41445 0.462 258 0.897 0.6425
##
## DayF = 62, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N 2.55840 0.487 258 5.254 <.0001
## A - N+P 2.24356 0.478 282 4.693 <.0001
## N - N+P -0.31484 0.478 282 -0.659 0.7876
##
## DayF = 75, MoteGen = G_50:
## contrast estimate SE df t.ratio p.value
## A - N 3.07633 0.675 455 4.558 <.0001
## A - N+P 2.81281 0.478 282 5.884 <.0001
## N - N+P -0.26351 0.669 468 -0.394 0.9180
##
## P value adjustment: tukey method for comparing a family of 3 estimates BW.emmc<-emmeans(LM_5, ~MoteGen|DayF|Nutrients)
contrast(BW.emmc, "tukey")## DayF = -28, Nutrients = Ambient:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 -0.3801 0.331 258 -1.148 0.8607
## G_48 - G_31 -0.8996 0.395 258 -2.280 0.2062
## G_48 - G_08 -0.2483 0.558 258 -0.445 0.9978
## G_48 - G_07 -0.5887 0.331 258 -1.778 0.4816
## G_48 - G_50 -0.9486 0.426 258 -2.225 0.2298
## G_62 - G_31 -0.5196 0.402 258 -1.293 0.7888
## G_62 - G_08 0.1318 0.563 258 0.234 0.9999
## G_62 - G_07 -0.2087 0.340 258 -0.614 0.9899
## G_62 - G_50 -0.5685 0.433 258 -1.313 0.7777
## G_31 - G_08 0.6514 0.603 258 1.081 0.8888
## G_31 - G_07 0.3109 0.402 258 0.774 0.9717
## G_31 - G_50 -0.0489 0.483 258 -0.101 1.0000
## G_08 - G_07 -0.3404 0.563 258 -0.604 0.9907
## G_08 - G_50 -0.7003 0.624 258 -1.122 0.8718
## G_07 - G_50 -0.3599 0.433 258 -0.831 0.9615
##
## DayF = 28, Nutrients = Ambient:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.0851 0.338 271 0.252 0.9999
## G_48 - G_31 -1.0567 0.395 258 -2.678 0.0833
## G_48 - G_08 1.7609 0.558 258 3.155 0.0219
## G_48 - G_07 0.4640 0.331 258 1.402 0.7261
## G_48 - G_50 -0.1719 0.426 258 -0.403 0.9986
## G_62 - G_31 -1.1418 0.408 267 -2.800 0.0606
## G_62 - G_08 1.6759 0.567 263 2.953 0.0397
## G_62 - G_07 0.3790 0.347 271 1.093 0.8838
## G_62 - G_50 -0.2569 0.438 266 -0.586 0.9919
## G_31 - G_08 2.8176 0.603 258 4.674 0.0001
## G_31 - G_07 1.5207 0.402 258 3.784 0.0026
## G_31 - G_50 0.8848 0.483 258 1.831 0.4479
## G_08 - G_07 -1.2969 0.563 258 -2.302 0.1968
## G_08 - G_50 -1.9328 0.624 258 -3.097 0.0261
## G_07 - G_50 -0.6359 0.433 258 -1.469 0.6846
##
## DayF = 62, Nutrients = Ambient:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.3379 0.331 258 1.021 0.9107
## G_48 - G_31 -0.4649 0.395 258 -1.178 0.8470
## G_48 - G_08 2.0845 0.558 258 3.735 0.0031
## G_48 - G_07 1.5730 0.338 271 4.651 0.0001
## G_48 - G_50 1.2373 0.426 258 2.903 0.0458
## G_62 - G_31 -0.8029 0.402 258 -1.998 0.3464
## G_62 - G_08 1.7466 0.563 258 3.101 0.0258
## G_62 - G_07 1.2351 0.347 271 3.563 0.0057
## G_62 - G_50 0.8994 0.433 258 2.077 0.3026
## G_31 - G_08 2.5495 0.603 258 4.229 0.0005
## G_31 - G_07 2.0380 0.408 267 4.997 <.0001
## G_31 - G_50 1.7022 0.483 258 3.522 0.0067
## G_08 - G_07 -0.5115 0.567 263 -0.901 0.9459
## G_08 - G_50 -0.8472 0.624 258 -1.358 0.7521
## G_07 - G_50 -0.3357 0.438 266 -0.766 0.9730
##
## DayF = 75, Nutrients = Ambient:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.2734 0.331 258 0.826 0.9625
## G_48 - G_31 0.0357 0.395 258 0.091 1.0000
## G_48 - G_08 2.0726 0.558 258 3.714 0.0034
## G_48 - G_07 1.6449 0.331 258 4.968 <.0001
## G_48 - G_50 1.4158 0.426 258 3.321 0.0130
## G_62 - G_31 -0.2377 0.402 258 -0.591 0.9915
## G_62 - G_08 1.7991 0.563 258 3.194 0.0194
## G_62 - G_07 1.3714 0.340 258 4.038 0.0010
## G_62 - G_50 1.1423 0.433 258 2.638 0.0919
## G_31 - G_08 2.0368 0.603 258 3.379 0.0108
## G_31 - G_07 1.6091 0.402 258 4.004 0.0011
## G_31 - G_50 1.3800 0.483 258 2.855 0.0522
## G_08 - G_07 -0.4277 0.563 258 -0.759 0.9739
## G_08 - G_50 -0.6568 0.624 258 -1.053 0.8994
## G_07 - G_50 -0.2291 0.433 258 -0.529 0.9950
##
## DayF = -28, Nutrients = Nutrients:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.1113 0.234 258 0.476 0.9970
## G_48 - G_31 -1.2915 0.276 258 -4.684 0.0001
## G_48 - G_08 -1.2345 0.380 323 -3.251 0.0159
## G_48 - G_07 -0.9089 0.244 258 -3.730 0.0032
## G_48 - G_50 -1.1804 0.294 258 -4.013 0.0011
## G_62 - G_31 -1.4029 0.270 258 -5.187 <.0001
## G_62 - G_08 -1.3458 0.376 325 -3.581 0.0053
## G_62 - G_07 -1.0202 0.238 258 -4.292 0.0004
## G_62 - G_50 -1.2917 0.289 258 -4.466 0.0002
## G_31 - G_08 0.0571 0.403 316 0.142 1.0000
## G_31 - G_07 0.3826 0.279 258 1.372 0.7435
## G_31 - G_50 0.1111 0.324 258 0.343 0.9994
## G_08 - G_07 0.3256 0.382 323 0.853 0.9572
## G_08 - G_50 0.0541 0.416 312 0.130 1.0000
## G_07 - G_50 -0.2715 0.297 258 -0.914 0.9427
##
## DayF = 28, Nutrients = Nutrients:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.4596 0.234 258 1.963 0.3665
## G_48 - G_31 -0.4347 0.276 258 -1.577 0.6148
## G_48 - G_08 0.7819 0.340 258 2.302 0.1970
## G_48 - G_07 1.1321 0.244 258 4.646 0.0001
## G_48 - G_50 0.7067 0.294 258 2.402 0.1592
## G_62 - G_31 -0.8943 0.270 258 -3.307 0.0136
## G_62 - G_08 0.3223 0.335 258 0.961 0.9297
## G_62 - G_07 0.6725 0.238 258 2.829 0.0561
## G_62 - G_50 0.2471 0.289 258 0.854 0.9567
## G_31 - G_08 1.2166 0.366 258 3.327 0.0128
## G_31 - G_07 1.5668 0.279 258 5.620 <.0001
## G_31 - G_50 1.1414 0.324 258 3.525 0.0066
## G_08 - G_07 0.3502 0.342 258 1.023 0.9098
## G_08 - G_50 -0.0752 0.380 258 -0.198 1.0000
## G_07 - G_50 -0.4254 0.297 258 -1.432 0.7074
##
## DayF = 62, Nutrients = Nutrients:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.2874 0.234 258 1.228 0.8229
## G_48 - G_31 1.0888 0.276 258 3.949 0.0014
## G_48 - G_08 1.6050 0.340 258 4.725 0.0001
## G_48 - G_07 1.6719 0.244 258 6.861 <.0001
## G_48 - G_50 2.1501 0.303 274 7.103 <.0001
## G_62 - G_31 0.8013 0.270 258 2.963 0.0387
## G_62 - G_08 1.3175 0.335 258 3.928 0.0015
## G_62 - G_07 1.3845 0.238 258 5.825 <.0001
## G_62 - G_50 1.8627 0.298 275 6.253 <.0001
## G_31 - G_08 0.5162 0.366 258 1.412 0.7200
## G_31 - G_07 0.5831 0.279 258 2.092 0.2949
## G_31 - G_50 1.0613 0.332 272 3.200 0.0190
## G_08 - G_07 0.0670 0.342 258 0.196 1.0000
## G_08 - G_50 0.5452 0.386 268 1.411 0.7205
## G_07 - G_50 0.4782 0.305 274 1.565 0.6221
##
## DayF = 75, Nutrients = Nutrients:
## contrast estimate SE df t.ratio p.value
## G_48 - G_62 0.1627 0.234 258 0.695 0.9824
## G_48 - G_31 0.8966 0.276 258 3.252 0.0162
## G_48 - G_08 1.2937 0.340 258 3.809 0.0024
## G_48 - G_07 1.0067 0.246 265 4.088 0.0008
## G_48 - G_50 1.9287 0.343 343 5.626 <.0001
## G_62 - G_31 0.7339 0.270 258 2.714 0.0760
## G_62 - G_08 1.1310 0.335 258 3.372 0.0110
## G_62 - G_07 0.8440 0.240 265 3.511 0.0069
## G_62 - G_50 1.7660 0.339 345 5.216 <.0001
## G_31 - G_08 0.3971 0.366 258 1.086 0.8867
## G_31 - G_07 0.1100 0.281 263 0.391 0.9988
## G_31 - G_50 1.0320 0.369 332 2.800 0.0600
## G_08 - G_07 -0.2870 0.344 262 -0.834 0.9609
## G_08 - G_50 0.6350 0.419 315 1.517 0.6536
## G_07 - G_50 0.9220 0.347 345 2.656 0.0870
##
## P value adjustment: tukey method for comparing a family of 6 estimates BW.emmc<-emmeans(LM_5, ~Nutrients|MoteGen|DayF)
contrast(BW.emmc, "tukey")## MoteGen = G_48, DayF = -28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 0.0266 0.284 258 0.094 0.9254
##
## MoteGen = G_62, DayF = -28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 0.5180 0.289 258 1.791 0.0744
##
## MoteGen = G_31, DayF = -28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients -0.3652 0.389 258 -0.940 0.3482
##
## MoteGen = G_08, DayF = -28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients -0.9595 0.612 283 -1.567 0.1182
##
## MoteGen = G_07, DayF = -28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients -0.2935 0.297 258 -0.988 0.3239
##
## MoteGen = G_50, DayF = -28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients -0.2052 0.433 258 -0.474 0.6359
##
## MoteGen = G_48, DayF = 28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 0.4411 0.284 258 1.552 0.1218
##
## MoteGen = G_62, DayF = 28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 0.8156 0.297 275 2.743 0.0065
##
## MoteGen = G_31, DayF = 28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.0631 0.389 258 2.736 0.0067
##
## MoteGen = G_08, DayF = 28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients -0.5379 0.588 258 -0.914 0.3614
##
## MoteGen = G_07, DayF = 28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.1092 0.297 258 3.734 0.0002
##
## MoteGen = G_50, DayF = 28:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.3197 0.433 258 3.048 0.0025
##
## MoteGen = G_48, DayF = 62:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.5485 0.284 258 5.449 <.0001
##
## MoteGen = G_62, DayF = 62:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.4980 0.289 258 5.180 <.0001
##
## MoteGen = G_31, DayF = 62:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 3.1022 0.389 258 7.983 <.0001
##
## MoteGen = G_08, DayF = 62:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.0689 0.588 258 1.817 0.0704
##
## MoteGen = G_07, DayF = 62:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.6474 0.305 274 5.402 <.0001
##
## MoteGen = G_50, DayF = 62:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 2.4613 0.439 266 5.609 <.0001
##
## MoteGen = G_48, DayF = 75:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 2.3396 0.284 258 8.233 <.0001
##
## MoteGen = G_62, DayF = 75:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 2.2289 0.289 258 7.707 <.0001
##
## MoteGen = G_31, DayF = 75:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 3.2005 0.389 258 8.235 <.0001
##
## MoteGen = G_08, DayF = 75:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.5607 0.588 258 2.653 0.0085
##
## MoteGen = G_07, DayF = 75:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 1.7014 0.299 263 5.687 <.0001
##
## MoteGen = G_50, DayF = 75:
## contrast estimate SE df t.ratio p.value
## Ambient - Nutrients 2.8525 0.467 303 6.103 <.0001 BW.emmc<-emmeans(LM_5, ~MoteGen*DayF*Nutrients)
BW_groups_genotypes<-cld(BW.emmc)
#BW_groups_genotypes<-as.data.frame(BW_groups_genotypes[complete.cases(BW_groups_genotypes),])
BW_groups_genotypes<-BW_groups_genotypes[order(BW_groups_genotypes$MoteGen,
BW_groups_genotypes$DayF,
BW_groups_genotypes$Nutrients),]
BW_groups_genotypes #write.csv(BW_groups_genotypes, "Outputs/BW_groups_genotypes.csv", row.names = F)