Load necesary packages

library(tidyverse)
library(dplyr)
library(car)
library(readr)
demo_fiona_NO_T <- read_csv("C:/Users/17872/Documents/tesis 2023/demo fiona NO T.csv",col_types = cols(Period = col_factor(levels = c("Pre-Fiona","Post-Fiona")), Location = col_factor(levels = c("ENO","PTC", "BHL"))))

I have unbalanced design

with(demo_fiona_NO_T, table(Period,Location))
##             Location
## Period       ENO PTC BHL
##   Pre-Fiona   15  15  15
##   Post-Fiona  15  15  10

Calculate summary statistics for Thalasia Cover

demo_fiona_NO_T %>%
  group_by(Location,Period) %>%
  summarise(
    count = n(),
    mean = round(mean(`Th.Cover`, na.rm = TRUE), 2),
    median = round(median(`Th.Cover`, na.rm = TRUE), 2),
    sd = round(sd(Th.Cover, na.rm = TRUE), 2),
    cv = round(sd/mean, 2),
  ) %>%
  ungroup()
## # A tibble: 6 × 7
##   Location Period     count  mean median    sd    cv
##   <fct>    <fct>      <int> <dbl>  <dbl> <dbl> <dbl>
## 1 ENO      Pre-Fiona     15  53.8     80 41.5   0.77
## 2 ENO      Post-Fiona    15  37.2     20 38.7   1.04
## 3 PTC      Pre-Fiona     15  84.1     85  7.14  0.08
## 4 PTC      Post-Fiona    15  61.5     62  9.79  0.16
## 5 BHL      Pre-Fiona     15  50.7     70 33.1   0.65
## 6 BHL      Post-Fiona    10  30       35 20.6   0.69

Box Plot and Line Plot for Thalasia cover

library(ggpubr)
library(cowplot)
TCB1 <- ggboxplot(demo_fiona_NO_T, x = 'Location', y = 'Th.Cover', fill = 'Period')+
    scale_fill_grey()+
  ylab(element_blank())+
   # ylab("Cover (%)")+
    theme_classic()
TCB1

TCB2 <- ggboxplot(demo_fiona_NO_T, x = 'Period', y = 'Th.Cover', fill = 'Period')+
    scale_fill_grey()+
    ylab("Cover(%)")+
    theme_classic()
TCB2

TCL <- ggline(
  demo_fiona_NO_T, x = 'Location', y = 'Th.Cover', color = 'Period',
  add = c('mean_se', 'dotplot'), size = 1,
  palette = c('#257ABA', '#C9B826'))
TCL

Goal is to perform an Unbalanced Two way anova type III sum of squares (SS)

First lets do a two way anova on original data

aov.org <- aov(
  demo_fiona_NO_T$Th.Cover~ Period * Location, data = demo_fiona_NO_T,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(aov.org, type = 'III')
## Anova Table (Type III tests)
## 
## Response: demo_fiona_NO_T$Th.Cover
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)     232190  1 276.7976 < 2.2e-16 ***
## Period            8271  1   9.8598 0.0023766 ** 
## Location         17211  2  10.2588 0.0001095 ***
## Period:Location    140  2   0.0834 0.9200876    
## Residuals        66269 79                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova on log transformed data

aov.log <- aov(
  log(demo_fiona_NO_T$`Th.Cover`+1) ~ Period*Location, data = demo_fiona_NO_T,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(aov.log, type = 'III')
## Anova Table (Type III tests)
## 
## Response: log(demo_fiona_NO_T$Th.Cover + 1)
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)     991.72  1 584.0618 < 2.2e-16 ***
## Period            3.20  1   1.8848 0.1736751    
## Location         32.67  2   9.6200 0.0001823 ***
## Period:Location   0.13  2   0.0385 0.9622285    
## Residuals       134.14 79                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova on rank transformed data

aov.rnk <- aov(
  rank(demo_fiona_NO_T$`Th.Cover`) ~Period *Location, data = demo_fiona_NO_T,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(aov.rnk, type = 'III')
## Anova Table (Type III tests)
## 
## Response: rank(demo_fiona_NO_T$Th.Cover)
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)     146025  1 324.3387 < 2.2e-16 ***
## Period            7848  1  17.4315 7.606e-05 ***
## Location          7055  2   7.8354 0.0007865 ***
## Period:Location    927  2   1.0294 0.3619797    
## Residuals        35568 79                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Anova with rank transformation indicate statistical difference by period and location. No significance in their interaction.

