# load libraries
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.5.2
## -- Attaching packages ----------------------------------------------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0     v purrr   0.2.5
## v tibble  2.1.1     v dplyr   0.7.7
## v tidyr   0.8.1     v stringr 1.3.1
## v readr   1.1.1     v forcats 0.3.0
## Warning: package 'ggplot2' was built under R version 3.5.2
## Warning: package 'tibble' was built under R version 3.5.3
## -- Conflicts -------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(readxl)
library(broom)
library(knitr)
source('theme_javier.R')

# 01 read and preparae data----- 
ci <- read_excel('data/Length and CI.xlsx')
dotchart(ci$CI)

# remove one outlier
head(arrange(ci, CI))# sample 32 CI<1
## # A tibble: 6 x 14
##   Sample Boat_no Boat_weight Site  Density   Rep Type  Mussel_no Length
##    <dbl>   <dbl>       <dbl> <chr> <chr>   <dbl> <chr>     <dbl>  <dbl>
## 1     32      32       0.897 K     C           4 <NA>          2   70.3
## 2    370     220       0.897 S     C           3 <NA>         10   51.7
## 3     91      91       0.897 K     H           2 <NA>          1   62.7
## 4    223      73       0.897 K     C           2 <NA>          3   78.9
## 5    222      72       0.897 K     C           2 <NA>          2   79.1
## 6    361     211       0.897 S     C           3 <NA>          1   56.8
## # ... with 5 more variables: Flesh_DW_boat <dbl>, Shell_DW <dbl>,
## #   Flesh_DW <dbl>, CI <dbl>, Notes <chr>
ci <- 
  filter(ci,CI>1) %>% 
  mutate(Site = fct_recode(Site, Koutou = "K", Shipwreck = "S"),
         Density = fct_recode(Density, Control = "C", Low = "L", High = "H"),
         Density = fct_relevel(Density, "Control", "Low"),
         Type = if_else(is.na(Type) & Site=="Koutou", "U", Type),
         Type = if_else(is.na(Type), "E", Type))
  
dotchart(ci$CI)

# 02 Summary stats condion index ----------
ci %>% 
  drop_na(CI) %>% 
  group_by(Site,    Density,    Rep, Type) %>% 
  summarise(n = n())
## # A tibble: 37 x 5
## # Groups:   Site, Density, Rep [?]
##    Site   Density   Rep Type      n
##    <fct>  <fct>   <dbl> <chr> <int>
##  1 Koutou Control     1 U        10
##  2 Koutou Control     2 U        10
##  3 Koutou Control     3 U        10
##  4 Koutou Control     4 U         9
##  5 Koutou Control     5 U        10
##  6 Koutou Low         1 U        10
##  7 Koutou Low         2 U        10
##  8 Koutou Low         3 U        10
##  9 Koutou Low         4 U        10
## 10 Koutou Low         5 U        10
## # ... with 27 more rows
ci %>% 
  drop_na(CI) %>% 
  group_by(Site) %>% 
  summarise(mean = mean(CI),
            n = n(),
            sd = sd(CI),
            se = sd/sqrt(n))
## # A tibble: 2 x 5
##   Site       mean     n    sd    se
##   <fct>     <dbl> <int> <dbl> <dbl>
## 1 Koutou     8.65   149  1.67 0.137
## 2 Shipwreck  8.38   229  1.84 0.122
# 03 Boxplot of mussel size-----
ggplot(ci, aes(x = Density, y = Length)) +
  geom_boxplot() +
  facet_grid( Type ~ Site) +
  theme_bw()

# 04 Boxplot of condition index ---------
ci_plot2 <- 
  ci %>% 
  drop_na(CI) %>% 
  mutate(Type = fct_recode(Type, Trasplanted = "E", Undisturbed = "U")) %>% 
ggplot( aes(x = Density, y = CI)) +  
  geom_boxplot() +
  facet_wrap( ~ paste(Site, Type, sep =" ")) +
  theme_javier() +
  labs(y = 'Condition index', x = '')

