A-Maine-Tank-2022_Bivalve-Height_V7

this script will be the same as 6 except with changes/additions based on discussion with Brittany to be more like Cameron 2022 methods. Would also like finalize plots.

Set up

read in relevant files

I took the dfs Diana sent me and read them in here.

remove dead and lost

We only had 66 removed. This seems a bit high to me and I wonder if the scallops that were removed weren’t marked in coloration? or if something else strange happened. But rolling with this for now. From Diana there should be 66 deaths across these groups.

## [1] 477  15

## [1] 406  15

## [1] 35

look at data

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

strange specimens

based on our conversation with Brittany, it seems safe to assume that measurements where growth is below 0 are mostly a reflection of sampling error and not actual loss in shell max height over the experiment. There I would like to make the assumption that any - values here are really a reflection of 0 growth. In some sense we are flooring percent change in height to 0. At this point all specimens marked as dead or lost have been removed so this shouldn’t effect those and all other specimens with missing data are NA for change in height so they wont be included in this either.

Another approach to this would be to add the value of the most negative number to the whole column so that the lowest % change is 0 % change, but this makes less sense to me in the context of this data set.

height.all.w.negatives<-height.all
height.all$`percentage change height` <- replace(height.all$`percentage change height`, which(height.all$`percentage change height` < 0), 0)

height.all$prop.change.height<-height.all$`percentage change height`/100
height.all %>% 
  ggplot(aes(prop.change.height, group = species, fill=species)) +
    geom_boxplot()

## Warning: Removed 35 rows containing non-finite values (`stat_boxplot()`).

height.all %>% 
  ggplot(aes(x=temp.average, y=prop.change.height,color = species, fill=species)) +
    geom_point()+
    geom_smooth(method = loess) #sweet no longer any below 0 

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 35 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 35 rows containing missing values (`geom_point()`).

look for trends

Scallop doesn’t have all temp pH groups covered, will not be able to model I don’t think… Also mercenaria 12 temp seems to almost be showing two trends, I wonder if this is explained by tank?

In this version it is updated so each species shows geom_smooth with each temp/pH grouping its own line. pH scale is switched so light color = more acidic. Remove SE coloring as plot was too busy.Noting here that it could be nice to have this plot in Altair in the future, so you could hover over points, but thats a rabbit hole for another day.

## Loading required package: ggpp

## 
## Attaching package: 'ggpp'

## The following object is masked from 'package:ggplot2':
## 
##     annotate

## Loading required package: viridisLite

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

summary table

summary_by_species <- describeBy(height.all, group = "species")

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

summary_by_species

