this script will be the same as 4 except with changes/additions based on discussion with Diana (reduce to one comparison). Would also like to combine all final species plots. Also removed scallops.

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] 389  15

## [1] 331  15

## [1] 3

look at data

## `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 3 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 3 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 3 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

summary_by_species

## 
##  Descriptive statistics by group 
## 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  1.00  0.00   1.00    1.00  0.00  1.00   1.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  2.00  0.00   2.00    2.00  0.00 2.00   2.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  3.00  0.00   3.00    3.00  0.00  3.00   3.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

#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

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/ 3 levels "juv.arctica",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ 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?

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

simulationOutput1<-simulateResiduals(height.1)

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

shapiro.test(resid(height.1))

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

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?

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

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

all.v.temp<-anova(all.m, temp.m)

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

all.v.temp

## Data: mercenaria
## Models:
## temp.m: prop.change.height ~ temp.average + (1 | tank)
## all.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##        npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## temp.m    4 -251.43 -240.45 129.72  -259.43                     
## all.m     6 -248.26 -231.79 130.13  -260.26 0.8306  2     0.6601

summary(all.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: -239.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6900 -0.2228  0.0242  0.5214  2.7104 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.006773 0.08230 
##  Residual             0.004525 0.06727 
## Number of obs: 115, groups:  tank, 16
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)                -0.51824    0.17530 10.84260  -2.956  0.01325 * 
## temp.average                0.07669    0.01897 10.95874   4.043  0.00195 **
## pH.normalized              -0.22052    0.43074 10.85520  -0.512  0.61894   
## temp.average:pH.normalized  0.03020    0.04804 10.97435   0.629  0.54242   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr pH.nrm
## temp.averag -0.976              
## pH.normalzd -0.818  0.811       
## tmp.vrg:pH.  0.787 -0.819 -0.975

summary(temp.m)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: mercenaria
## 
## REML criterion at convergence: -246
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7086 -0.2398  0.0233  0.5195  2.7099 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.006016 0.07756 
##  Residual             0.004526 0.06728 
## Number of obs: 115, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  -0.58596    0.09504 12.71658  -6.165 3.75e-05 ***
## temp.average  0.08597    0.01025 12.80303   8.384 1.48e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.976

#both these give essentially the same answer for coefficent of temp