print(ci_plot2)

ggsave(ci_plot2, 
       filename = 'figures/ci_plots2.tiff',
       device = 'tiff',
       compression = 'lzw',
       width = 7,
       height = 2.5,
       dpi = 600)

# anova condition index ---------
# undistrubed mussels
und_ci <- filter(ci,Type=="U")
m1 <- aov(sqrt(CI)~Site*Density, und_ci)
kable(tidy(m1))
term df sumsq meansq statistic p.value
Site 1 0.0000259 0.0000259 0.0003198 0.9857455
Density 2 0.5736569 0.2868285 3.5448399 0.0303160
Site:Density 2 0.5623719 0.2811859 3.4751055 0.0324443
Residuals 253 20.4713331 0.0809144 NA NA
kable(tidy(TukeyHSD(m1)))
term comparison estimate conf.low conf.high adj.p.value
Site Shipwreck-Koutou -0.0006395 -0.0710608 0.0697818 0.9857455
Density Low-Control -0.0383865 -0.1386382 0.0618652 0.6390819
Density High-Control -0.1149013 -0.2182213 -0.0115812 0.0250606
Density High-Low -0.0765147 -0.1795627 0.0265333 0.1886245
Site:Density Shipwreck:Control-Koutou:Control 0.0527241 -0.1213435 0.2267917 0.9533454
Site:Density Koutou:Low-Koutou:Control 0.0449818 -0.1192228 0.2091865 0.9696023
Site:Density Shipwreck:Low-Koutou:Control -0.0892733 -0.2633408 0.0847943 0.6819471
Site:Density Koutou:High-Koutou:Control -0.1190230 -0.2832277 0.0451816 0.3002917
Site:Density Shipwreck:High-Koutou:Control -0.0447147 -0.2340824 0.1446529 0.9842293
Site:Density Koutou:Low-Shipwreck:Control -0.0077422 -0.1810257 0.1655412 0.9999951
Site:Density Shipwreck:Low-Shipwreck:Control -0.1419974 -0.3246542 0.0406595 0.2266302
Site:Density Koutou:High-Shipwreck:Control -0.1717471 -0.3450306 0.0015363 0.0536005
Site:Density Shipwreck:High-Shipwreck:Control -0.0974388 -0.2947307 0.0998531 0.7159242
Site:Density Shipwreck:Low-Koutou:Low -0.1342551 -0.3075386 0.0390284 0.2300066
Site:Density Koutou:High-Koutou:Low -0.1640049 -0.3273781 -0.0006317 0.0484954
Site:Density Shipwreck:High-Koutou:Low -0.0896966 -0.2783437 0.0989506 0.7476521
Site:Density Koutou:High-Shipwreck:Low -0.0297498 -0.2030333 0.1435337 0.9963886
Site:Density Shipwreck:High-Shipwreck:Low 0.0445585 -0.1527334 0.2418505 0.9871012
Site:Density Shipwreck:High-Koutou:High 0.0743083 -0.1143388 0.2629555 0.8680293
# transplanted mussels
tr_ci <- filter(ci,Type!="U")
m2 <- aov(sqrt(CI)~Density, tr_ci)
kable(tidy(m2))
term df sumsq meansq statistic p.value
Density 2 1.126326 0.5631628 6.282525 0.002567
Residuals 116 10.398190 0.0896396 NA NA
kable(tidy(TukeyHSD(m2)))
term comparison estimate conf.low conf.high adj.p.value
Density Low-Control 0.0582341 -0.0839308 0.2003991 0.5956610
Density High-Control 0.2845258 0.0929568 0.4760948 0.0017393
Density High-Low 0.2262917 0.0347226 0.4178607 0.0161716
# anova size 
m3 <- aov(sqrt(Length)~Site*Density, ci)
tidy(m3)
## # A tibble: 4 x 6
##   term            df  sumsq meansq statistic p.value
##   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
## 1 Site             1   3.96  3.96      4.81   0.0289
## 2 Density          2   2.91  1.46      1.77   0.172 
## 3 Site:Density     2   1.11  0.554     0.672  0.511 
## 4 Residuals      372 306.    0.823    NA     NA
kable(tidy(TukeyHSD(m3)))
term comparison estimate conf.low conf.high adj.p.value
Site Shipwreck-Koutou -0.2094790 -0.3972542 -0.0217038 0.0288794
Density Low-Control -0.2042166 -0.4598497 0.0514164 0.1459540
Density High-Control -0.1091099 -0.3898793 0.1716596 0.6315009
Density High-Low 0.0951068 -0.1852452 0.3754588 0.7043134
Site:Density Shipwreck:Control-Koutou:Control -0.3503983 -0.8118198 0.1110231 0.2517983
Site:Density Koutou:Low-Koutou:Control -0.3731248 -0.8955741 0.1493246 0.3184894
Site:Density Shipwreck:Low-Koutou:Control -0.4617909 -0.9232123 -0.0003695 0.0496832
Site:Density Koutou:High-Koutou:Control -0.2093309 -0.7317802 0.3131185 0.8609330
Site:Density Shipwreck:High-Koutou:Control -0.4006701 -0.9257514 0.1244113 0.2467900
Site:Density Koutou:Low-Shipwreck:Control -0.0227264 -0.4811505 0.4356976 0.9999919
Site:Density Shipwreck:Low-Shipwreck:Control -0.1113926 -0.4988316 0.2760464 0.9630649
Site:Density Koutou:High-Shipwreck:Control 0.1410675 -0.3173566 0.5994915 0.9507829
Site:Density Shipwreck:High-Shipwreck:Control -0.0502717 -0.5116931 0.4111497 0.9996034
Site:Density Shipwreck:Low-Koutou:Low -0.0886662 -0.5470902 0.3697579 0.9937749
Site:Density Koutou:High-Koutou:Low 0.1637939 -0.3560101 0.6835979 0.9456853
Site:Density Shipwreck:High-Koutou:Low -0.0275453 -0.5499946 0.4949040 0.9999890
Site:Density Koutou:High-Shipwreck:Low 0.2524601 -0.2059640 0.7108841 0.6138875
Site:Density Shipwreck:High-Shipwreck:Low 0.0611209 -0.4003005 0.5225423 0.9989735
Site:Density Shipwreck:High-Koutou:High -0.1913392 -0.7137885 0.3311101 0.9008231
# Initial mussel size data---------------
size <- 
  read_csv('data/UW-Caliper 0111_Pyura_4_17.CSV') %>% 
  filter(mm>0) 
## Parsed with column specification:
## cols(
##   Button = col_integer(),
##   Date = col_character(),
##   Site = col_character(),
##   Time = col_time(format = ""),
##   mm = col_double()
## )
summary(size$mm)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   30.50   44.38   47.25   47.99   52.00   66.50
size %>% 
  group_by(Site) %>% 
  summarise(mean = mean(mm),
            n = n(),
            sd = sd(mm),
            se = sd/sqrt(n))
## # A tibble: 2 x 5
##   Site   mean     n    sd    se
##   <chr> <dbl> <int> <dbl> <dbl>
## 1 K      52.4    31  7.51 1.35 
## 2 S      45.6    57  6.07 0.804
size_freq_dist <- 
  ggplot(size, aes(x = mm,  stat(count),fill = Site)) +
  geom_density(alpha = .3) +
  labs(x = 'Mussel size (mm)', y = 'Number of individuals') + 
  theme_javier() +
  theme(legend.position = c(.9,.7))

print(size_freq_dist)

ggsave(size_freq_dist, 
       filename = 'figures/size_freq_dist.tiff',
       device = 'tiff',
       compression = 'lzw',
       width = 4,
       height = 2.5,
       dpi = 600)