## 
##  Descriptive statistics by group 
## species: ad.arctica
##                          vars  n  mean   sd median trimmed  mad   min   max
## tank*                       1 10  8.10 5.47   6.50    8.00 7.41  1.00 16.00
## Sample_ID_20220415*         2 10  5.50 3.03   5.50    5.50 3.71  1.00 10.00
## Max_heigh_mm                3 10 53.34 3.94  55.02   53.74 2.49 46.61 56.81
## Max_height_mm               4 10 53.81 3.87  55.31   54.17 2.62 47.46 57.21
## percentage change height    5 10  0.93 0.59   0.75    0.93 0.51  0.00  1.82
## ph*                         6 10  2.60 1.26   3.00    2.62 1.48  1.00  4.00
## temp*                       7 10  2.30 0.67   2.00    2.38 0.74  1.00  3.00
## died*                       8  0   NaN   NA     NA     NaN   NA   Inf  -Inf
## species*                    9 10  1.00 0.00   1.00    1.00 0.00  1.00  1.00
## temp.average               10 10  9.88 1.94   9.18   10.08 2.22  6.14 12.06
## temp.stdev                 11 10  0.49 0.13   0.44    0.47 0.14  0.33  0.75
## pH.average.YSI             12 10  7.72 0.26   7.82    7.73 0.26  7.36  8.00
## pH.stdev.YSI               13 10  0.07 0.01   0.07    0.07 0.01  0.04  0.08
## pH.average                 14 10  7.72 0.25   7.81    7.72 0.28  7.39  8.01
## pH.stdev                   15 10  0.06 0.01   0.06    0.06 0.01  0.05  0.08
## prop.change.height         16 10  0.01 0.01   0.01    0.01 0.01  0.00  0.02
##                          range  skew kurtosis   se
## tank*                    15.00  0.19    -1.74 1.73
## Sample_ID_20220415*       9.00  0.00    -1.56 0.96
## Max_heigh_mm             10.20 -0.65    -1.40 1.25
## Max_height_mm             9.75 -0.58    -1.53 1.22
## percentage change height  1.82  0.17    -1.41 0.19
## ph*                       3.00 -0.20    -1.77 0.40
## temp*                     2.00 -0.31    -1.14 0.21
## died*                     -Inf    NA       NA   NA
## species*                  0.00   NaN      NaN 0.00
## temp.average              5.92 -0.32    -1.12 0.61
## temp.stdev                0.42  0.59    -1.05 0.04
## pH.average.YSI            0.64 -0.28    -1.79 0.08
## pH.stdev.YSI              0.04 -0.77    -0.16 0.00
## pH.average                0.62 -0.27    -1.79 0.08
## pH.stdev                  0.03  0.44    -1.08 0.00
## prop.change.height        0.02  0.17    -1.41 0.00
## ------------------------------------------------------------ 
## species: juv.arctica
##                          vars  n  mean    sd median trimmed   mad   min    max
## tank*                       1 66  8.35  4.65   8.00    8.31  5.93  1.00  16.00
## Sample_ID_20220415*         2 66 33.50 19.20  33.50   33.50 24.46  1.00  66.00
## Max_heigh_mm                3 63 12.94  2.12  12.66   12.97  2.15  6.64  17.35
## Max_height_mm               4 66 21.60  3.04  21.59   21.68  3.22 14.87  28.67
## percentage change height    5 63 68.27 20.16  65.30   67.28 21.07 30.89 127.86
## ph*                         6 66  2.59  1.14   3.00    2.61  1.48  1.00   4.00
## temp*                       7 66  1.98  0.69   2.00    1.98  0.00  1.00   3.00
## died*                       8  0   NaN    NA     NA     NaN    NA   Inf   -Inf
## species*                    9 66  2.00  0.00   2.00    2.00  0.00  2.00   2.00
## temp.average               10 66  9.01  1.94   9.03    8.99  0.58  6.14  12.06
## temp.stdev                 11 66  0.51  0.16   0.44    0.49  0.16  0.33   0.82
## pH.average.YSI             12 66  7.72  0.23   7.82    7.73  0.28  7.36   8.02
## pH.stdev.YSI               13 66  0.07  0.02   0.07    0.07  0.01  0.04   0.14
## pH.average                 14 66  7.74  0.23   7.81    7.74  0.28  7.39   8.05
## pH.stdev                   15 66  0.07  0.02   0.06    0.06  0.01  0.05   0.14
## prop.change.height         16 63  0.68  0.20   0.65    0.67  0.21  0.31   1.28
##                          range  skew kurtosis   se
## tank*                    15.00  0.06    -1.29 0.57
## Sample_ID_20220415*      65.00  0.00    -1.25 2.36
## Max_heigh_mm             10.71 -0.18    -0.20 0.27
## Max_height_mm            13.80 -0.16    -0.28 0.37
## percentage change height 96.97  0.52    -0.25 2.54
## ph*                       3.00 -0.10    -1.43 0.14
## temp*                     2.00  0.02    -0.93 0.08
## died*                     -Inf    NA       NA   NA
## species*                  0.00   NaN      NaN 0.00
## temp.average              5.92  0.08    -0.93 0.24
## temp.stdev                0.49  0.68    -0.99 0.02
## pH.average.YSI            0.66 -0.23    -1.41 0.03
## pH.stdev.YSI              0.10  2.07     5.18 0.00
## pH.average                0.66 -0.25    -1.35 0.03
## pH.stdev                  0.09  2.22     5.12 0.00
## prop.change.height        0.97  0.52    -0.25 0.03
## ------------------------------------------------------------ 
## species: mercenaria
##                          vars   n  mean    sd median trimmed   mad  min    max
## tank*                       1 115  8.07  4.57   8.00    8.01  5.93 1.00  16.00
## Sample_ID_20220415*         2 115 58.00 33.34  58.00   58.00 43.00 1.00 115.00
## Max_heigh_mm                3 115 12.02  1.20  11.86   11.97  1.10 9.60  15.21
## Max_height_mm               4 115 14.13  2.78  13.59   13.93  2.67 9.69  20.75
## percentage change height    5 115 18.43 19.76  13.02   15.79 19.31 0.00  71.26
## ph*                         6 115  2.50  1.13   2.00    2.51  1.48 1.00   4.00
## temp*                       7 115  1.93  0.73   2.00    1.91  1.48 1.00   3.00
## died*                       8   0   NaN    NA     NA     NaN    NA  Inf   -Inf
## species*                    9 115  3.00  0.00   3.00    3.00  0.00 3.00   3.00
## temp.average               10 115  8.88  2.07   9.03    8.83  3.75 6.14  12.06
## temp.stdev                 11 115  0.53  0.16   0.44    0.52  0.16 0.33   0.82
## pH.average.YSI             12 115  7.71  0.23   7.63    7.71  0.33 7.36   8.02
## pH.stdev.YSI               13 115  0.07  0.02   0.07    0.07  0.01 0.04   0.14
## pH.average                 14 115  7.72  0.23   7.67    7.72  0.39 7.39   8.05
## pH.stdev                   15 115  0.07  0.02   0.06    0.07  0.01 0.05   0.14
## prop.change.height         16 115  0.18  0.20   0.13    0.16  0.19 0.00   0.71
##                           range  skew kurtosis   se
## tank*                     15.00  0.05    -1.19 0.43
## Sample_ID_20220415*      114.00  0.00    -1.23 3.11
## Max_heigh_mm               5.61  0.43    -0.56 0.11
## Max_height_mm             11.06  0.58    -0.45 0.26
## percentage change height  71.26  0.97    -0.08 1.84
## ph*                        3.00  0.04    -1.41 0.11
## temp*                      2.00  0.11    -1.16 0.07
## died*                      -Inf    NA       NA   NA
## species*                   0.00   NaN      NaN 0.00
## temp.average               5.92  0.13    -1.14 0.19
## temp.stdev                 0.49  0.47    -1.22 0.02
## pH.average.YSI             0.66 -0.09    -1.43 0.02
## pH.stdev.YSI               0.10  1.95     4.38 0.00
## pH.average                 0.66 -0.12    -1.38 0.02
## pH.stdev                   0.09  2.05     4.19 0.00
## prop.change.height         0.71  0.97    -0.08 0.02
## ------------------------------------------------------------ 
## species: mya
##                          vars   n  mean    sd median trimmed   mad   min    max
## tank*                       1 150  8.40  4.56   8.00    8.39  5.93  1.00  16.00
## Sample_ID_20220415*         2 150 75.50 43.45  75.50   75.50 55.60  1.00 150.00
## Max_heigh_mm                3 150 15.32  1.22  15.29   15.28  1.14 12.07  20.94
## Max_height_mm               4 150 19.17  2.16  19.27   19.25  1.91 12.13  23.93
## percentage change height    5 150 25.43 12.29  25.87   25.42 14.31  0.00  54.09
## ph*                         6 150  2.50  1.12   2.50    2.50  0.74  1.00   4.00
## temp*                       7 150  1.95  0.69   2.00    1.93  0.00  1.00   3.00
## died*                       8   0   NaN    NA     NA     NaN    NA   Inf   -Inf
## species*                    9 150  4.00  0.00   4.00    4.00  0.00  4.00   4.00
## temp.average               10 150  8.92  1.95   9.03    8.88  0.58  6.14  12.06
## temp.stdev                 11 150  0.52  0.16   0.44    0.51  0.16  0.33   0.82
## pH.average.YSI             12 150  7.71  0.23   7.72    7.71  0.23  7.36   8.02
## pH.stdev.YSI               13 150  0.07  0.02   0.07    0.07  0.01  0.04   0.14
## pH.average                 14 150  7.72  0.23   7.74    7.72  0.22  7.39   8.05
## pH.stdev                   15 150  0.07  0.02   0.06    0.06  0.01  0.05   0.14
## prop.change.height         16 150  0.25  0.12   0.26    0.25  0.14  0.00   0.54
##                           range  skew kurtosis   se
## tank*                     15.00  0.00    -1.23 0.37
## Sample_ID_20220415*      149.00  0.00    -1.22 3.55
## Max_heigh_mm               8.87  0.73     2.51 0.10
## Max_height_mm             11.80 -0.36     0.09 0.18
## percentage change height  54.09  0.01    -0.66 1.00
## ph*                        3.00  0.00    -1.37 0.09
## temp*                      2.00  0.07    -0.93 0.06
## died*                      -Inf    NA       NA   NA
## species*                   0.00   NaN      NaN 0.00
## temp.average               5.92  0.09    -0.93 0.16
## temp.stdev                 0.49  0.59    -1.12 0.01
## pH.average.YSI             0.66 -0.14    -1.39 0.02
## pH.stdev.YSI               0.10  2.01     4.75 0.00
## pH.average                 0.66 -0.15    -1.33 0.02
## pH.stdev                   0.09  2.17     4.73 0.00
## prop.change.height         0.54  0.01    -0.66 0.01
## ------------------------------------------------------------ 
## species: scallop
##                          vars  n  mean    sd median trimmed   mad   min    max
## tank*                       1 65  8.51  4.33   9.00    8.49  5.93  1.00  16.00
## Sample_ID_20220415*         2 65 33.00 18.91  33.00   33.00 23.72  1.00  65.00
## Max_heigh_mm                3 65 32.42  3.33  32.81   32.36  3.68 25.01  40.44
## Max_height_mm               4 33 54.40  7.46  55.57   54.61  7.16 34.81  69.96
## percentage change height    5 33 67.52 18.09  68.72   68.34 17.47 28.26 103.84
## ph*                         6 65  2.46  1.12   2.00    2.45  1.48  1.00   4.00
## temp*                       7 65  1.88  0.72   2.00    1.85  1.48  1.00   3.00
## died*                       8  0   NaN    NA     NA     NaN    NA   Inf   -Inf
## species*                    9 65  5.00  0.00   5.00    5.00  0.00  5.00   5.00
## temp.average               10 65  8.71  1.99   9.00    8.63  3.71  6.14  12.06
## temp.stdev                 11 65  0.54  0.17   0.55    0.54  0.25  0.33   0.82
## pH.average.YSI             12 65  7.70  0.23   7.63    7.71  0.30  7.36   8.02
## pH.stdev.YSI               13 65  0.07  0.02   0.07    0.07  0.01  0.04   0.14
## pH.average                 14 65  7.72  0.23   7.67    7.72  0.33  7.39   8.05
## pH.stdev                   15 65  0.07  0.02   0.06    0.06  0.01  0.05   0.14
## prop.change.height         16 33  0.68  0.18   0.69    0.68  0.17  0.28   1.04
##                          range  skew kurtosis   se
## tank*                    15.00  0.01    -1.24 0.54
## Sample_ID_20220415*      64.00  0.00    -1.26 2.35
## Max_heigh_mm             15.43  0.13    -0.35 0.41
## Max_height_mm            35.15 -0.34     0.03 1.30
## percentage change height 75.58 -0.36    -0.35 3.15
## ph*                       3.00  0.03    -1.39 0.14
## temp*                     2.00  0.18    -1.09 0.09
## died*                     -Inf    NA       NA   NA
## species*                  0.00   NaN      NaN 0.00
## temp.average              5.92  0.19    -1.08 0.25
## temp.stdev                0.49  0.31    -1.51 0.02
## pH.average.YSI            0.66 -0.10    -1.41 0.03
## pH.stdev.YSI              0.10  2.20     5.61 0.00
## pH.average                0.66 -0.13    -1.36 0.03
## pH.stdev                  0.09  2.18     4.96 0.00
## prop.change.height        0.76 -0.36    -0.35 0.03