Lets inspect the grand means and the group means

model.tables(aov.org, type = "means", se = TRUE)
## Tables of means
## Grand mean
##          
## 54.21176 
## 
##  Period 
##     Pre-Fiona Post-Fiona
##         62.84       44.5
## rep     45.00       40.0
## 
##  Location 
##       ENO   PTC   BHL
##     46.04 73.31 41.11
## rep 30.00 30.00 25.00
## 
##  Period:Location 
##             Location
## Period       ENO   PTC   BHL  
##   Pre-Fiona  53.80 84.07 50.67
##   rep        15.00 15.00 15.00
##   Post-Fiona 37.20 61.47 30.00
##   rep        15.00 15.00 10.00

Lets inspect the residuals. Rank transformation are not expected to be normally distributed.

res.org = aov.org$resid
res.log = aov.log$resid
res.rnk = aov.rnk$resid
qqnorm(
  res.org, pch = 20, main = "Original Data",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(res.org)

qqnorm(
  res.log, pch = 20, main = "Log-Transformed",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(res.log)

qqnorm(
  res.rnk, pch = 20, main = "Rank-Transformed",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(res.rnk)

#tests not normal 
shapiro.test(res.org)
## 
##  Shapiro-Wilk normality test
## 
## data:  res.org
## W = 0.96608, p-value = 0.02455
leveneTest(demo_fiona_NO_T$`Th.Cover`~ Period*Location, data = demo_fiona_NO_T)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value    Pr(>F)    
## group  5  4.8256 0.0006678 ***
##       79                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Now, let’s inspect the residuals versus the fitted values as a visual check for homogeneity of variances.

plot(aov.org, 1, main = "Original Data")

plot(aov.log, 1, main = "Log-Transformed")

plot(aov.rnk, 1, main = "Rank-Transformed")

Post Hoc test to understand Location

library(multcomp)
summary(glht(aov.rnk, linfct = mcp(Location = "Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = rank(demo_fiona_NO_T$Th.Cover) ~ Period * Location, 
##     data = demo_fiona_NO_T, contrasts = list(Location = "contr.sum", 
##         Period = "contr.sum"))
## 
## Linear Hypotheses:
##                Estimate Std. Error t value Pr(>|t|)    
## PTC - ENO == 0   14.333      5.479   2.616 0.028437 *  
## BHL - ENO == 0   -8.058      5.811  -1.387 0.352463    
## BHL - PTC == 0  -22.392      5.811  -3.853 0.000697 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Thalasia Density Converted to square meter.

demo_fiona_NO_T %>%
  group_by(Location, Period) %>%
  summarise(
    count = n(),
    mean = round(mean(`Th.conv`, na.rm = TRUE), 2),
    median = round(median(`Th.conv`, na.rm = TRUE), 2),
    sd = round(sd(Th.conv, na.rm = TRUE), 2),
    cv = round(sd/mean, 2),
  ) %>%
  ungroup()
## # A tibble: 6 × 7
##   Location Period     count  mean median    sd    cv
##   <fct>    <fct>      <int> <dbl>  <dbl> <dbl> <dbl>
## 1 ENO      Pre-Fiona     15  241.   267.  117.  0.48
## 2 ENO      Post-Fiona    15  183.   200   122.  0.67
## 3 PTC      Pre-Fiona     15  567.   589.  101.  0.18
## 4 PTC      Post-Fiona    15  506.   511.  120.  0.24
## 5 BHL      Pre-Fiona     15  289.   311.  160.  0.55
## 6 BHL      Post-Fiona    10  200    167.  158.  0.79

Box plot

TDB1 <- ggboxplot(
  demo_fiona_NO_T, x = 'Location', y = 'Th.den.conv', fill = 'Period')+
    scale_fill_grey()+
  ylab(element_blank())+
    #ylab("Density (shoots/m^2)")+
    theme_classic()
TDB1

TDB2 <- ggboxplot(
  demo_fiona_NO_T, x = 'Period', y = 'Th.den.conv', fill = 'Period')+
    scale_fill_grey()+
    ylab("Density(shoots/m^2)")+
    theme_classic()
TDB2

TDL <- ggline(
  demo_fiona_NO_T, x = 'Location', y = 'Th.den.conv', color = 'Period',
  add = c('mean_se', 'dotplot'), size = 1,
  palette = c('#257ABA', '#C9B826'))
TDL

 # scale_fill_grey(start = 0.4,end = .9)+
  #ylab("Cover (%)")
  #theme(axis.title.y = element_blank())

Anovas

TDaov.org <- aov(
  demo_fiona_NO_T$Th.den.conv~ Location * Period, data = demo_fiona_NO_T,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(TDaov.org, type = 'III')
## Anova Table (Type III tests)
## 
## Response: demo_fiona_NO_T$Th.den.conv
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)     597.58  1 542.2996 < 2.2e-16 ***
## Location        122.60  2  55.6308 8.354e-16 ***
## Period            6.61  1   5.9954   0.01656 *  
## Location:Period   0.24  2   0.1075   0.89824    
## Residuals        87.05 79                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA procedure on log-transformed data
TDaov.log <- aov(
  log(demo_fiona_NO_T$`Th.den.conv`+1) ~ Location * Period, data = demo_fiona_NO_T,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(TDaov.log, type = 'III')
## Anova Table (Type III tests)
## 
## Response: log(demo_fiona_NO_T$Th.den.conv + 1)
##                  Sum Sq Df  F value    Pr(>F)    
## (Intercept)     118.462  1 918.3373 < 2.2e-16 ***
## Location          9.638  2  37.3559 3.814e-12 ***
## Period            0.759  1   5.8870   0.01753 *  
## Location:Period   0.092  2   0.3558   0.70176    
## Residuals        10.191 79                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA procedure on rank-transformed data
TDaov.rnk <- aov(
  rank(demo_fiona_NO_T$`Th.conv`) ~ Location * Period, data = demo_fiona_NO_T,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(TDaov.rnk, type = 'III')
## Anova Table (Type III tests)
## 
## Response: rank(demo_fiona_NO_T$Th.conv)
##                 Sum Sq Df  F value    Pr(>F)    
## (Intercept)     146373  1 540.1743 < 2.2e-16 ***
## Location         28369  2  52.3455 3.348e-15 ***
## Period            1890  1   6.9741  0.009964 ** 
## Location:Period    155  2   0.2856  0.752336    
## Residuals        21407 79                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.tables(TDaov.org, type = "means", se = TRUE)
## Tables of means
## Grand mean
##          
## 2.744471 
## 
##  Location 
##        ENO    PTC    BHL
##      1.719  4.347  2.052
## rep 30.000 30.000 25.000
## 
##  Period 
##     Pre-Fiona Post-Fiona
##         3.002      2.454
## rep    45.000     40.000
## 
##  Location:Period 
##         Period
## Location Pre-Fiona Post-Fiona
##      ENO  1.956     1.482    
##      rep 15.000    15.000    
##      PTC  4.596     4.098    
##      rep 15.000    15.000    
##      BHL  2.340     1.620    
##      rep 15.000    10.000

Normality and variance is normal. Best Choice is to use original data.

###norm

TDres.org = TDaov.org$resid
TDres.log = TDaov.log$resid
TDres.rnk = TDaov.rnk$resid

shapiro.test(TDres.org)
## 
##  Shapiro-Wilk normality test
## 
## data:  TDres.org
## W = 0.97833, p-value = 0.163
leveneTest(Th.den.conv ~ Period*Location, data = demo_fiona_NO_T)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  5  0.4155 0.8367
##       79
qqnorm(
  TDres.org, pch = 20, main = "Original Data",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(TDres.org)

#
qqnorm(
  TDres.log, pch = 20, main = "Log-Transformed",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(TDres.log)

#
qqnorm(
  TDres.rnk, pch = 20, main = "Rank-Transformed",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(TDres.rnk)

#
plot(TDaov.org, 1, main = "Original Data")

plot(TDaov.log, 1, main = "Log-Transformed")

plot(TDaov.rnk, 1, main = "Rank-Transformed")

Post HOC test identifies PTC as significantly different to ENO and BHL. BHL and ENO are nos significantly different.

library(multcomp)
summary(glht(TDaov.org, linfct = mcp(Location = "Tukey")))
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: aov(formula = demo_fiona_NO_T$Th.den.conv ~ Location * Period, 
##     data = demo_fiona_NO_T, contrasts = list(Location = "contr.sum", 
##         Period = "contr.sum"))
## 
## Linear Hypotheses:
##                Estimate Std. Error t value Pr(>|t|)    
## PTC - ENO == 0   2.6280     0.2710   9.696   <1e-04 ***
## BHL - ENO == 0   0.2610     0.2875   0.908    0.637    
## BHL - PTC == 0  -2.3670     0.2875  -8.234   <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Thalasia canopy length

Tcanopy <- read_csv("C:/Users/17872/Documents/tesis 2023/cannopy plots thalassia.csv",col_types = cols(Location = col_factor(levels = c("ENO","PTC", "BHL")), Period = col_factor(levels = c("Pre-Fiona","Post-Fiona"))))

head(Tcanopy)
## # A tibble: 6 × 3
##   Location Height Period   
##   <fct>     <dbl> <fct>    
## 1 BHL          23 Pre-Fiona
## 2 BHL          34 Pre-Fiona
## 3 BHL          15 Pre-Fiona
## 4 BHL          20 Pre-Fiona
## 5 BHL          22 Pre-Fiona
## 6 BHL          15 Pre-Fiona

T canopy length by period

Tcanopy %>%
  group_by(Period) %>%
  summarise(
    count = n(),
    mean = round(mean(`Height`, na.rm = TRUE), 2),
    median = round(median(`Height`, na.rm = TRUE), 2),
    sd = round(sd(Height, na.rm = TRUE), 2),
    cv = round(sd/mean, 2),
  ) %>%
  ungroup()
## # A tibble: 2 × 6
##   Period     count  mean median    sd    cv
##   <fct>      <int> <dbl>  <dbl> <dbl> <dbl>
## 1 Pre-Fiona    159  18.6     19  7.54  0.41
## 2 Post-Fiona   159  12.2     11  5.37  0.44

T canopy length by period and location

Tcanopy %>%
  group_by(Location,Period) %>%
  summarise(
    count = n(),
    mean = round(mean(`Height`, na.rm = TRUE), 2),
    median = round(median(`Height`, na.rm = TRUE), 2),
    sd = round(sd(Height, na.rm = TRUE), 2),
    cv = round(sd/mean, 2),
  ) %>%
  ungroup()
## # A tibble: 6 × 7
##   Location Period     count  mean median    sd    cv
##   <fct>    <fct>      <int> <dbl>  <dbl> <dbl> <dbl>
## 1 ENO      Pre-Fiona     46  26.6     27  4.26  0.16
## 2 ENO      Post-Fiona    46  15.4     14  5.73  0.37
## 3 PTC      Pre-Fiona     64  16.0     16  5.26  0.33
## 4 PTC      Post-Fiona    64  10.3     10  2.88  0.28
## 5 BHL      Pre-Fiona     49  14.5     13  6.72  0.46
## 6 BHL      Post-Fiona    49  11.7     10  6.2   0.53

Graficas

TcanB1 <- ggboxplot(
  Tcanopy, x = 'Location', y = 'Height', fill = 'Period')+
    scale_fill_grey()+
  ylab(element_blank())+
   # ylab("Canopy Height (cm)")+
    theme_classic()
TcanB1

TcanB2 <- ggboxplot(
  Tcanopy, x = 'Period', y = 'Height', fill = 'Period')+
    scale_fill_grey()+
    ylab("Canopy(cm)")+
    theme_classic()
TcanB2

TcanL <- ggline(
  Tcanopy, x = 'Location', y = 'Height', color = 'Period',
  add = c('mean_se', 'dotplot'), size = 1,
  palette = c('#257ABA', '#C9B826'))
TcanL

Anovas

Tcanaov.org <- aov(
  Tcanopy$Height~ Period * Location, data = Tcanopy,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(Tcanaov.org, type = 'III')
## Anova Table (Type III tests)
## 
## Response: Tcanopy$Height
##                 Sum Sq  Df  F value    Pr(>F)    
## (Intercept)      77385   1 2815.528 < 2.2e-16 ***
## Period            3351   1  121.932 < 2.2e-16 ***
## Location          4085   2   74.319 < 2.2e-16 ***
## Period:Location    874   2   15.897 2.665e-07 ***
## Residuals         8575 312                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA procedure on log-transformed data
Tcanaov.log <- aov(
  log(Height) ~ Period * Location, data = Tcanopy,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(Tcanaov.log, type = 'III')
## Anova Table (Type III tests)
## 
## Response: log(Height)
##                  Sum Sq  Df    F value    Pr(>F)    
## (Intercept)     2179.85   1 16563.1007 < 2.2e-16 ***
## Period            12.76   1    96.9472 < 2.2e-16 ***
## Location          15.18   2    57.6693 < 2.2e-16 ***
## Period:Location    1.93   2     7.3481 0.0007615 ***
## Residuals         41.06 312                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA procedure on rank-transformed data
Tcanaov.rnk <- aov(
  rank(Height) ~ Period * Location, data = Tcanopy,
  contrasts = list(
    Location = 'contr.sum',
    Period = 'contr.sum'
  )
)
Anova(Tcanaov.rnk, type = 'III')
## Anova Table (Type III tests)
## 
## Response: rank(Height)
##                  Sum Sq  Df   F value  Pr(>F)    
## (Intercept)     8306515   1 1717.8708 < 2e-16 ***
## Period           546368   1  112.9943 < 2e-16 ***
## Location         556393   2   57.5339 < 2e-16 ***
## Period:Location   50279   2    5.1991 0.00601 ** 
## Residuals       1508631 312                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model.tables(Tcanaov.org, type = "means", se = TRUE)
## Tables of means
## Grand mean
##          
## 15.41509 
## 
##  Period 
##     Pre-Fiona Post-Fiona
##         18.61      12.22
## rep    159.00     159.00
## 
##  Location 
##       ENO    PTC   BHL
##     21.03  13.16 13.09
## rep 92.00 128.00 98.00
## 
##  Period:Location 
##             Location
## Period       ENO   PTC   BHL  
##   Pre-Fiona  26.65 15.98 14.49
##   rep        46.00 64.00 49.00
##   Post-Fiona 15.41 10.34 11.68
##   rep        46.00 64.00 49.00

Nomality checks. T canopy data is not normal. Best choice is Rank Transformed

###norm

Tcanres.org = Tcanaov.org$resid
Tcanres.log = Tcanaov.log$resid
Tcanres.rnk = Tcanaov.rnk$resid

shapiro.test(Tcanres.org)
## 
##  Shapiro-Wilk normality test
## 
## data:  Tcanres.org
## W = 0.94613, p-value = 2.233e-09
leveneTest(Height ~ Period*Location, data = Tcanopy)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value    Pr(>F)    
## group   5  4.7746 0.0003211 ***
##       312                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqnorm(
  Tcanres.org, pch = 20, main = "Original Data",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(Tcanres.org)

#
qqnorm(
  TDres.log, pch = 20, main = "Log-Transformed",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(TDres.log)

#
qqnorm(
  Tcanres.rnk, pch = 20, main = "Rank-Transformed",
  cex.lab = 1, cex.axis = 0.7, cex.main = 1
)
qqline(Tcanres.rnk)

#
plot(Tcanaov.org, 1, main = "Original Data")

plot(Tcanaov.log, 1, main = "Log-Transformed")

plot(Tcanaov.rnk, 1, main = "Rank-Transformed")

Lets combine the three main box plots. Two columns. 1 By period 2 by period an location.

comb<- ggarrange(TCB2,TCB1,TDB2,TDB1,TcanB2,TcanB1,
                        ncol = 2, nrow = 3,common.legend = TRUE, legend="bottom")
comb

  annotate_figure(comb, top = text_grob("Thalasia testudinum", face = "italic", size = 14))