# Load all libraries and sources required to run the script
library(tidyverse)
library(ggthemes)
library(ggplot2)
library(lme4)
library(lmerTest)
#library(multcomp)
#library(multcompView)
#library(emmeans)
#library(effects)
ggthe_bw<-theme_bw()+
theme(panel.grid= element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)
# Data
data<-read.csv("Data/Respiration_Adults.csv", header = TRUE)
summary(data)
## Species Core Timepoint Date Tank
## Past:126 A8 : 6 Exposure:110 05/08/17:110 : 6
## Ssid:160 D11 : 6 Initial :110 05/13/17:110 bottom:141
## D17 : 6 Recovery: 66 05/31/17: 66 top :139
## D25 : 6
## D26 : 6
## D27 : 6
## (Other):250
## Trea Treatment Plate Well
## CHS :52 Reef High Shaded:52 Min. :1.000 A1 : 13
## PHS :46 Port High Shaded:46 1st Qu.:2.000 A2 : 13
## PL :42 Low Port :42 Median :3.000 A3 : 13
## CL :41 Low Reef :41 Mean :2.769 A4 : 13
## PH :38 High Port :38 3rd Qu.:4.000 A5 : 13
## NS :37 Control :37 Max. :5.000 A6 : 13
## (Other):30 (Other) :30 (Other):208
## MO2..nmol.min. Corrected.vol.vol FINAL.MO2..nmol.min.
## #VALUE! : 4 : 10 Min. : 0.1829
## : 3 #VALUE! : 4 1st Qu.: 1.7895
## 0.228756043479332: 1 0.669929188315137: 2 Median : 2.6771
## 0.725635126765311: 1 0.738986035620189: 2 Mean : 3.5116
## 0.884822958431114: 1 0.747149255769029: 2 3rd Qu.: 4.1449
## 0.94714543 : 1 0.747787815891798: 2 Max. :19.3248
## (Other) :275 (Other) :264 NA's :24
## Notes X
## :250 :285
## no coral data : 6 NS=No sediment: 1
## no exposure : 6
## dead? : 5
## bad curve? : 4
## starts below 80%: 2
## (Other) : 13
data$Core<-paste(data$Species, data$Core, sep = "-")
data$Tank<-factor(data$Tank, ordered = F)
data$Plate<-factor(data$Plate, ordered = F)
data$Well<-factor(data$Well, ordered = F)
data$Timepoint<-factor(data$Timepoint,
levels=c("Initial", "Exposure","Recovery"))
data<-data[data$Treatment!="Reef High Shaded", ]
data<-data[data$Treatment!="Port High Shaded", ]
data<-data[data$Species!="NA", ]
data$Treatment<-factor(data$Treatment,
# levels=c("Control", "Low Reef","High Reef",
# "Reef High Shaded",
# "Low Port","High Port", "Port High Shaded"))
levels=c("Control", "Low Reef","High Reef",
"Low Port","High Port"))
summary(data)
## Species Core Timepoint Date Tank
## Past: 73 Length:188 Initial :72 05/08/17:72 : 6
## Ssid:115 Class :character Exposure:70 05/13/17:70 bottom:98
## Mode :character Recovery:46 05/31/17:46 top :84
##
##
##
##
## Trea Treatment Plate Well MO2..nmol.min.
## PL :42 Control :37 1:45 C1 : 12 #VALUE! : 2
## CL :41 Low Reef :41 2:46 C3 : 12 : 1
## PH :38 High Reef:24 3:53 B6 : 11 0.228756043479332: 1
## NS :37 Low Port :42 4:21 D3 : 10 0.94714543 : 1
## CH :24 High Port:38 5:23 B4 : 9 0.979184820129077: 1
## : 6 NA's : 6 C2 : 9 1.11640013689544 : 1
## (Other): 0 (Other):125 (Other) :181
## Corrected.vol.vol FINAL.MO2..nmol.min. Notes
## : 6 Min. : 0.1829 :162
## #VALUE! : 4 1st Qu.: 1.8413 no exposure : 6
## 0.669929188315137: 2 Median : 2.7221 dead? : 5
## 0.738986035620189: 2 Mean : 3.6222 no coral data: 5
## 0.749824003596305: 2 3rd Qu.: 4.0938 bad curve? : 2
## 0.751842024117296: 2 Max. :19.3248 bad curve : 1
## (Other) :170 NA's :16 (Other) : 7
## X
## :187
## NS=No sediment: 1
##
##
##
##
##
wide.data<-select(data, c("Species", "Core", "Timepoint", "Tank", "Treatment", "FINAL.MO2..nmol.min."))
wide.data<- reshape(wide.data, idvar = "Core", timevar = "Timepoint", direction = "wide")
head(wide.data)
wide.data$change<-wide.data$FINAL.MO2..nmol.min..Exposure-wide.data$FINAL.MO2..nmol.min..Initial
wide.data<-wide.data[which(wide.data$Species.Exposure!="NA"),]
# YII.Wide$treat_rec<-YII.Wide$YII.after_exposure-YII.Wide$YII.after_recovery
O2<- ggplot(data, aes (Treatment, FINAL.MO2..nmol.min.)) +
geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
stat_summary(fun=mean, geom="point") +
#geom_point(shape=21)+
geom_jitter(alpha=0.5, shape=21)+
theme(legend.position = "bottom")+
scale_y_continuous(# limits = c(0, 15),
# breaks = seq(0,1,0.2),
# expand = c(0.01, 0.01),
name=("O2 [nmol/min]")) +
ggthe_bw +
facet_grid(Species~Timepoint)
O2
delta_O2<- ggplot(wide.data, aes (Treatment.Exposure, change)) +
geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, colour="blue")+
stat_summary(fun=mean, geom="point", colour="blue") +
#geom_point(shape=21)+
geom_point(alpha=0.5, shape=21)+
theme(legend.position = "bottom")+
scale_y_continuous(#limits = c(0, 15),
# breaks = seq(0,1,0.2),
# expand = c(0.01, 0.01),
name=("delta O2 [nmol/min]")) +
ggthe_bw + geom_abline(slope = 0, intercept = 0)+
facet_grid(~Species.Exposure)
delta_O2
log_delta_O2<- ggplot(wide.data, aes (Treatment.Exposure, log(change))) +
geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2, colour="blue")+
stat_summary(fun=mean, geom="point", colour="blue") +
#geom_point(shape=21)+
geom_point(alpha=0.5, shape=21)+
theme(legend.position = "bottom")+
scale_y_continuous(#limits = c(0, 15),
# breaks = seq(0,1,0.2),
# expand = c(0.01, 0.01),
name=("log delta O2 [nmol/min]")) +
ggthe_bw + geom_abline(slope = 0, intercept = 0)+
facet_grid(~Species.Exposure)
log_delta_O2
## model 1: LMER for both species
fit1<-lmerTest::lmer(change ~Treatment.Exposure *
Species.Exposure + (1|Tank.Exposure), data=wide.data)
isSingular(fit1)
## [1] TRUE
anova(fit1)
ranova(fit1)
summary(fit1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: change ~ Treatment.Exposure * Species.Exposure + (1 | Tank.Exposure)
## Data: wide.data
##
## REML criterion at convergence: 248.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1129 -0.2867 0.0899 0.4347 2.0981
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank.Exposure (Intercept) 0.000 0.000
## Residual 6.572 2.564
## Number of obs: 59, groups: Tank.Exposure, 2
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -0.73522 1.14649 49.00000
## Treatment.ExposureLow Reef 1.40018 1.55235 49.00000
## Treatment.ExposureHigh Reef 1.33928 2.14488 49.00000
## Treatment.ExposureLow Port 0.25224 1.46149 49.00000
## Treatment.ExposureHigh Port -0.13807 1.55235 49.00000
## Species.ExposureSsid 0.25165 1.50111 49.00000
## Treatment.ExposureLow Reef:Species.ExposureSsid 0.03967 2.07064 49.00000
## Treatment.ExposureHigh Reef:Species.ExposureSsid -2.65712 2.61798 49.00000
## Treatment.ExposureLow Port:Species.ExposureSsid 0.14972 2.00343 49.00000
## Treatment.ExposureHigh Port:Species.ExposureSsid -0.02515 2.10809 49.00000
## t value Pr(>|t|)
## (Intercept) -0.641 0.524
## Treatment.ExposureLow Reef 0.902 0.371
## Treatment.ExposureHigh Reef 0.624 0.535
## Treatment.ExposureLow Port 0.173 0.864
## Treatment.ExposureHigh Port -0.089 0.929
## Species.ExposureSsid 0.168 0.868
## Treatment.ExposureLow Reef:Species.ExposureSsid 0.019 0.985
## Treatment.ExposureHigh Reef:Species.ExposureSsid -1.015 0.315
## Treatment.ExposureLow Port:Species.ExposureSsid 0.075 0.941
## Treatment.ExposureHigh Port:Species.ExposureSsid -0.012 0.991
##
## Correlation of Fixed Effects:
## (Intr) Tr.ELR Tr.EHR Tr.ELP Tr.EHP Spc.ES T.ELR: T.EHR: T.ELP:
## Trtmnt.ExLR -0.739
## Trtmnt.ExHR -0.535 0.395
## Trtmnt.ExLP -0.784 0.579 0.419
## Trtmnt.ExHP -0.739 0.545 0.395 0.579
## Spcs.ExpsrS -0.764 0.564 0.408 0.599 0.564
## Tr.ELR:S.ES 0.554 -0.750 -0.296 -0.434 -0.409 -0.725
## Tr.EHR:S.ES 0.438 -0.323 -0.819 -0.344 -0.323 -0.573 0.416
## Tr.ELP:S.ES 0.572 -0.423 -0.306 -0.729 -0.423 -0.749 0.543 0.430
## Tr.EHP:S.ES 0.544 -0.402 -0.291 -0.427 -0.736 -0.712 0.516 0.408 0.534
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
par(mfrow=c(2,2))
plot(fit1)
par(mfrow=c(1,1))
#step(fit1)
## model 2: LMER for Past
fit_Past<-lmerTest::lmer(change ~Treatment.Exposure + (1|Tank.Exposure),
data=wide.data[wide.data$Species.Exposure=="Past",])
isSingular(fit_Past)
## [1] TRUE
anova(fit_Past)
ranova(fit_Past)
summary(fit_Past)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: change ~ Treatment.Exposure + (1 | Tank.Exposure)
## Data: wide.data[wide.data$Species.Exposure == "Past", ]
##
## REML criterion at convergence: 106.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4367 -0.1933 0.1040 0.4958 1.3914
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank.Exposure (Intercept) 0.000 0.000
## Residual 5.116 2.262
## Number of obs: 27, groups: Tank.Exposure, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.7352 1.0115 22.0000 -0.727 0.475
## Treatment.ExposureLow Reef 1.4002 1.3696 22.0000 1.022 0.318
## Treatment.ExposureHigh Reef 1.3393 1.8924 22.0000 0.708 0.487
## Treatment.ExposureLow Port 0.2522 1.2894 22.0000 0.196 0.847
## Treatment.ExposureHigh Port -0.1381 1.3696 22.0000 -0.101 0.921
##
## Correlation of Fixed Effects:
## (Intr) Tr.ELR Tr.EHR Tr.ELP
## Trtmnt.ExLR -0.739
## Trtmnt.ExHR -0.535 0.395
## Trtmnt.ExLP -0.784 0.579 0.419
## Trtmnt.ExHP -0.739 0.545 0.395 0.579
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
par(mfrow=c(2,2))
plot(fit_Past)
par(mfrow=c(1,1))
#step(fit_Past)
## model 3: LMER for Ssid
fit_Ssid<-lmerTest::lmer(change ~Treatment.Exposure + (1|Tank.Exposure),
data=wide.data[wide.data$Species.Exposure=="Ssid",])
isSingular(fit_Ssid)
## [1] TRUE
anova(fit_Ssid)
ranova(fit_Ssid)
summary(fit_Ssid)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: change ~ Treatment.Exposure + (1 | Tank.Exposure)
## Data: wide.data[wide.data$Species.Exposure == "Ssid", ]
##
## REML criterion at convergence: 141.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7853 -0.4305 0.0789 0.3817 1.9310
##
## Random effects:
## Groups Name Variance Std.Dev.
## Tank.Exposure (Intercept) 0.000 0.000
## Residual 7.759 2.785
## Number of obs: 32, groups: Tank.Exposure, 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.4836 1.0528 27.0000 -0.459 0.650
## Treatment.ExposureLow Reef 1.4399 1.4889 27.0000 0.967 0.342
## Treatment.ExposureHigh Reef -1.3178 1.6310 27.0000 -0.808 0.426
## Treatment.ExposureLow Port 0.4020 1.4889 27.0000 0.270 0.789
## Treatment.ExposureHigh Port -0.1632 1.5497 27.0000 -0.105 0.917
##
## Correlation of Fixed Effects:
## (Intr) Tr.ELR Tr.EHR Tr.ELP
## Trtmnt.ExLR -0.707
## Trtmnt.ExHR -0.645 0.456
## Trtmnt.ExLP -0.707 0.500 0.456
## Trtmnt.ExHP -0.679 0.480 0.439 0.480
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
plot(fit_Ssid)
#step(fit_Ssid)
# Pairwise comparisons
# Day specific comparisons
# Sw.emmc<-emmeans(fit4, ~Treatment)
# Sw_groups<-cld(Sw.emmc)
# Sw_groups
## model 1: LMER for both species
fit1<-lm(change ~Treatment.Exposure * Species.Exposure,
data=wide.data)
anova(fit1)
summary(fit1)
##
## Call:
## lm(formula = change ~ Treatment.Exposure * Species.Exposure,
## data = wide.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5439 -0.7349 0.2305 1.1144 5.3787
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.73522 1.14649 -0.641
## Treatment.ExposureLow Reef 1.40018 1.55235 0.902
## Treatment.ExposureHigh Reef 1.33928 2.14488 0.624
## Treatment.ExposureLow Port 0.25224 1.46149 0.173
## Treatment.ExposureHigh Port -0.13807 1.55235 -0.089
## Species.ExposureSsid 0.25165 1.50111 0.168
## Treatment.ExposureLow Reef:Species.ExposureSsid 0.03967 2.07064 0.019
## Treatment.ExposureHigh Reef:Species.ExposureSsid -2.65712 2.61798 -1.015
## Treatment.ExposureLow Port:Species.ExposureSsid 0.14972 2.00343 0.075
## Treatment.ExposureHigh Port:Species.ExposureSsid -0.02515 2.10809 -0.012
## Pr(>|t|)
## (Intercept) 0.524
## Treatment.ExposureLow Reef 0.371
## Treatment.ExposureHigh Reef 0.535
## Treatment.ExposureLow Port 0.864
## Treatment.ExposureHigh Port 0.929
## Species.ExposureSsid 0.868
## Treatment.ExposureLow Reef:Species.ExposureSsid 0.985
## Treatment.ExposureHigh Reef:Species.ExposureSsid 0.315
## Treatment.ExposureLow Port:Species.ExposureSsid 0.941
## Treatment.ExposureHigh Port:Species.ExposureSsid 0.991
##
## Residual standard error: 2.564 on 49 degrees of freedom
## (11 observations deleted due to missingness)
## Multiple R-squared: 0.09554, Adjusted R-squared: -0.07058
## F-statistic: 0.5751 on 9 and 49 DF, p-value: 0.8108
par(mfrow=c(2,2))
plot(fit1)
par(mfrow=c(1,1))
#step(fit1)
## model 2: LMER for Past
fit_Past<-lm(change ~Treatment.Exposure,
data=wide.data[wide.data$Species.Exposure=="Past",])
anova(fit_Past)
summary(fit_Past)
##
## Call:
## lm(formula = change ~ Treatment.Exposure, data = wide.data[wide.data$Species.Exposure ==
## "Past", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.5114 -0.4373 0.2353 1.1214 3.1471
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7352 1.0115 -0.727 0.475
## Treatment.ExposureLow Reef 1.4002 1.3696 1.022 0.318
## Treatment.ExposureHigh Reef 1.3393 1.8924 0.708 0.487
## Treatment.ExposureLow Port 0.2522 1.2894 0.196 0.847
## Treatment.ExposureHigh Port -0.1381 1.3696 -0.101 0.921
##
## Residual standard error: 2.262 on 22 degrees of freedom
## (5 observations deleted due to missingness)
## Multiple R-squared: 0.08458, Adjusted R-squared: -0.08186
## F-statistic: 0.5081 on 4 and 22 DF, p-value: 0.7303
par(mfrow=c(2,2))
plot(fit_Past)
par(mfrow=c(1,1))
#step(fit_Past)
## model 3: LMER for Ssid
fit_Ssid<-lm(change ~Treatment.Exposure,
data=wide.data[wide.data$Species.Exposure=="Ssid",])
anova(fit_Ssid)
summary(fit_Ssid)
##
## Call:
## lm(formula = change ~ Treatment.Exposure, data = wide.data[wide.data$Species.Exposure ==
## "Ssid", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.5439 -1.1991 0.2197 1.0633 5.3787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4836 1.0528 -0.459 0.650
## Treatment.ExposureLow Reef 1.4399 1.4889 0.967 0.342
## Treatment.ExposureHigh Reef -1.3178 1.6310 -0.808 0.426
## Treatment.ExposureLow Port 0.4020 1.4889 0.270 0.789
## Treatment.ExposureHigh Port -0.1632 1.5497 -0.105 0.917
##
## Residual standard error: 2.785 on 27 degrees of freedom
## (6 observations deleted due to missingness)
## Multiple R-squared: 0.1013, Adjusted R-squared: -0.0319
## F-statistic: 0.7605 on 4 and 27 DF, p-value: 0.5601
par(mfrow=c(2,2))
plot(fit_Ssid)
par(mfrow=c(1,1))
#step(fit_Ssid)
## model 1: LMER for both species
fit1<-lm(log(change) ~Treatment.Exposure * Species.Exposure,
data=wide.data)
anova(fit1)
#ranova(fit1)
summary(fit1)
##
## Call:
## lm(formula = log(change) ~ Treatment.Exposure * Species.Exposure,
## data = wide.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.51721 -0.44240 0.09378 0.70734 1.66252
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -1.1564 0.6565 -1.761
## Treatment.ExposureLow Reef 1.4986 0.8685 1.725
## Treatment.ExposureHigh Reef 1.3862 1.3131 1.056
## Treatment.ExposureLow Port 0.8482 0.9285 0.914
## Treatment.ExposureHigh Port 1.5354 1.0381 1.479
## Species.ExposureSsid 1.0999 0.9285 1.185
## Treatment.ExposureLow Reef:Species.ExposureSsid -1.1863 1.2283 -0.966
## Treatment.ExposureHigh Reef:Species.ExposureSsid -1.1862 1.6082 -0.738
## Treatment.ExposureLow Port:Species.ExposureSsid -1.4724 1.2714 -1.158
## Treatment.ExposureHigh Port:Species.ExposureSsid -2.1627 1.4680 -1.473
## Pr(>|t|)
## (Intercept) 0.0943 .
## Treatment.ExposureLow Reef 0.1007
## Treatment.ExposureHigh Reef 0.3044
## Treatment.ExposureLow Port 0.3724
## Treatment.ExposureHigh Port 0.1555
## Species.ExposureSsid 0.2508
## Treatment.ExposureLow Reef:Species.ExposureSsid 0.3462
## Treatment.ExposureHigh Reef:Species.ExposureSsid 0.4698
## Treatment.ExposureLow Port:Species.ExposureSsid 0.2612
## Treatment.ExposureHigh Port:Species.ExposureSsid 0.1571
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.137 on 19 degrees of freedom
## (41 observations deleted due to missingness)
## Multiple R-squared: 0.2319, Adjusted R-squared: -0.1319
## F-statistic: 0.6374 on 9 and 19 DF, p-value: 0.752
par(mfrow=c(2,2))
plot(fit1)
par(mfrow=c(1,1))
#step(fit1)
## model 2: LMER for Past
fit_Past<-lm(log(change) ~Treatment.Exposure,
data=wide.data[wide.data$Species.Exposure=="Past",])
anova(fit_Past)
summary(fit_Past)
##
## Call:
## lm(formula = log(change) ~ Treatment.Exposure, data = wide.data[wide.data$Species.Exposure ==
## "Past", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.51721 -0.44240 0.09378 0.70734 1.66252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.1564 0.7481 -1.546 0.161
## Treatment.ExposureLow Reef 1.4986 0.9896 1.514 0.168
## Treatment.ExposureHigh Reef 1.3862 1.4961 0.927 0.381
## Treatment.ExposureLow Port 0.8482 1.0579 0.802 0.446
## Treatment.ExposureHigh Port 1.5354 1.1828 1.298 0.230
##
## Residual standard error: 1.296 on 8 degrees of freedom
## (19 observations deleted due to missingness)
## Multiple R-squared: 0.2627, Adjusted R-squared: -0.1059
## F-statistic: 0.7127 on 4 and 8 DF, p-value: 0.606
par(mfrow=c(2,2))
plot(fit_Past)
par(mfrow=c(1,1))
#step(fit_Past)
## model 3: LMER for Ssid
fit_Ssid<-lm(log(change) ~Treatment.Exposure,
data=wide.data[wide.data$Species.Exposure=="Ssid",])
anova(fit_Ssid)
summary(fit_Ssid)
##
## Call:
## lm(formula = log(change) ~ Treatment.Exposure, data = wide.data[wide.data$Species.Exposure ==
## "Ssid", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.6553 -0.4468 0.1069 0.5068 1.2405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05645 0.58098 -0.097 0.924
## Treatment.ExposureLow Reef 0.31226 0.76857 0.406 0.692
## Treatment.ExposureHigh Reef 0.19998 0.82163 0.243 0.812
## Treatment.ExposureLow Port -0.62421 0.76857 -0.812 0.434
## Treatment.ExposureHigh Port -0.62726 0.91861 -0.683 0.509
##
## Residual standard error: 1.006 on 11 degrees of freedom
## (22 observations deleted due to missingness)
## Multiple R-squared: 0.191, Adjusted R-squared: -0.1032
## F-statistic: 0.6491 on 4 and 11 DF, p-value: 0.6392
par(mfrow=c(2,2))
plot(fit_Ssid)
par(mfrow=c(1,1))
#step(fit_Ssid)
# Pairwise comparisons
# Day specific comparisons
# Sw.emmc<-emmeans(fit4, ~Treatment)
# Sw_groups<-cld(Sw.emmc)
# Sw_groups
# Creates bibliography
#knitr::write_bib(c(.packages()), "packages.bib")
Arnold, Jeffrey B. 2021. Ggthemes: Extra Themes, Scales and Geoms for Ggplot2. https://github.com/jrnold/ggthemes.
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.
Bates, Douglas, and Martin Maechler. 2021. Matrix: Sparse and Dense Matrix Classes and Methods. http://Matrix.R-forge.R-project.org/.
Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2021. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.
Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.
Kuznetsova, Alexandra, Per B. Brockhoff, and Rune H. B. Christensen. 2017. “lmerTest Package: Tests in Linear Mixed Effects Models.” Journal of Statistical Software 82 (13): 1–26. https://doi.org/10.18637/jss.v082.i13.
Kuznetsova, Alexandra, Per Bruun Brockhoff, and Rune Haubo Bojesen Christensen. 2019. LmerTest: Tests in Linear Mixed Effects Models. https://github.com/runehaubo/lmerTestR.
Müller, Kirill, and Hadley Wickham. 2021. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
———. 2019. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.
———. 2021a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.
———. 2021b. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.
———. 2021c. Tidyverse: Easily Install and Load the Tidyverse. https://CRAN.R-project.org/package=tidyverse.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey Dunnington. 2021. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, and Jim Hester. 2020. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.