#summary_by_species <- height.all %>%
 # group_by(species) %>%
 # summarise(
  #  mean_height = mean(prop.change.height),
  #  median_height = median(prop.change.height),
  #  min_height = min(prop.change.height),
  #  max_height = max(prop.change.height),
   # N=n ()
  ##)

#write.csv(summary_by_species,'02_output/01_modified_data/height_summary_by_species.csv', row.names=F)

add a pH where minimum pH is 0 (to avoid weird intercepts later on)

Not totally sure if this is necessary for temperature too - did not do it here.

## [1] 7.39

Individual Species Groups

plot scallop and adult arctica growth

#species.labs <- c(  `ad.arctica` = "Adult Arctica islandica",
                    #`scallop` = "Placopecten magellanicus",
                   # )
#species.labs <- c("Arctica islandica","Mercenaria mercenaria", "Mya arenaria")

height.scal.plot.1 <- ggplot(scallop,aes(x=ph, y=prop.change.height, color=temp)) +
  geom_smooth(aes(group=interaction(species,temp)),method="lm", se=FALSE)+
  geom_point(aes(x=ph, y=prop.change.height,color=temp), size=3, position = position_dodge(width = 0.5)) + 
  xlab("pH Treatment")+
  ylab("% ∆ in Shell Height")+
  labs(x="pH Treatment", y="% ∆ in Shell Height", color = "Temp (C)")+
  scale_color_viridis(discrete = TRUE, option="plasma")+ 
  #stat_poly_eq(formula = x ~ y) +
  facet_wrap(~species, ncol = 1,
             labeller = as_labeller(species.labs))+
    theme_classic()+
    theme(
      strip.text.x = element_text(
        size = 12, color = "black", face = "italic"),
     panel.background = element_blank(),
     panel.border = element_rect(colour = "black", fill=NA, size=1),
     #strip.background = element_blank(),
     #panel.grid.major = element_blank(),
     #panel.grid.minor = element_blank(),
     legend.position = "bottom",
     
     axis.text.x = element_text(size= 10, family = "Helvetica"),
     axis.text.y = element_text(size= 10, family = "Helvetica"),
     axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
     axis.ticks.y = element_blank(),
     axis.title.y = element_blank(),
     
    )
height.scal.plot.1

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 32 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 32 rows containing missing values (`geom_point()`).

height.scal.plot.2 <- ggplot(scallop,aes(x=temp, y=prop.change.height, color=ph)) +
  geom_smooth(aes(group=interaction(species,ph)),method="lm", se=FALSE)+
  geom_point(aes(x=temp, y=prop.change.height,color=ph), size=3, position = position_dodge(width = 0.5)) + 
  labs(x="Temp Treatment (C)", y="% ∆ in Shell Height", color = "pH")+
  scale_color_viridis(discrete = TRUE, direction = -1, option="cividis")+ 
  #stat_poly_eq(formula = x ~ y) +
  facet_wrap(~species, ncol = 1,
             labeller = as_labeller(species.labs))+
    theme_classic()+
    theme(
      strip.text.x = element_text(
        size = 12, color = "black", face = "italic"),
     panel.background = element_blank(),
     panel.border = element_rect(colour = "black", fill=NA, size=1),
     #strip.background = element_blank(),
     #panel.grid.major = element_blank(),
     #panel.grid.minor = element_blank(),
     legend.position = "bottom",
     axis.text.x = element_text(size= 10, family = "Helvetica"),
     axis.text.y = element_text(size= 10, family = "Helvetica"),
     axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
     axis.ticks.y = element_blank(),
     axis.title.y = element_blank(),
     
    )
height.scal.plot.2

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 32 rows containing non-finite values (`stat_smooth()`).
## Removed 32 rows containing missing values (`geom_point()`).

height.ad.plot.1 <- ggplot(ad.arctica,aes(x=ph, y=prop.change.height, color=temp)) +
  geom_smooth(aes(group=interaction(species,temp)),method="lm", se=FALSE)+
  geom_point(aes(x=ph, y=prop.change.height,color=temp), size=3, position = position_dodge(width = 0.5)) + 
  xlab("pH Treatment")+
  ylab("% ∆ in Shell Height")+
  labs(x="pH Treatment", y="% ∆ in Shell Height", color = "Temp (C)")+
  scale_color_viridis(discrete = TRUE, option="plasma")+ 
  #stat_poly_eq(formula = x ~ y) +
  facet_wrap(~species, ncol = 1,
             labeller = as_labeller(species.labs))+
    theme_classic()+
    theme(
      strip.text.x = element_text(
        size = 12, color = "black", face = "italic"),
     panel.background = element_blank(),
     panel.border = element_rect(colour = "black", fill=NA, size=1),
     #strip.background = element_blank(),
     #panel.grid.major = element_blank(),
     #panel.grid.minor = element_blank(),
     legend.position = "bottom",
     
     axis.text.x = element_text(size= 10, family = "Helvetica"),
     axis.text.y = element_text(size= 10, family = "Helvetica"),
     axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
     axis.ticks.y = element_blank(),
     axis.title.y = element_blank(),
     
    )
height.ad.plot.1

## `geom_smooth()` using formula = 'y ~ x'