tab_model(all.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	-0.52	-0.87 – -0.17	0.004
temp average	0.08	0.04 – 0.11	<0.001
pH normalized	-0.22	-1.07 – 0.63	0.610
temp average × pH normalized	0.03	-0.07 – 0.13	0.531
Random Effects
σ²	0.00
τ₀₀ _tank	0.01
ICC	0.60
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.739 / 0.896

tab_model(temp.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	-0.59	-0.77 – -0.40	<0.001
temp average	0.09	0.07 – 0.11	<0.001
Random Effects
σ²	0.00
τ₀₀ _tank	0.01
ICC	0.57
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.750 / 0.893

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.

effect_top_model.temp<-effects::effect(term="temp.average", mod=all.m)

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

effect_top_model.temp<-as.data.frame(effect_top_model.temp)

effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  geom_point(data=mercenaria, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) +
  labs(x="Temperature (C)", y="Proportional Change in Height", color = "pH")+
  theme_classic()+
  scale_color_viridis(direction = -1, option = "cividis")
effect_temp

#mercenaria.model.plot<-ggarrange(effect_temp,
             # labels = c("A"),
              #ncol = 1, nrow = 1)
ggsave("02_output/02_plots/height_mercenaria_linear.png", effect_temp, width = 5, height = 4, dpi = 300)

effect_temp_swap <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
 # geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  geom_point(data=mercenaria, mapping = aes(x=pH.average, y=prop.change.height, color=temp.average), size=3) +
  labs(x="pH", y="Proportional Change in Height", color = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma")
effect_temp_swap

mercenaria.model.plot<-ggarrange(effect_temp_swap,effect_temp,
             labels = c("A"),
              ncol = 2, nrow = 1)
mercenaria.model.plot<-annotate_figure(mercenaria.model.plot, top = text_grob("Mercenaria mercenaria", 
               color = "black", face = 3, size = 22))

ggsave("02_output/02_plots/height_mercenaria_linear_both_swap.png", mercenaria.model.plot, width = 14, height = 4, dpi = 300)

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 <- all.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_none_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/ 3 levels "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.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 none, which seems to make sense based on visual inspection of the data. residuals etc look pretty good but I will try a couple other families. I assume this will be best.

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

simulationOutput1<-simulateResiduals(height.1)

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$ph)

plotResiduals(height.1, mya$temp)

shapiro.test(resid(height.1))

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

#plot_model(top_model1, 
                  # axis.labels=c("Temp"),
                   #show.values=TRUE, show.p=TRUE,
                  # title="Effect of Temp on mya Growth")+
 # theme_classic()

make models and likelihood ratio test

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

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

all.v.temp<-anova(all.m, temp.m)

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

all.v.temp

## Data: mya
## Models:
## temp.m: prop.change.height ~ temp.average + (1 | tank)
## all.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##        npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## temp.m    4 -217.63 -205.59 112.82  -225.63                     
## all.m     6 -213.91 -195.85 112.96  -225.91 0.2784  2       0.87

summary(all.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: -201.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.27851 -0.57911  0.02336  0.68283  2.21693 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.001316 0.03627 
##  Residual             0.012441 0.11154 
## Number of obs: 150, groups:  tank, 16
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)  
## (Intercept)                 0.069097   0.105664 12.317557   0.654   0.5252  
## temp.average                0.021707   0.011536 12.862748   1.882   0.0827 .
## pH.normalized              -0.016417   0.259894 12.331665  -0.063   0.9506  
## temp.average:pH.normalized -0.001191   0.029289 12.958119  -0.041   0.9682  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr pH.nrm
## temp.averag -0.976              
## pH.normalzd -0.819  0.813       
## tmp.vrg:pH.  0.787 -0.821 -0.976

summary(temp.m)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: mya
## 
## REML criterion at convergence: -210.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2953 -0.6103  0.0161  0.6556  2.2966 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  tank     (Intercept) 0.0009807 0.03132 
##  Residual             0.0124407 0.11154 
## Number of obs: 150, groups:  tank, 16
## 
## Fixed effects:
##               Estimate Std. Error        df t value Pr(>|t|)   
## (Intercept)   0.060265   0.056452 14.160598   1.068  0.30359   
## temp.average  0.021708   0.006141 14.695341   3.535  0.00309 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.977

#both these give essentially the same answer for coefficent of temp

tab_model(all.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.07	-0.14 – 0.28	0.514
temp average	0.02	-0.00 – 0.04	0.062
pH normalized	-0.02	-0.53 – 0.50	0.950
temp average × pH normalized	-0.00	-0.06 – 0.06	0.968
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.10
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.116 / 0.201

tab_model(temp.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.06	-0.05 – 0.17	0.287
temp average	0.02	0.01 – 0.03	0.001
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.07
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.118 / 0.182

This looks very similar to mercenaria, interaction not significant. temperature significant, tempp only not sig different from all.

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

effect_top_model.temp<-effects::effect(term="temp.average", mod=all.m)

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

effect_top_model.temp<-as.data.frame(effect_top_model.temp)

effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  geom_point(data=mya, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) +
  labs(x="Temperature (C)", y="Proportional Change in Height", color = "pH")+
  theme_classic()+
  scale_color_viridis(direction = -1, option = "cividis")
effect_temp

#mya.model.plot<-ggarrange(effect_temp,
             # labels = c("A"),
              #ncol = 1, nrow = 1)
ggsave("02_output/02_plots/height_mya_none_linear.png", effect_temp, width = 5, height = 4, dpi = 300)

effect_none_pH <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  #geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  geom_point(data=mya, mapping = aes(x=pH.average, y=prop.change.height, color=temp.average), size=3) +
  labs(x="pH", y="Proportional Change in Height", color = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option = "plasma")
effect_none_pH

mya.model.plot<-ggarrange(effect_none_pH,effect_temp,
             labels = c("B"),
              ncol = 2, nrow = 1)
mya.model.plot<-annotate_figure(mya.model.plot, top = text_grob("Mya arenaria", 
               color = "black", face = 3, size = 22))

ggsave("02_output/02_plots/height_mya_linear_both_swap.png", mya.model.plot, width = 14, height = 4, dpi = 300)

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/ 3 levels "juv.arctica",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ 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, which seems to make sense based on visual inspection of the data & the biology of arctica. residuals etc look pretty good but I will try a couple other families. I assume this will be best.

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

simulationOutput1<-simulateResiduals(height.1)

plot(height.1)

shapiro.test(resid(height.1))

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

#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

make models and likelihood ratio test

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

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

all.v.temp<-anova(all.m, temp.m)

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

all.v.temp

## Data: juv.arctica
## Models:
## temp.m: prop.change.height ~ temp.average + (1 | tank)
## all.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##        npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## temp.m    4 -46.763 -38.190 27.381  -54.763                     
## all.m     6 -44.256 -31.397 28.128  -56.256 1.4935  2     0.4739

summary(all.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: -35.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5903 -0.6584  0.0693  0.5605  3.7929 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  tank     (Intercept) 0.0009975 0.03158 
##  Residual             0.0247960 0.15747 
## Number of obs: 63, groups:  tank, 16
## 
## Fixed effects:
##                            Estimate Std. Error       df t value Pr(>|t|)  
## (Intercept)                 0.27199    0.17774 14.33755   1.530   0.1477  
## temp.average                0.04530    0.01945 15.03784   2.330   0.0342 *
## pH.normalized              -0.47274    0.43363 14.00792  -1.090   0.2940  
## temp.average:pH.normalized  0.05378    0.04879 14.59288   1.102   0.2882  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr pH.nrm
## temp.averag -0.974              
## pH.normalzd -0.826  0.819       
## tmp.vrg:pH.  0.795 -0.830 -0.975

summary(temp.m)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: -41.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8062 -0.6487  0.0213  0.6443  3.8253 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  tank     (Intercept) 0.0007874 0.02806 
##  Residual             0.0246356 0.15696 
## Number of obs: 63, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   0.11391    0.09806 16.53624   1.162    0.262    
## temp.average  0.06312    0.01063 16.99098   5.937 1.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.977

#both these give essentially the same answer for coefficent of temp

tab_model(all.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.27	-0.08 – 0.63	0.131
temp average	0.05	0.01 – 0.08	0.023
pH normalized	-0.47	-1.34 – 0.40	0.280
temp average × pH normalized	0.05	-0.04 – 0.15	0.275
Random Effects
σ²	0.02
τ₀₀ _tank	0.00
ICC	0.04
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.388 / 0.412

tab_model(temp.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.11	-0.08 – 0.31	0.250
temp average	0.06	0.04 – 0.08	<0.001
Random Effects
σ²	0.02
τ₀₀ _tank	0.00
ICC	0.03
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.383 / 0.402

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

effect_top_model.temp<-effects::effect(term="temp.average", mod=all.m)

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

effect_top_model.temp<-as.data.frame(effect_top_model.temp)

effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  geom_point(data=juv.arctica, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) +
  labs(x="Temperature (C)", y="Proportional Change in Height", color="pH")+
  theme_classic()+
  scale_color_viridis(direction = -1, option = "cividis")
effect_temp

#juv.arctica.model.plot<-ggarrange(effect_temp,
             # labels = c("A"),
              #ncol = 1, nrow = 1)
ggsave("02_output/02_plots/height_juv.arctica_temp_linear.png", effect_temp, width = 5, height = 4, dpi = 300)

effect_temp_swap <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  #geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  geom_point(data=juv.arctica, mapping = aes(x=pH.average, y=prop.change.height, color=temp.average), size=3) +
  labs(x="pH", y="Proportional Change in Height", color="Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option = "plasma")
effect_temp_swap

juv.arctica.model.plot<-ggarrange(effect_temp_swap,effect_temp,
             labels = c("C"),
              ncol = 2, nrow = 1)
juv.arctica.model.plot<-annotate_figure(juv.arctica.model.plot, top = text_grob("Arctica islandica", 
               color = "black", face = 3, size = 22))
ggsave("02_output/02_plots/height_juv.arctica_linear_both_swap.png", juv.arctica.model.plot, width = 14, height = 4, dpi = 300)

Testing this way shows the same thing, just in a clearer way and without worrying about all that arbitary ‘best fit’ that didn’t make sense in our context. Juvenile arctica show the same pattern, temp is significant, pH and interaction are not.

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 = 10, height = 12, dpi = 300, bg = "white")

A-Maine-Tank-2022_Bivalve-Height_V6

McMahon

2023-09-27

Set up

read in relevant files

remove dead and lost

look at data

strange specimens

look for trends

summary table

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

Individual Species Groups

mercenaria

check structure

linear model

‘normal’ linear mixed effect model

check the normality assumptions of full model

make models and likelihood ratio test

formal plot

another plot, ph

mya

check structure

linear model

‘normal’ linear mixed effect model

make models and likelihood ratio test

formal plot

juvenile arctica

check structure

linear model

‘normal’ linear mixed effect model

make models and likelihood ratio test

formal plot

all combined model plots