General project set-up

# 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")
          )

1. Data

# 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

2. Exploratory plots

Plot all

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

Plot change

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

3. LMMs

## 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

4. ANOVAs RAW DATA

## 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)

5. ANOVAs log data

## 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

Packages used

# 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.