height.ad.plot.2 <- ggplot(ad.arctica,aes(x=temp, y=prop.change.height, color=ph)) +
  geom_smooth(aes(group=interaction(species,ph)),method="lm", se=FALSE)+
  geom_point(aes(x=temp, y=prop.change.height,color=ph), size=3, position = position_dodge(width = 0.5)) + 
  labs(x="Temp Treatment (C)", y="% ∆ in Shell Height", color = "pH")+
  scale_color_viridis(discrete = TRUE, direction = -1, option="cividis")+ 
  #stat_poly_eq(formula = x ~ y) +
  facet_wrap(~species, ncol = 1,
             labeller = as_labeller(species.labs))+
    theme_classic()+
    theme(
      strip.text.x = element_text(
        size = 12, color = "black", face = "italic"),
     panel.background = element_blank(),
     panel.border = element_rect(colour = "black", fill=NA, size=1),
     #strip.background = element_blank(),
     #panel.grid.major = element_blank(),
     #panel.grid.minor = element_blank(),
     legend.position = "bottom",
     axis.text.x = element_text(size= 10, family = "Helvetica"),
     axis.text.y = element_text(size= 10, family = "Helvetica"),
     axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
     axis.ticks.y = element_blank(),
     axis.title.y = element_blank(),
     
    )
height.ad.plot.2

## `geom_smooth()` using formula = 'y ~ x'

height_plot_combined<-ggarrange(height.scal.plot.1, height.scal.plot.2, height.ad.plot.1, height.ad.plot.2, ncol = 2, nrow = 2)

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 32 rows containing non-finite values (`stat_smooth()`).
## Removed 32 rows containing missing values (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 32 rows containing non-finite values (`stat_smooth()`).
## Removed 32 rows containing missing values (`geom_point()`).

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

height_plot_combined

ggsave("02_output/02_plots/height.scal.ad.plot.png", height_plot_combined, width = 6, height = 8, dpi = 300, bg="white")

mercenaria

check structure

Non normal, majorly skewed many at 0

## 'data.frame':    115 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ Sample_ID_20220415      : chr  "b47" "b39" "b38" "b43" ...
##  $ Max_heigh_mm            : num  11.4 11.4 12.4 11.9 10.2 ...
##  $ Max_height_mm           : num  17.4 17.2 18.3 17.2 15 ...
##  $ percentage change height: num  53.6 49.9 47.5 44.7 46 ...
##  $ ph                      : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 4 ...
##  $ temp                    : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 2 ...
##  $ died                    : chr  NA NA NA NA ...
##  $ species                 : Factor w/ 5 levels "ad.arctica","juv.arctica",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ temp.average            : num  12.1 12.1 12.1 12.1 12.1 ...
##  $ temp.stdev              : num  0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.63 ...
##  $ pH.average.YSI          : num  7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 8 ...
##  $ pH.stdev.YSI            : num  0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.07 ...
##  $ pH.average              : num  7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 8 ...
##  $ pH.stdev                : num  0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.05 ...
##  $ prop.change.height      : num  0.536 0.499 0.475 0.447 0.46 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mercenaria$prop.change.height
## W = 0.85068, p-value = 2.09e-09

##    vars   n mean  sd median trimmed  mad min  max range skew kurtosis   se
## X1    1 115 0.18 0.2   0.13    0.16 0.19   0 0.71  0.71 0.97    -0.08 0.02

## [1] 0

## [1] 115  17

linear model

‘normal’ linear mixed effect model

first we will try this and see how it looks as per Brittany’s suggestion. This will not specify a family, most similar to what I did way back when with CCA.

best fit includes temp, observed vs expected is off QQ line, residual vs predicted is wacky. We can see if other models look better or worse.

check the normality assumptions of full model

looks a bit funky/kurtose, see difference in variance in temperature groups. Also some strange patterning in residuals.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = mercenaria, na.action = na.fail)
simulationOutput1<-simulateResiduals(height.1)
plot(simulationOutput1)

shapiro.test(resid(height.1)) # ! non normal

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(height.1)
## W = 0.8921, p-value = 1.298e-07

plot(height.1)

#plot(top_model1)
testZeroInflation(height.1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided

plotResiduals(height.1, mercenaria$temp)

plotResiduals(height.1, mercenaria$ph)

#plotResiduals(top_model1, mercenaria$temp)

make models and likelihood ratio test

when talking with Diana we decided that 1) there is evidence that temperature effects growth (see long lived arctica records etc.) but that we were not really testing that here. We were more interested in whether there was a relationship with pH. I will therefore test against a null model of … 1) temperature + tank 2) tank alone are we asking vs full model? or main effects?

full.m<-lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = mercenaria, na.action = na.fail)
main.m<-lmer(prop.change.height~temp.average+pH.normalized+(1|tank), data = mercenaria, na.action = na.fail)

full.v.main<-anova(full.m, main.m, test="Chisq")

## refitting model(s) with ML (instead of REML)

full.v.main # no dif - drops AIC, drop interaction 

## Data: mercenaria
## Models:
## main.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
## full.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##        npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## main.m    5 -249.72 -235.99 129.86  -259.72                     
## full.m    6 -248.26 -231.79 130.13  -260.26 0.5481  1     0.4591

summary(main.m)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##    Data: mercenaria
## 
## REML criterion at convergence: -243.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7026 -0.2263  0.0282  0.5222  2.7092 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.006423 0.08014 
##  Residual             0.004525 0.06727 
## Number of obs: 115, groups:  tank, 16
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   -0.60480    0.10563 11.76761  -5.725 0.000103 ***
## temp.average   0.08646    0.01061 11.84902   8.147 3.39e-06 ***
## pH.normalized  0.04371    0.09307 11.98456   0.470 0.646995    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.937       
## pH.normalzd -0.375  0.094

confint(main.m)

## Computing profile confidence intervals ...

##                     2.5 %      97.5 %
## .sig01         0.04756438  0.10969799
## .sigma         0.05891418  0.07801634
## (Intercept)   -0.80252479 -0.40725191
## temp.average   0.06661002  0.10631516
## pH.normalized -0.13048709  0.21771482

drop1(main.m,test="Chisq")

## Single term deletions using Satterthwaite's method:
## 
## Model:
## prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##                 Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## temp.average  0.300366 0.300366     1 11.849 66.3752 3.389e-06 ***
## pH.normalized 0.000998 0.000998     1 11.985  0.2206     0.647    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(main.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	-0.60	-0.81 – -0.40	<0.001
temp average	0.09	0.07 – 0.11	<0.001
pH normalized	0.04	-0.14 – 0.23	0.639
Random Effects
σ2	0.00
τ00 tank	0.01
ICC	0.59
N tank	16
Observations	115
Marginal R2 / Conditional R2	0.744 / 0.894

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more. cant make a plot with both since pH isnt in model.

plot(predictorEffects(full.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

plot(predictorEffects(main.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

Another way to plot

values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=full.m, xlevels=list(temp.average=c(6, 9, 12), pH.normalized=c(0, 0.2, 0.4, 0.6)))

## Warning in term == terms: longer object length is not a multiple of shorter
## object length

## Warning in term == names: longer object length is not a multiple of shorter
## object length

## NOTE: temp.averagepH.normalized is not a high-order term in the model

effect_interact.both<-as.data.frame(effect_interact.both)
effect_interact.both

##    temp.average pH.normalized         fit         se       lower       upper
## 1             6           0.0 -0.05806880 0.06882930 -0.19445865 0.078321057
## 2             9           0.0  0.17201470 0.03815258  0.09641282 0.247616585
## 3            12           0.0  0.40209820 0.06818960  0.26697595 0.537220440
## 4             6           0.2 -0.06592969 0.04614203 -0.15736321 0.025503824
## 5             9           0.2  0.18227502 0.02488040  0.13297286 0.231577189
## 6            12           0.2  0.43047974 0.04325440  0.34476826 0.516191225
## 7             6           0.4 -0.07379059 0.04047314 -0.15399082 0.006409638
## 8             9           0.4  0.19253535 0.02284410  0.14726824 0.237802454
## 9            12           0.4  0.45886129 0.04215161  0.37533505 0.542387524
## 10            6           0.6 -0.08165149 0.05713525 -0.19486879 0.031565815
## 11            9           0.6  0.20279567 0.03411855  0.13518748 0.270403861
## 12           12           0.6  0.48724283 0.06608531  0.35629038 0.618195276

original_values <- c(0, 0.2, 0.4, 0.6)
custom_labels_x <- c("7.4", "7.6", "7.8", "8.0")

effect_temp <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=pH.normalized, y=fit), color="black") +
  geom_line(data = effect_interact.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_interact.both, aes(x=pH.normalized, ymin=lower, ymax=upper, group = temp.average, fill=temp.average), alpha= 0.3) +
  geom_point(data=mercenaria, mapping = aes(x=pH.normalized, y=prop.change.height, color=temp.average), size=3) +
  labs(x="pH", y="Proportional Change in Height", fill = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = "none", breaks = c(6, 9, 12))+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)+
  #scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica"),
    axis.title.y = element_text(size=15,family = "Helvetica"),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp

effect_temp_swap <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=temp.average, y=fit), color="black") +
   geom_line(data = effect_interact.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
   geom_ribbon(data= effect_interact.both, aes(x=temp.average, ymin=lower, ymax=upper, group = pH.normalized, fill=pH.normalized), alpha= 0.3) +
   geom_point(data=mercenaria, mapping = aes(x=temp.average, y=prop.change.height, color=pH.normalized), size=3) +
   labs(x="Temperature (C)", y="Proportional Change in Height", color = "pH")+
  theme_classic()+
  scale_color_viridis( option = "cividis", breaks = original_values, labels = custom_labels_x, direction = -1)+
  scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
    axis.title.y = element_text(size=15,family = "Helvetica",hjust = 0.5),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp_swap

mercenaria.model.plot<-ggarrange(effect_temp_swap,
             labels = c("A"),
             font.label = list(size = 20, color = "black", face = "bold", family = "Helvetica"),
             ncol = 1, nrow = 1)

mercenaria.model.plot<-annotate_figure(mercenaria.model.plot, top = text_grob("M. mercenaria", 
               color = "black", face = 3, size = 20, family = "Helvetica"))

ggsave("02_output/02_plots/height_mercenaria_linear.png", mercenaria.model.plot, width = 6, height = 5, dpi = 300, bg = 'white')

as expected, temperature was a significant predictor in the full model and a full model of all was not significantly different than a model which only included temperature. pH had no influence.

another plot, ph

based on Tredennick which bases it off of https://biologyforfun.wordpress.com/2015/06/17/confidence-intervals-for-prediction-in-glmms/, should keep in mind that confidence or prediction intervals might underestimate uncertainty around our estimates.

mercenaria_vec <- seq(min(mercenaria$temp.average),max(mercenaria$temp.average))

upper_ph <- quantile(mercenaria$pH.normalized, 0.9) 
lower_ph <- quantile(mercenaria$pH.normalized, 0.1)

newdat <- tibble(temp.average = mercenaria_vec,
               Tank = mercenaria$tank[1],
               pH.normalized = lower_ph,
               ph_level = "low") %>%
  bind_rows(
    tibble(temp.average = mercenaria_vec,
           Tank = mercenaria$tank[1],
           pH.normalized = upper_ph,
           ph_level = "high")
  ) %>%
  as.data.frame()

m <- full.m
allout <- tibble()  # empty object
for(do_level in unique(newdat$ph_level)) {
  tmpdat <- newdat %>%
    filter(ph_level == do_level)
  
  mm <- model.matrix(~pH.normalized*temp.average, tmpdat)
  tmpdat$y <- mm%*%fixef(m) 
  # predict(m, newdat, re.form=NA) would give the same results
  pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm))
  
  outdat <- data.frame(tmpdat,
                       lower_ci = tmpdat$y-1.96*sqrt(pvar1),
                       upper_ci = tmpdat$y+1.96*sqrt(pvar1))
  allout <- bind_rows(allout, outdat)
}

mercenaria.int.ci.plot<-ggplot(allout, aes(x = temp.average, y = y,
                   color = ph_level, fill = ph_level)) +
  geom_hline(aes(yintercept = 0), linetype = 3) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), color = NA, alpha = 0.25) +
  geom_line() +
  scale_color_manual(values = rev(c("coral","dodgerblue")), 
                     labels = c("High pH", "Low pH"),
                     name = NULL) +
  scale_fill_manual(values = rev(c("coral","dodgerblue")),
                    labels = c("High pH", "Low pH"),
                    name = NULL) +
  #scale_y_continuous(breaks = c(0,1)) +
  #scale_x_continuous(breaks = c(6,9,12)) +
  labs(x = "Mean Temperature", y = "Proportional Change in Height") +
  theme_classic() +
  theme(legend.position = c(0.8, 0.85))
mercenaria.int.ci.plot

ggsave("02_output/02_plots/height_mercenaria_int.ci.plot.png", height=3.2, width=4, units = "in")

okay so this plot shows there is no interaction, but is otherwise pretty unhelpful, I am also a bit concerned that the axis scales don’t make sense. I have a feeling all our plots would look like this/the same bc it was only siginificant in one case, therfore I am not sure how helpful it is.

mya

check structure

This distribution looks pretty normal.

## 'data.frame':    150 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ Sample_ID_20220415      : chr  "r39" "r40" "r37" "r38" ...
##  $ Max_heigh_mm            : num  13.8 16.5 17.1 13.8 14.3 ...
##  $ Max_height_mm           : num  19.9 21.6 22.2 15.6 15 ...
##  $ percentage change height: num  44.62 31.04 30.25 12.68 4.76 ...
##  $ ph                      : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 4 ...
##  $ temp                    : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 2 ...
##  $ died                    : chr  NA NA NA NA ...
##  $ species                 : Factor w/ 5 levels "ad.arctica","juv.arctica",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ temp.average            : num  12.1 12.1 12.1 12.1 12.1 ...
##  $ temp.stdev              : num  0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.63 ...
##  $ pH.average.YSI          : num  7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 8 ...
##  $ pH.stdev.YSI            : num  0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.07 ...
##  $ pH.average              : num  7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 8 ...
##  $ pH.stdev                : num  0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.05 ...
##  $ prop.change.height      : num  0.4462 0.3104 0.3025 0.1268 0.0476 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mya$prop.change.height
## W = 0.98778, p-value = 0.2124

##    vars   n mean   sd median trimmed  mad min  max range skew kurtosis   se
## X1    1 150 0.25 0.12   0.26    0.25 0.14   0 0.54  0.54 0.01    -0.66 0.01

linear model

‘normal’ linear mixed effect model

first we will try this and see how it looks as per Brittany’s suggestion. This will not specify a family, most similar to what I did way back when with CCA.

best fit includes temp, observed vs expected is off QQ line, residual vs predicted is wacky. We can see if other models look better or worse.

check the normality assumptions of full model

looks a bit funky/kurtose, see difference in variance in temperature groups. Also some strange patterning in residuals.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = mya, na.action = na.fail)
simulationOutput1<-simulateResiduals(height.1)
plot(simulationOutput1)

shapiro.test(resid(height.1)) # ! non normal

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(height.1)
## W = 0.98946, p-value = 0.3221

plot(height.1)

#plot(top_model1)
testZeroInflation(height.1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = Inf, p-value < 2.2e-16
## alternative hypothesis: two.sided

plotResiduals(height.1, mya$temp)

plotResiduals(height.1, mya$ph)

#plotResiduals(top_model1, mya$temp)

make models and likelihood ratio test

when talking with Diana we decided that 1) there is evidence that temperature effects growth (see long lived arctica records etc.) but that we were not really testing that here. We were more interested in whether there was a relationship with pH. I will therefore test against a null model of … 1) temperature + tank 2) tank alone are we asking vs full model? or main effects?

full.m<-lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = mya, na.action = na.fail)
main.m<-lmer(prop.change.height~temp.average+pH.normalized+(1|tank), data = mya, na.action = na.fail)

full.v.main<-anova(full.m, main.m, test="Chisq")

## refitting model(s) with ML (instead of REML)

full.v.main # no dif - drops AIC, drop interaction 

## Data: mya
## Models:
## main.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
## full.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##        npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)
## main.m    5 -215.91 -200.86 112.96  -225.91                    
## full.m    6 -213.91 -195.85 112.96  -225.91 0.001  1      0.975

summary(main.m)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##    Data: mya
## 
## REML criterion at convergence: -206.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.29429 -0.58316  0.00947  0.67847  2.23781 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.001121 0.03348 
##  Residual             0.012438 0.11152 
## Number of obs: 150, groups:  tank, 16
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)    0.072040   0.062757 13.237098   1.148  0.27132   
## temp.average   0.021370   0.006351 13.690657   3.365  0.00475 **
## pH.normalized -0.026586   0.054867 13.028571  -0.485  0.63604   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.937       
## pH.normalzd -0.377  0.097

confint(main.m)

## Computing profile confidence intervals ...

##                      2.5 %     97.5 %
## .sig01         0.000000000 0.05515095
## .sigma         0.099386885 0.12629520
## (Intercept)   -0.045817179 0.19027507
## temp.average   0.009407677 0.03329598
## pH.normalized -0.129821572 0.07651674

drop1(main.m,test="Chisq")

## Single term deletions using Satterthwaite's method:
## 
## Model:
## prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##                Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## temp.average  0.14083 0.14083     1 13.691 11.3226 0.004754 **
## pH.normalized 0.00292 0.00292     1 13.029  0.2348 0.636040   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(main.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.07	-0.05 – 0.20	0.253
temp average	0.02	0.01 – 0.03	0.001
pH normalized	-0.03	-0.14 – 0.08	0.629
Random Effects
σ2	0.01
τ00 tank	0.00
ICC	0.08
N tank	16
Observations	150
Marginal R2 / Conditional R2	0.118 / 0.191

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more. cant make a plot with both since pH isnt in model.

plot(predictorEffects(full.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

plot(predictorEffects(main.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

Another way to plot

values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=full.m, xlevels=list(temp.average=c(6, 9, 12), pH.normalized=c(0, 0.2, 0.4, 0.6)))

## Warning in term == terms: longer object length is not a multiple of shorter
## object length

## Warning in term == names: longer object length is not a multiple of shorter
## object length

## NOTE: temp.averagepH.normalized is not a high-order term in the model

effect_interact.both<-as.data.frame(effect_interact.both)
effect_interact.both

##    temp.average pH.normalized       fit         se     lower     upper
## 1             6           0.0 0.1993400 0.04092660 0.1184549 0.2802251
## 2             9           0.0 0.2644614 0.02287098 0.2192605 0.3096624
## 3            12           0.0 0.3295828 0.04203060 0.2465158 0.4126498
## 4             6           0.2 0.1946269 0.02746439 0.1403478 0.2489060
## 5             9           0.2 0.2590335 0.01484378 0.2296970 0.2883699
## 6            12           0.2 0.3234401 0.02652063 0.2710262 0.3758540
## 7             6           0.4 0.1899138 0.02404763 0.1423874 0.2374402
## 8             9           0.4 0.2536056 0.01362957 0.2266688 0.2805423
## 9            12           0.4 0.3172973 0.02582734 0.2662536 0.3683411
## 10            6           0.6 0.1852007 0.03386701 0.1182678 0.2521336
## 11            9           0.6 0.2481776 0.02047836 0.2077053 0.2886499
## 12           12           0.6 0.3111545 0.04071480 0.2306880 0.3916211

original_values <- c(0, 0.2, 0.4, 0.6)
custom_labels_x <- c("7.4", "7.6", "7.8", "8.0")

effect_temp <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=pH.normalized, y=fit), color="black") +
  geom_line(data = effect_interact.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_interact.both, aes(x=pH.normalized, ymin=lower, ymax=upper, group = temp.average, fill=temp.average), alpha= 0.3) +
  geom_point(data=mya, mapping = aes(x=pH.normalized, y=prop.change.height, color=temp.average), size=3) +
  labs(x="pH", y="Proportional Change in Height", fill = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = "none", breaks = c(6, 9, 12))+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)+
  #scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica"),
    axis.title.y = element_text(size=15,family = "Helvetica"),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp

effect_temp_swap <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=temp.average, y=fit), color="black") +
   geom_line(data = effect_interact.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
   geom_ribbon(data= effect_interact.both, aes(x=temp.average, ymin=lower, ymax=upper, group = pH.normalized, fill=pH.normalized), alpha= 0.3) +
   geom_point(data=mya, mapping = aes(x=temp.average, y=prop.change.height, color=pH.normalized), size=3) +
   labs(x="Temperature (C)", y="Proportional Change in Height", color = "pH")+
  theme_classic()+
  scale_color_viridis( option = "cividis", breaks = original_values, labels = custom_labels_x, direction = -1)+
  scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
    axis.title.y = element_text(size=15,family = "Helvetica",hjust = 0.5),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp_swap

mya.model.plot<-ggarrange(effect_temp_swap,
             labels = c("B"),
             font.label = list(size = 20, color = "black", face = "bold", family = "Helvetica"),
             ncol = 1, nrow = 1)

mya.model.plot<-annotate_figure(mya.model.plot, top = text_grob("M. arenaria", 
               color = "black", face = 3, size = 20, family = "Helvetica"))

ggsave("02_output/02_plots/height_mya_linear.png", mya.model.plot, width = 6, height = 5, dpi = 300, bg = 'white')

as expected, temperature was a significant predictor in the full model and a full model of all was not significantly different than a model which only included temperature. pH had no influence.

another plot, ph

based on Tredennick which bases it off of https://biologyforfun.wordpress.com/2015/06/17/confidence-intervals-for-prediction-in-glmms/, should keep in mind that confidence or prediction intervals might underestimate uncertainty around our estimates.

mya_vec <- seq(min(mya$temp.average),max(mya$temp.average))

upper_ph <- quantile(mya$pH.normalized, 0.9) 
lower_ph <- quantile(mya$pH.normalized, 0.1)

newdat <- tibble(temp.average = mya_vec,
               Tank = mya$tank[1],
               pH.normalized = lower_ph,
               ph_level = "low") %>%
  bind_rows(
    tibble(temp.average = mya_vec,
           Tank = mya$tank[1],
           pH.normalized = upper_ph,
           ph_level = "high")
  ) %>%
  as.data.frame()

m <- full.m
allout <- tibble()  # empty object
for(do_level in unique(newdat$ph_level)) {
  tmpdat <- newdat %>%
    filter(ph_level == do_level)
  
  mm <- model.matrix(~pH.normalized*temp.average, tmpdat)
  tmpdat$y <- mm%*%fixef(m) 
  # predict(m, newdat, re.form=NA) would give the same results
  pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm))
  
  outdat <- data.frame(tmpdat,
                       lower_ci = tmpdat$y-1.96*sqrt(pvar1),
                       upper_ci = tmpdat$y+1.96*sqrt(pvar1))
  allout <- bind_rows(allout, outdat)
}

mya.int.ci.plot<-ggplot(allout, aes(x = temp.average, y = y,
                   color = ph_level, fill = ph_level)) +
  geom_hline(aes(yintercept = 0), linetype = 3) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), color = NA, alpha = 0.25) +
  geom_line() +
  scale_color_manual(values = rev(c("coral","dodgerblue")), 
                     labels = c("High pH", "Low pH"),
                     name = NULL) +
  scale_fill_manual(values = rev(c("coral","dodgerblue")),
                    labels = c("High pH", "Low pH"),
                    name = NULL) +
  #scale_y_continuous(breaks = c(0,1)) +
  #scale_x_continuous(breaks = c(6,9,12)) +
  labs(x = "Mean Temperature", y = "Proportional Change in Height") +
  theme_classic() +
  theme(legend.position = c(0.8, 0.85))
mya.int.ci.plot

ggsave("02_output/02_plots/height_mya_int.ci.plot.png", height=3.2, width=4, units = "in")

okay so this plot shows there is no interaction, but is otherwise pretty unhelpful, I am also a bit concerned that the axis scales don’t make sense. I have a feeling all our plots would look like this/the same bc it was only significant in one case, therefore I am not sure how helpful it is.

juvenile arctica

check structure

This distribution looks pretty normal, although it looks like there is one outlier above 1 (g69). Could possibly make max growth = 1 but I think I will leave for now.

## 'data.frame':    66 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Sample_ID_20220415      : chr  "g32" "g29" "g30" "g31" ...
##  $ Max_heigh_mm            : num  15 13.3 14.2 16.2 14.4 ...
##  $ Max_height_mm           : num  27.7 25.5 23.9 28.7 23.3 ...
##  $ percentage change height: num  84.4 92.3 67.6 77.3 61.5 ...
##  $ ph                      : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 4 4 4 4 4 4 ...
##  $ temp                    : Factor w/ 3 levels "6","9","12": 3 3 3 3 2 2 2 2 1 1 ...
##  $ died                    : chr  NA NA NA NA ...
##  $ species                 : Factor w/ 5 levels "ad.arctica","juv.arctica",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ temp.average            : num  12.06 12.06 12.06 12.06 9.19 ...
##  $ temp.stdev              : num  0.42 0.42 0.42 0.42 0.63 0.63 0.63 0.63 0.78 0.78 ...
##  $ pH.average.YSI          : num  7.57 7.57 7.57 7.57 8 8 8 8 8.02 8.02 ...
##  $ pH.stdev.YSI            : num  0.04 0.04 0.04 0.04 0.07 0.07 0.07 0.07 0.06 0.06 ...
##  $ pH.average              : num  7.59 7.59 7.59 7.59 8 8 8 8 8.05 8.05 ...
##  $ pH.stdev                : num  0.07 0.07 0.07 0.07 0.05 0.05 0.05 0.05 0.05 0.05 ...
##  $ prop.change.height      : num  0.844 0.923 0.676 0.773 0.615 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.61 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  juv.arctica$prop.change.height
## W = 0.96809, p-value = 0.1012

##    vars  n mean  sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 63 0.68 0.2   0.65    0.67 0.21 0.31 1.28  0.97 0.52    -0.25 0.03

## [1] 0

## [1] 4

linear model

‘normal’ linear mixed effect model

first we will try this and see how it looks as per Brittany’s suggestion. This will not specify a family, most similar to what I did way back when with CCA.

best fit includes temp, observed vs expected is off QQ line, residual vs predicted is wacky. We can see if other models look better or worse.

check the normality assumptions of full model

looks a bit funky/kurtose, see difference in variance in temperature groups. Also some strange patterning in residuals.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = juv.arctica, na.action = na.fail)
simulationOutput1<-simulateResiduals(height.1)
plot(simulationOutput1)

shapiro.test(resid(height.1)) # ! non normal

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(height.1)
## W = 0.94464, p-value = 0.006836

plot(height.1)

#plot(top_model1)
testZeroInflation(height.1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(height.1, juv.arctica$temp)

plotResiduals(height.1, juv.arctica$ph)

#plotResiduals(top_model1, juv.arctica$temp)

make models and likelihood ratio test

when talking with Diana we decided that 1) there is evidence that temperature effects growth (see long lived arctica records etc.) but that we were not really testing that here. We were more interested in whether there was a relationship with pH. I will therefore test against a null model of … 1) temperature + tank 2) tank alone are we asking vs full model? or main effects?

full.m<-lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = juv.arctica, na.action = na.fail)
main.m<-lmer(prop.change.height~temp.average+pH.normalized+(1|tank), data = juv.arctica, na.action = na.fail)

full.v.main<-anova(full.m, main.m, test="Chisq")

## refitting model(s) with ML (instead of REML)

full.v.main # no dif - drops AIC, drop interaction 

## Data: juv.arctica
## Models:
## main.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
## full.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##        npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## main.m    5 -44.763 -34.047 27.381  -54.763                     
## full.m    6 -44.256 -31.397 28.128  -56.256 1.4933  1     0.2217

summary(main.m)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: -38.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8291 -0.6244  0.0138  0.6179  3.7785 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.001295 0.03598 
##  Residual             0.024664 0.15705 
## Number of obs: 63, groups:  tank, 16
## 
## Fixed effects:
##                Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)    0.116780   0.109861 14.584785   1.063    0.305    
## temp.average   0.063107   0.011044 14.788089   5.714 4.33e-05 ***
## pH.normalized -0.008728   0.097179 13.835307  -0.090    0.930    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.930       
## pH.normalzd -0.380  0.077

confint(main.m)

## Computing profile confidence intervals ...

##                     2.5 %     97.5 %
## .sig01         0.00000000 0.08530094
## .sigma         0.12999150 0.18858436
## (Intercept)   -0.09030652 0.32445215
## temp.average   0.04228462 0.08395225
## pH.normalized -0.19677250 0.17751909

drop1(main.m,test="Chisq")

## Single term deletions using Satterthwaite's method:
## 
## Model:
## prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##                Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## temp.average  0.80529 0.80529     1 14.788 32.6497 4.332e-05 ***
## pH.normalized 0.00020 0.00020     1 13.835  0.0081    0.9297    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(main.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.12	-0.10 – 0.34	0.292
temp average	0.06	0.04 – 0.09	<0.001
pH normalized	-0.01	-0.20 – 0.19	0.929
Random Effects
σ2	0.02
τ00 tank	0.00
ICC	0.05
N tank	16
Observations	63
Marginal R2 / Conditional R2	0.378 / 0.409

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more. cant make a plot with both since pH isnt in model.

plot(predictorEffects(full.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

plot(predictorEffects(main.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

Another way to plot

values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=full.m, xlevels=list(temp.average=c(6, 9, 12), pH.normalized=c(0, 0.2, 0.4, 0.6)))

## Warning in term == terms: longer object length is not a multiple of shorter
## object length

## Warning in term == names: longer object length is not a multiple of shorter
## object length

## NOTE: temp.averagepH.normalized is not a high-order term in the model

effect_interact.both<-as.data.frame(effect_interact.both)
effect_interact.both

##    temp.average pH.normalized       fit         se     lower     upper
## 1             6           0.0 0.5438041 0.06922427 0.4052866 0.6823215
## 2             9           0.0 0.6797093 0.04002847 0.5996125 0.7598061
## 3            12           0.0 0.8156146 0.07224533 0.6710520 0.9601772
## 4             6           0.2 0.5137870 0.04631160 0.4211177 0.6064563
## 5             9           0.2 0.6819577 0.02587944 0.6301730 0.7337423
## 6            12           0.2 0.8501283 0.04523869 0.7596059 0.9406507
## 7             6           0.4 0.4837700 0.03978918 0.4041520 0.5633879
## 8             9           0.4 0.6842060 0.02204619 0.6400917 0.7283203
## 9            12           0.4 0.8846420 0.04139455 0.8018117 0.9674723
## 10            6           0.6 0.4537529 0.05574265 0.3422121 0.5652937
## 11            9           0.6 0.6864543 0.03242137 0.6215793 0.7513294
## 12           12           0.6 0.9191557 0.06496380 0.7891635 1.0491480

original_values <- c(0, 0.2, 0.4, 0.6)
custom_labels_x <- c("7.4", "7.6", "7.8", "8.0")

effect_temp <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=pH.normalized, y=fit), color="black") +
  geom_line(data = effect_interact.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_interact.both, aes(x=pH.normalized, ymin=lower, ymax=upper, group = temp.average, fill=temp.average), alpha= 0.3) +
  geom_point(data=juv.arctica, mapping = aes(x=pH.normalized, y=prop.change.height, color=temp.average), size=3) +
  labs(x="pH", y="Proportional Change in Height", fill = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = "none", breaks = c(6, 9, 12))+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)+
  #scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica"),
    axis.title.y = element_text(size=15,family = "Helvetica"),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp

effect_temp_swap <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=temp.average, y=fit), color="black") +
   geom_line(data = effect_interact.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
   geom_ribbon(data= effect_interact.both, aes(x=temp.average, ymin=lower, ymax=upper, group = pH.normalized, fill=pH.normalized), alpha= 0.3) +
   geom_point(data=juv.arctica, mapping = aes(x=temp.average, y=prop.change.height, color=pH.normalized), size=3) +
   labs(x="Temperature (C)", y="Proportional Change in Height", color = "pH")+
  theme_classic()+
  scale_color_viridis( option = "cividis", breaks = original_values, labels = custom_labels_x, direction = -1)+
  scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
    axis.title.y = element_text(size=15,family = "Helvetica",hjust = 0.5),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp_swap

juv.arctica.model.plot<-ggarrange(effect_temp_swap,
             labels = c("C"),
             font.label = list(size = 20, color = "black", face = "bold", family = "Helvetica"),
             ncol = 1, nrow = 1)

juv.arctica.model.plot<-annotate_figure(juv.arctica.model.plot, top = text_grob("A. islandica", 
               color = "black", face = 3, size = 20, family = "Helvetica"))

ggsave("02_output/02_plots/height_juv.arctica_linear.png", juv.arctica.model.plot, width = 6, height = 5, dpi = 300, bg = 'white')

as expected, temperature was a significant predictor in the full model and a full model of all was not significantly different than a model which only included temperature. pH had no influence.

another plot, ph

based on Tredennick which bases it off of https://biologyforfun.wordpress.com/2015/06/17/confidence-intervals-for-prediction-in-glmms/, should keep in mind that confidence or prediction intervals might underestimate uncertainty around our estimates.

juv.arctica_vec <- seq(min(juv.arctica$temp.average),max(juv.arctica$temp.average))

upper_ph <- quantile(juv.arctica$pH.normalized, 0.9) 
lower_ph <- quantile(juv.arctica$pH.normalized, 0.1)

newdat <- tibble(temp.average = juv.arctica_vec,
               Tank = juv.arctica$tank[1],
               pH.normalized = lower_ph,
               ph_level = "low") %>%
  bind_rows(
    tibble(temp.average = juv.arctica_vec,
           Tank = juv.arctica$tank[1],
           pH.normalized = upper_ph,
           ph_level = "high")
  ) %>%
  as.data.frame()

m <- full.m
allout <- tibble()  # empty object
for(do_level in unique(newdat$ph_level)) {
  tmpdat <- newdat %>%
    filter(ph_level == do_level)
  
  mm <- model.matrix(~pH.normalized*temp.average, tmpdat)
  tmpdat$y <- mm%*%fixef(m) 
  # predict(m, newdat, re.form=NA) would give the same results
  pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm))
  
  outdat <- data.frame(tmpdat,
                       lower_ci = tmpdat$y-1.96*sqrt(pvar1),
                       upper_ci = tmpdat$y+1.96*sqrt(pvar1))
  allout <- bind_rows(allout, outdat)
}

juv.arctica.int.ci.plot<-ggplot(allout, aes(x = temp.average, y = y,
                   color = ph_level, fill = ph_level)) +
  geom_hline(aes(yintercept = 0), linetype = 3) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), color = NA, alpha = 0.25) +
  geom_line() +
  scale_color_manual(values = rev(c("coral","dodgerblue")), 
                     labels = c("High pH", "Low pH"),
                     name = NULL) +
  scale_fill_manual(values = rev(c("coral","dodgerblue")),
                    labels = c("High pH", "Low pH"),
                    name = NULL) +
  #scale_y_continuous(breaks = c(0,1)) +
  #scale_x_continuous(breaks = c(6,9,12)) +
  labs(x = "Mean Temperature", y = "Proportional Change in Height") +
  theme_classic() +
  theme(legend.position = c(0.8, 0.85))
juv.arctica.int.ci.plot

ggsave("02_output/02_plots/height_juv.arctica_int.ci.plot.png", height=3.2, width=4, units = "in")

okay so this plot shows there is no interaction, but is otherwise pretty unhelpful, I am also a bit concerned that the axis scales don’t make sense. I have a feeling all our plots would look like this/the same bc it was only significant in one case, therefore I am not sure how helpful it is.

all combined model plots

all.model.plot<-ggarrange(mercenaria.model.plot,mya.model.plot, juv.arctica.model.plot,
              ncol = 1, nrow = 3)

ggsave("02_output/02_plots/all_height_mod_plot.png", all.model.plot, width = 5, height = 13, dpi = 300, bg = "white")