Set up

read in relevant files

I took the sheet Diana sent me and saved a csv with each tab of that sheet then read those in here.

remove dead and lost

We had 59 removed. This is lower than height because some dead were not measured.

## [1] 425  18
## [1] 366  18

look at data

## Looking for trends

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

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

Not totally sure if this is necessary for temperature too, but it makes sense to me.

## [1] 7.39
## [1] 6.14

Plot

## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:ggpubr':
## 
##     mutate
##    temp  pH  N      mean       sd       se        ci
## 1     6 7.4 25 124.31220 30.77729 6.155459 10.531264
## 2     6 7.6 27 111.38259 27.62450 5.316337  9.067639
## 3     6 7.8 24 102.63350 25.12788 5.129208  8.790803
## 4     6   8 27 107.49989 24.26398 4.669605  7.964562
## 5     9 7.4 46 127.29761 34.42511 5.075706  8.524280
## 6     9 7.6 45 114.15244 29.07736 4.334596  7.283119
## 7     9 7.8 42 107.31557 26.61527 4.106825  6.911286
## 8     9   8 51 105.50114 26.45628 3.704620  6.208592
## 9    12 7.4 18 123.82133 29.88144 7.043122 12.252263
## 10   12 7.6 22 122.27141 35.53425 7.575928 13.036225
## 11   12 7.8 21  99.05114 26.22239 5.722195  9.869175
## 12   12   8 18 110.24456 29.29434 6.904743 12.011537

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Individual Species Groups

mercenaria

check structure

Looks pretty normal to me!

## 'data.frame':    116 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ File path      : chr  NA "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\mercenaria\\4th round\\b41.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\mercenaria\\b38.jpeg" NA ...
##  $ file name      : chr  "b38.jpeg" "b42" "b39" "b47" ...
##  $ shell          : chr  "b38" "b42" "b39" "b47" ...
##  $ N              : num  456633 274746 316471 329233 343592 ...
##  $ mean           : num  146 140 142 126 138 ...
##  $ stdev          : num  24.1 16.8 22.7 20.9 21.1 ...
##  $ mode           : num  139 144 128 134 145 126 133 137 130 127 ...
##  $ 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 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "mercenaria" "mercenaria" "mercenaria" "mercenaria" ...
##  $ 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 ...
##  $ pH.normalized  : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 3.05 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mercenaria$mean
## W = 0.9899, p-value = 0.5504
##    vars   n   mean    sd median trimmed   mad   min    max range  skew kurtosis
## X1    1 116 121.68 11.03 122.41  121.71 12.29 96.69 146.08 49.39 -0.04    -0.67
##      se
## X1 1.02

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 with AIC includes temp and ph and interaction but none are significant??? resid vs predicted - quantile deviation significant, but it doesn’t seem like much trend to me. 7.8 group looks a bit odd, temp groups seem fine. CI on pH is huge…

The top 2 models with AIC are quite similar, when ranking by BIC interaction is not included and the model looks more similar to my results from the method from V1.

color.1 <- lmer(mean~temp.average*pH.normalized+(1|Tank), data = mercenaria, na.action = na.fail)

simulationOutput1<-simulateResiduals(color.1)

dd1<-dredge(color.1, rank=BIC) #BIC! more conservative
## Warning in dredge(color.1, rank = BIC): comparing models fitted by REML
## Fixed term is "(Intercept)"
dd1
## Global model call: lmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = mercenaria, na.action = na.fail)
## ---
## Model selection table 
##    (Int) pH.nrm tmp.avr pH.nrm:tmp.avr df   logLik   BIC delta weight
## 4 102.90 -19.63   2.779                 5 -388.293 800.4  0.00  0.510
## 8 108.60 -37.04   2.132          1.995  6 -386.010 800.5  0.19  0.465
## 2 128.70 -22.03                         4 -394.073 807.2  6.81  0.017
## 3  94.55          2.978                 4 -394.895 808.8  8.45  0.007
## 1 121.40                                3 -399.851 814.0 13.61  0.001
## Models ranked by BIC(x) 
## Random terms (all models): 
##   1 | Tank
top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1
## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ pH.normalized + temp.average + (1 | Tank)
##    Data: mercenaria
## REML criterion at convergence: 776.5863
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 5.190   
##  Residual             6.445   
## Number of obs: 116, groups:  Tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       102.895        -19.633          2.779
simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

#plot(top_model1)
tab_model(top_model1)
  mean
Predictors Estimates CI p
(Intercept) 102.90 88.70 – 117.09 <0.001
pH normalized -19.63 -32.18 – -7.08 0.002
temp average 2.78 1.35 – 4.21 <0.001
Random Effects
σ2 41.54
Ï„00 Tank 26.94
ICC 0.39
N Tank 16
Observations 116
Marginal R2 / Conditional R2 0.457 / 0.670 
testZeroInflation(top_model1)

## 
##  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(top_model1, mercenaria$pH)

plotResiduals(top_model1, mercenaria$temp)

plot_model(top_model1, 
                   axis.labels=c("Temp", "pH"),
                   show.values=TRUE, show.p=TRUE,
                   title="Effect of Temp & pH on Mercenaria Shell Coloration")+
  theme_classic()

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Temperature"),
                  dv.labels= "Effect of Temp & pH on Mercenaria Shell Coloration")
  Effect of Temp & pH on Mercenaria Shell Coloration
Predictors Estimates CI p
(Intercept) 102.90 88.70 – 117.09 <0.001
pH (normalized) -19.63 -32.18 – -7.08 0.002
Temperature 2.78 1.35 – 4.21 <0.001
Random Effects
σ2 41.54
Ï„00 Tank 26.94
ICC 0.39
N Tank 16
Observations 116
Marginal R2 / Conditional R2 0.457 / 0.670
formal plot

i think we want three lines here because temp and pH slopes - this seems pretty good but I would eventually like to be able to have just 4 lines for our 4 groups and to label them based on their non normalized numbers but this is tricky to do because pH in the model is CONTINUOUS so hard to plot effect size at discrete numbers?

values_of_pH <- data.frame(pH.normalized = c(0, 0.2, 0.4, 0.6))

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=top_model1, at=values_of_pH)
## NOTE: temp.averagepH.normalized does not appear in the model
effect_top_model.both<-as.data.frame(effect_top_model.both)
effect_top_model.both
##    temp.average pH.normalized      fit       se     lower    upper
## 1           6.1           0.0 119.8501 3.420371 113.07371 126.6265
## 2           7.6           0.0 124.0193 2.806780 118.45859 129.5801
## 3           9.1           0.0 128.1886 2.529103 123.17796 133.1992
## 4          11.0           0.0 133.4696 2.806733 127.90897 139.0303
## 5          12.0           0.0 136.2491 3.188253 129.93261 142.5656
## 6           6.1           0.2 115.9236 2.744892 110.48545 121.3617
## 7           7.6           0.2 120.0928 1.991343 116.14762 124.0380
## 8           9.1           0.2 124.2621 1.654004 120.98519 127.5390
## 9          11.0           0.2 129.5431 2.129844 125.32351 133.7627
## 10         12.0           0.2 132.3226 2.644126 127.08412 137.5611
## 11          6.1           0.3 113.9603 2.579877 108.84913 119.0715
## 12          7.6           0.3 118.1296 1.792414 114.57848 121.6807
## 13          9.1           0.3 122.2988 1.452239 119.42167 125.1760
## 14         11.0           0.3 127.5799 2.017192 123.58344 131.5763
## 15         12.0           0.3 130.3594 2.570643 125.26645 135.4523
## 16          6.1           0.5 110.0338 2.702888 104.67892 115.3887
## 17          7.6           0.5 114.2031 2.028428 110.18439 118.2218
## 18          9.1           0.5 118.3723 1.806296 114.79372 121.9509
## 19         11.0           0.5 123.6534 2.354122 118.98942 128.3173
## 20         12.0           0.5 126.4329 2.872105 120.74270 132.1230
## 21          6.1           0.7 106.1073 3.341407  99.48739 112.7273
## 22          7.6           0.7 110.2766 2.868075 104.59440 115.9587
## 23          9.1           0.7 114.4458 2.761517 108.97475 119.9169
## 24         11.0           0.7 119.7269 3.197542 113.39195 126.0618
## 25         12.0           0.7 122.5064 3.619300 115.33587 129.6768
custom_labels <- c("7.4", "7.6", "7.8", "8.0")

effect_temp <- ggplot() + 
  geom_line(data = effect_top_model.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
  geom_ribbon(data = effect_top_model.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 = mean, color = pH.normalized), size = 3) +
  labs(x = "Average Temperature (C)", y = "Mean Shell Coloration", color = "Average pH") +
  theme_classic() +
  scale_color_viridis(direction = -1, option = "cividis", breaks = c(0, 0.2, 0.4, 0.6), labels = custom_labels) +
  scale_fill_viridis(direction = -1, option = "cividis", breaks = c(0, 0.2, 0.4, 0.6), labels = custom_labels, guide="none")

effect_temp

#mercenaria.model.plot<-ggarrange(effect_temp,
             # labels = c("A"),
              #ncol = 1, nrow = 1)
ggsave("02_output/02_plots/color_mercenaria_linear.png", effect_temp, width = 7, height = 4, dpi = 300)
formal plot both
values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=top_model1, at=values_of_temp)
## NOTE: temp.averagepH.normalized does not appear in the model
effect_top_model.both<-as.data.frame(effect_top_model.both)
effect_top_model.both
##    temp.average pH.normalized      fit       se     lower    upper
## 1           6.1           0.0 119.8501 3.420371 113.07371 126.6265
## 2           7.6           0.0 124.0193 2.806780 118.45859 129.5801
## 3           9.1           0.0 128.1886 2.529103 123.17796 133.1992
## 4          11.0           0.0 133.4696 2.806733 127.90897 139.0303
## 5          12.0           0.0 136.2491 3.188253 129.93261 142.5656
## 6           6.1           0.2 115.9236 2.744892 110.48545 121.3617
## 7           7.6           0.2 120.0928 1.991343 116.14762 124.0380
## 8           9.1           0.2 124.2621 1.654004 120.98519 127.5390
## 9          11.0           0.2 129.5431 2.129844 125.32351 133.7627
## 10         12.0           0.2 132.3226 2.644126 127.08412 137.5611
## 11          6.1           0.3 113.9603 2.579877 108.84913 119.0715
## 12          7.6           0.3 118.1296 1.792414 114.57848 121.6807
## 13          9.1           0.3 122.2988 1.452239 119.42167 125.1760
## 14         11.0           0.3 127.5799 2.017192 123.58344 131.5763
## 15         12.0           0.3 130.3594 2.570643 125.26645 135.4523
## 16          6.1           0.5 110.0338 2.702888 104.67892 115.3887
## 17          7.6           0.5 114.2031 2.028428 110.18439 118.2218
## 18          9.1           0.5 118.3723 1.806296 114.79372 121.9509
## 19         11.0           0.5 123.6534 2.354122 118.98942 128.3173
## 20         12.0           0.5 126.4329 2.872105 120.74270 132.1230
## 21          6.1           0.7 106.1073 3.341407  99.48739 112.7273
## 22          7.6           0.7 110.2766 2.868075 104.59440 115.9587
## 23          9.1           0.7 114.4458 2.761517 108.97475 119.9169
## 24         11.0           0.7 119.7269 3.197542 113.39195 126.0618
## 25         12.0           0.7 122.5064 3.619300 115.33587 129.6768
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_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=temp.average), size=3) +
  labs(x="Average pH", y="Mean Shell Coloration", fill = "Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = FALSE)+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)
effect_temp
## Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
## ggplot2 3.3.4.
## ℹ Please use "none" instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=pH.normalized), size=3) +
  labs(x="Average Temperature (C)", y="Mean Shell Coloration", color = "Average pH")+
  theme_classic()+
  scale_color_viridis( option = "cividis", breaks = original_values, labels = custom_labels_x, direction = -1)+
  scale_fill_viridis( option= "cividis", guide = FALSE, direction = -1)
effect_temp_swap

mercenaria.model.plot<-ggarrange(effect_temp,effect_temp_swap,
             labels = c("A", "B"),
              ncol = 2, nrow = 1)
ggsave("02_output/02_plots/color_mercenaria_linear_both_swap.png", mercenaria.model.plot, width = 14, height = 4, dpi = 300)

Plot

This plot illustrates the trend in both pH and temperature for Mercenaria shell coloration.

##    temp  pH  N     mean        sd       se       ci
## 1     6 7.4  9 120.3091  5.836794 1.945598 3.617933
## 2     6 7.6  9 117.7914  5.050849 1.683616 3.130766
## 3     6 7.8  8 104.6461  4.402510 1.556522 2.948954
## 4     6   8  9 108.4222  7.452191 2.484064 4.619235
## 5     9 7.4 14 132.0149  8.297423 2.217580 3.927186
## 6     9 7.6 14 122.6619  8.737814 2.335279 4.135624
## 7     9 7.8 10 115.9668  5.505253 1.740914 3.191292
## 8     9   8 16 120.6646 10.680070 2.670018 4.680675
## 9    12 7.4  5 129.4374  3.401533 1.521212 3.242990
## 10   12 7.6  9 135.2718  6.861502 2.287167 4.253098
## 11   12 7.8  6 124.2305  7.934959 3.239433 6.527615
## 12   12   8  7 127.6544  5.246114 1.982845 3.853025

*** ## mya ### check structure Appears pretty normal.

## 'data.frame':    150 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r41.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r40.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r37.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r42.jpeg" ...
##  $ file name      : chr  "r41.jpeg" "r40.jpeg" "r37.jpeg" "r42.jpeg" ...
##  $ shell          : chr  "r41" "r40" "r37" "r42" ...
##  $ N              : num  661697 587852 417483 603232 286395 ...
##  $ mean           : num  130 140 158 139 118 ...
##  $ stdev          : num  45.9 46.1 39.6 46.9 51.9 ...
##  $ mode           : num  127 171 174 96 180 172 117 181 170 112 ...
##  $ 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 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "mya" "mya" "mya" "mya" ...
##  $ 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 ...
##  $ pH.normalized  : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 3.05 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mya$mean
## W = 0.98704, p-value = 0.1759
##    vars   n   mean    sd median trimmed  mad   min    max range skew kurtosis
## X1    1 150 132.61 15.32 129.99  132.29 15.1 96.47 172.74 76.27 0.21    -0.38
##      se
## X1 1.25

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 with AIC or BIC includes temp and ph and interaction but none are significant??? similar to what I saw above. this doesnt really make sense to me based on visual inspection of the data I therefore didn’t include interaction term in mode (temp.average+pH.normalized instead of temp.average*pH.normalized)

The top 2 models with AIC or BIC are quite similar, when ranking by BIC only pH is included - AIC pH and Temp but temp is not significant. Residuals look better for this model but I am not exactly sure how to interpret the temp term.

color.1 <- lmer(mean~temp.average+pH.normalized+(1|Tank), data = mya, na.action = na.fail)

simulationOutput1<-simulateResiduals(color.1)

dd1<-dredge(color.1, rank=AIC) #BIC = more conservative
## Warning in dredge(color.1, rank = AIC): comparing models fitted by REML
## Fixed term is "(Intercept)"
dd1
## Global model call: lmer(formula = mean ~ temp.average + pH.normalized + (1 | Tank), 
##     data = mya, na.action = na.fail)
## ---
## Model selection table 
##   (Int) pH.nrm tmp.avr df   logLik    AIC delta weight
## 4 162.3 -49.94 -1.4650  5 -545.348 1100.7  0.00  0.769
## 2 148.6 -48.71          4 -547.552 1103.1  2.41  0.231
## 3 140.7        -0.8963  4 -558.944 1125.9 25.19  0.000
## 1 132.6                 3 -560.532 1127.1 26.37  0.000
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank
top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1
## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ pH.normalized + temp.average + (1 | Tank)
##    Data: mya
## REML criterion at convergence: 1090.695
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 5.893   
##  Residual             8.780   
## Number of obs: 150, groups:  Tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       162.285        -49.937         -1.465
simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

#plot(top_model1)
tab_model(top_model1)
  mean
Predictors Estimates CI p
(Intercept) 162.28 145.89 – 178.68 <0.001
pH normalized -49.94 -64.30 – -35.57 <0.001
temp average -1.47 -3.11 – 0.18 0.081
Random Effects
σ2 77.08
Ï„00 Tank 34.73
ICC 0.31
N Tank 16
Observations 150
Marginal R2 / Conditional R2 0.540 / 0.683 
testZeroInflation(top_model1)

## 
##  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(top_model1, mya$pH)

plotResiduals(top_model1, mya$temp)

plot_model(top_model1, 
                   axis.labels=c("Temp", "pH"),
                   show.values=TRUE, show.p=TRUE,
                   title="Effect of Temp & pH on Mya Coloration")+
  theme_classic()

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp & pH on Mya Coloration")
  Effect of Temp & pH on Mya Coloration
Predictors Estimates CI p
(Intercept) 162.28 145.89 – 178.68 <0.001
pH (normalized) -49.94 -64.30 – -35.57 <0.001
Temp -1.47 -3.11 – 0.18 0.081
Random Effects
σ2 77.08
Ï„00 Tank 34.73
ICC 0.31
N Tank 16
Observations 150
Marginal R2 / Conditional R2 0.540 / 0.683
formal plot

i think we want one line as only pH significant

#values_of_pH <- data.frame(pH.normalized = c(0, 0.2, 0.4, 0.6))

effect_top_model<-effects::effect(term= c("pH.normalized"), mod=top_model1)
effect_top_model<-as.data.frame(effect_top_model)
effect_top_model
##   pH.normalized      fit       se    lower    upper
## 1           0.0 149.2078 2.909247 143.4585 154.9572
## 2           0.2 139.2203 1.896676 135.4721 142.9686
## 3           0.3 134.2266 1.657836 130.9503 137.5029
## 4           0.5 124.2391 2.054001 120.1799 128.2983
## 5           0.7 114.2516 3.148646 108.0291 120.4741
original_values <- c(0, 0.2, 0.4, 0.6)
custom_labels_x <- c("7.4", "7.6", "7.8", "8.0")

effect_pH <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model, aes(x = pH.normalized, y = fit), alpha = 0.7) +
  geom_ribbon(data= effect_top_model, aes(x=pH.normalized, ymin=lower, ymax=upper), alpha= 0.3) +
  geom_point(data=mya, mapping = aes(x=pH.normalized, y=mean, color=temp.average), size=3) +
  labs(x="Average pH", y="Mean Shell Coloration", color="Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option="plasma")+
  scale_fill_viridis(option ="plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)
effect_pH

#mercenaria.model.plot<-ggarrange(effect_temp,
             # labels = c("A"),
              #ncol = 1, nrow = 1)
ggsave("02_output/02_plots/color_mya_linear_ph.png", effect_pH, width = 7, height = 4, dpi = 300)
formal plot both
values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=top_model1, at=values_of_temp)
## NOTE: temp.averagepH.normalized does not appear in the model
effect_top_model.both<-as.data.frame(effect_top_model.both)
effect_top_model.both
##    temp.average pH.normalized      fit       se    lower    upper
## 1           6.1           0.0 153.3458 3.944119 145.5513 161.1403
## 2           7.6           0.0 151.1477 3.225761 144.7728 157.5225
## 3           9.1           0.0 148.9495 2.896430 143.2255 154.6735
## 4          11.0           0.0 146.1652 3.214795 139.8120 152.5184
## 5          12.0           0.0 144.6998 3.657704 137.4713 151.9282
## 6           6.1           0.2 143.3583 3.166702 137.1002 149.6165
## 7           7.6           0.2 141.1602 2.286697 136.6411 145.6792
## 8           9.1           0.2 138.9620 1.887930 135.2310 142.6930
## 9          11.0           0.2 136.1777 2.440316 131.3551 141.0003
## 10         12.0           0.2 134.7123 3.038648 128.7072 140.7173
## 11          6.1           0.3 138.3646 2.974588 132.4861 144.2431
## 12          7.6           0.3 136.1664 2.055460 132.1044 140.2285
## 13          9.1           0.3 133.9683 1.654077 130.6994 137.2371
## 14         11.0           0.3 131.1840 2.312891 126.6131 135.7548
## 15         12.0           0.3 129.7185 2.957174 123.8745 135.5626
## 16          6.1           0.5 128.3771 3.107561 122.2358 134.5184
## 17          7.6           0.5 126.1789 2.320603 121.5929 130.7650
## 18          9.1           0.5 123.9808 2.061014 119.9077 128.0538
## 19         11.0           0.5 121.1965 2.703059 115.8546 126.5383
## 20         12.0           0.5 119.7310 3.307007 113.1956 126.2664
## 21          6.1           0.7 118.3896 3.832787 110.8151 125.9641
## 22          7.6           0.7 116.1914 3.281773 109.7059 122.6770
## 23          9.1           0.7 113.9933 3.159769 107.7488 120.2377
## 24         11.0           0.7 111.2090 3.672622 103.9510 118.4669
## 25         12.0           0.7 109.7435 4.165605 101.5113 117.9757
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_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=temp.average), size=3) +
  labs(x="Average pH", y="Mean Shell Coloration", fill = "Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = FALSE)+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)
effect_temp

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.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=pH.normalized), size=3) +
  labs(x="Average Temperature (C)", y="Mean Shell Coloration", color = "Average pH")+
  theme_classic()+
  scale_color_viridis( option = "cividis", breaks = original_values, labels = custom_labels_x, direction = -1)+
  scale_fill_viridis( option= "cividis", guide = FALSE, direction = -1)
effect_temp_swap

mya.model.plot<-ggarrange(effect_temp,effect_temp_swap,
             labels = c("A", "B"),
              ncol = 2, nrow = 1)
ggsave("02_output/02_plots/color_mya_linear_both_swap.png", mya.model.plot, width = 14, height = 4, dpi = 300)

Plot

This plot depicts the strong relationship between pH and shell color in Mya, seemingly regardless of temperature.

##    temp  pH  N     mean        sd       se       ci
## 1     6 7.4 10 150.5219  5.412784 1.711673 3.137689
## 2     6 7.6 10 134.9325  5.662082 1.790508 3.282203
## 3     6 7.8  9 125.8381  6.966905 2.322302 4.318432
## 4     6   8 11 125.0058  5.317460 1.603274 2.905873
## 5     9 7.4 20 150.8367 14.428803 3.226378 5.578837
## 6     9 7.6 19 135.3885  9.583922 2.198703 3.812690
## 7     9 7.8 20 125.7649  7.166237 1.602419 2.770796
## 8     9   8 19 121.1604  8.098192 1.857853 3.221635
## 9    12 7.4  7 148.4751  6.220035 2.350952 4.568324
## 10   12 7.6  9 140.8783 11.784180 3.928060 7.304416
## 11   12 7.8  9 106.1281  8.849082 2.949694 5.485098
## 12   12   8  7 122.8366  4.856334 1.835522 3.566750

juv arctica

check structure

Appears pretty normal. Lower than others.

## 'data.frame':    66 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\juv Arctica\\g29.jpeg" NA NA NA ...
##  $ file name      : chr  "g29.jpeg" NA NA NA ...
##  $ shell          : chr  "g29" "g30" "g31" "g32" ...
##  $ N              : num  381123 339242 764948 602176 970122 ...
##  $ mean           : num  52.5 41 58.5 52.7 68.8 ...
##  $ stdev          : num  9.96 17.64 19.74 16.74 24.64 ...
##  $ mode           : num  47 28 46 51 68 47 46 50 47 58 ...
##  $ 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 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "juv.arctica" "juv.arctica" "juv.arctica" "juv.arctica" ...
##  $ 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 ...
##  $ pH.normalized  : num  0.2 0.2 0.2 0.2 0.61 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 3.05 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  juv.arctica$mean
## W = 0.9902, p-value = 0.8847
##    vars  n  mean   sd median trimmed  mad   min   max range  skew kurtosis   se
## X1    1 66 59.97 6.25  59.67   59.92 5.99 40.98 75.74 34.75 -0.03     0.43 0.77

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 with AIC includes temp and ph and interaction but none are significant similar to what I saw above, but makes sense here as we do not expect a relationship.

The top model with BIC is just pH but still not significant. Residuals look better for AIC model with all terms so moving forward with it.

color.1 <- lmer(mean~temp.average*pH.normalized+(1|Tank), data = juv.arctica, na.action = na.fail)

simulationOutput1<-simulateResiduals(color.1)

dd1<-dredge(color.1, rank=AIC) #BIC = more conservative
## Warning in dredge(color.1, rank = AIC): comparing models fitted by REML
## Fixed term is "(Intercept)"
dd1
## Global model call: lmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = juv.arctica, na.action = na.fail)
## ---
## Model selection table 
##   (Int) pH.nrm tmp.avr pH.nrm:tmp.avr df   logLik   AIC delta weight
## 8 67.23 -2.269 -0.5953        -0.3574  6 -207.039 426.1  0.00  0.494
## 4 68.29 -5.322 -0.7177                 5 -208.614 427.2  1.15  0.278
## 2 61.68 -4.936                         4 -210.057 428.1  2.04  0.178
## 3 65.99        -0.6649                 4 -211.918 431.8  5.76  0.028
## 1 60.00                                3 -213.146 432.3  6.21  0.022
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank
top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## mean ~ pH.normalized + temp.average + (1 | Tank) + pH.normalized:temp.average
##    Data: juv.arctica
## REML criterion at convergence: 414.0776
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 1.583   
##  Residual             5.976   
## Number of obs: 66, groups:  Tank, 16
## Fixed Effects:
##                (Intercept)               pH.normalized  
##                    67.2318                     -2.2687  
##               temp.average  pH.normalized:temp.average  
##                    -0.5953                     -0.3574
simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

#plot(top_model1)
tab_model(top_model1)
  mean
Predictors Estimates CI p
(Intercept) 67.23 53.09 – 81.37 <0.001
pH normalized -2.27 -36.79 – 32.25 0.896
temp average -0.60 -2.14 – 0.95 0.443
pH normalized × temp
average 		-0.36 -4.23 – 3.52 0.854
Random Effects
σ2 35.71
Ï„00 Tank 2.51
ICC 0.07
N Tank 16
Observations 66
Marginal R2 / Conditional R2 0.077 / 0.138 
testZeroInflation(top_model1)

## 
##  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(top_model1, juv.arctica$pH)

plotResiduals(top_model1, juv.arctica$temp)

plot_model(top_model1, 
                   axis.labels=c("pHxTemp","Temp", "pH"),
                   show.values=TRUE, show.p=TRUE,
                   title="Effect of Temp & pH on Juvenile Arctica Coloration")+
  theme_classic()

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Average Temperature","pHxTemperature"),
                  dv.labels= "Effect of Temp & pH on Juvenile Arctica Coloration")
  Effect of Temp & pH on Juvenile Arctica Coloration
Predictors Estimates CI p
(Intercept) 67.23 53.09 – 81.37 <0.001
pH (normalized) -2.27 -36.79 – 32.25 0.896
Average Temperature -0.60 -2.14 – 0.95 0.443
pHxTemperature -0.36 -4.23 – 3.52 0.854
Random Effects
σ2 35.71
Ï„00 Tank 2.51
ICC 0.07
N Tank 16
Observations 66
Marginal R2 / Conditional R2 0.077 / 0.138
formal plot

none our significant, this is to display that best fit line includes slope of 0 in CI, second plot is wack but I wont be fixing it.

#values_of_pH <- data.frame(pH.normalized = c(0, 0.2, 0.4, 0.6))

effect_top_model<-effects::effect(term= c("pH.normalized"), mod=top_model1)
## NOTE: pH.normalized is not a high-order term in the model
effect_top_model<-as.data.frame(effect_top_model)
effect_top_model
##   pH.normalized      fit        se    lower    upper
## 1           0.0 61.86720 1.5220857 58.82459 64.90980
## 2           0.2 60.76937 0.9883871 58.79361 62.74512
## 3           0.3 60.22045 0.8528985 58.51553 61.92537
## 4           0.5 59.12262 1.0341075 57.05547 61.18977
## 5           0.7 58.02479 1.5913851 54.84366 61.20592
original_values <- c(0, 0.2, 0.4, 0.6)
custom_labels_x <- c("7.4", "7.6", "7.8", "8.0")

effect_pH <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model, aes(x = pH.normalized, y = fit), alpha = 0.7) +
  geom_ribbon(data= effect_top_model, aes(x=pH.normalized, ymin=lower, ymax=upper), alpha= 0.3) +
  geom_point(data=juv.arctica, mapping = aes(x=pH.normalized, y=mean, color=temp.average), size=3) +
  labs(x="Average pH", y="Mean Shell Coloration", color="Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option = "plasma")+
  scale_fill_viridis(option = "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)
effect_pH

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

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

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized", "pH.normalized:temp.average"), mod=top_model1, at=values_of_temp)
effect_top_model.both<-as.data.frame(effect_top_model.both)
effect_top_model.both
##    temp.average pH.normalized      fit        se    lower    upper
## 1           6.1           0.0 63.60019 2.6706754 58.26159 68.93880
## 2           7.6           0.0 62.70718 1.8413610 59.02635 66.38800
## 3           9.1           0.0 61.81416 1.5259484 58.76383 64.86449
## 4          11.0           0.0 60.68300 2.1977740 56.28971 65.07629
## 5          12.0           0.0 60.08766 2.8060427 54.47846 65.69686
## 6           6.1           0.2 62.71043 1.7872351 59.13780 66.28306
## 7           7.6           0.2 61.71019 1.2336774 59.24411 64.17628
## 8           9.1           0.2 60.70996 0.9878956 58.73518 62.68473
## 9          11.0           0.2 59.44299 1.3787933 56.68682 62.19916
## 10         12.0           0.2 58.77617 1.7645888 55.24881 62.30353
## 11          6.1           0.3 62.26555 1.5486605 59.16982 65.36128
## 12          7.6           0.3 61.21170 1.0653220 59.08216 63.34125
## 13          9.1           0.3 60.15786 0.8527110 58.45331 61.86240
## 14         11.0           0.3 58.82299 1.1995414 56.42514 61.22084
## 15         12.0           0.3 58.12042 1.5378352 55.04634 61.19451
## 16          6.1           0.5 61.37579 1.7757275 57.82616 64.92542
## 17          7.6           0.5 60.21472 1.2063853 57.80319 62.62625
## 18          9.1           0.5 59.05366 1.0413501 56.97203 61.13529
## 19         11.0           0.5 57.58298 1.5872473 54.41012 60.75584
## 20         12.0           0.5 56.80894 2.0321646 52.74670 60.87118
## 21          6.1           0.7 60.48603 2.6527039 55.18335 65.78870
## 22          7.6           0.7 59.21774 1.7986725 55.62225 62.81324
## 23          9.1           0.7 57.94946 1.6067422 54.73763 61.16129
## 24         11.0           0.7 56.34297 2.5045700 51.33641 61.34953
## 25         12.0           0.7 55.49745 3.2006802 49.09938 61.89552
effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=temp.average), size=3) +
  labs(x="Average pH", y="Mean Shell Coloration", color="Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option = "plasma")+
  scale_fill_viridis(option = "plasma", guide=FALSE)+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)
effect_temp

#mercenaria.model.plot<-ggarrange(effect_temp,
             # labels = c("A"),
              #ncol = 1, nrow = 1)
ggsave("02_output/02_plots/color_juv.arctica_linear_both.png", effect_temp, width = 7, height = 4, dpi = 300)
formal plot both
values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=top_model1, at=values_of_temp)
## Warning in term == terms: longer object length is not a multiple of shorter
## object length
## NOTE: temp.averagepH.normalized does not appear in the model
## 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_top_model.both<-as.data.frame(effect_top_model.both)
effect_top_model.both
##    temp.average pH.normalized      fit        se    lower    upper
## 1           6.1           0.0 63.60019 2.6706754 58.26159 68.93880
## 2           7.6           0.0 62.70718 1.8413610 59.02635 66.38800
## 3           9.1           0.0 61.81416 1.5259484 58.76383 64.86449
## 4          11.0           0.0 60.68300 2.1977740 56.28971 65.07629
## 5          12.0           0.0 60.08766 2.8060427 54.47846 65.69686
## 6           6.1           0.2 62.71043 1.7872351 59.13780 66.28306
## 7           7.6           0.2 61.71019 1.2336774 59.24411 64.17628
## 8           9.1           0.2 60.70996 0.9878956 58.73518 62.68473
## 9          11.0           0.2 59.44299 1.3787933 56.68682 62.19916
## 10         12.0           0.2 58.77617 1.7645888 55.24881 62.30353
## 11          6.1           0.3 62.26555 1.5486605 59.16982 65.36128
## 12          7.6           0.3 61.21170 1.0653220 59.08216 63.34125
## 13          9.1           0.3 60.15786 0.8527110 58.45331 61.86240
## 14         11.0           0.3 58.82299 1.1995414 56.42514 61.22084
## 15         12.0           0.3 58.12042 1.5378352 55.04634 61.19451
## 16          6.1           0.5 61.37579 1.7757275 57.82616 64.92542
## 17          7.6           0.5 60.21472 1.2063853 57.80319 62.62625
## 18          9.1           0.5 59.05366 1.0413501 56.97203 61.13529
## 19         11.0           0.5 57.58298 1.5872473 54.41012 60.75584
## 20         12.0           0.5 56.80894 2.0321646 52.74670 60.87118
## 21          6.1           0.7 60.48603 2.6527039 55.18335 65.78870
## 22          7.6           0.7 59.21774 1.7986725 55.62225 62.81324
## 23          9.1           0.7 57.94946 1.6067422 54.73763 61.16129
## 24         11.0           0.7 56.34297 2.5045700 51.33641 61.34953
## 25         12.0           0.7 55.49745 3.2006802 49.09938 61.89552
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_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=temp.average), size=3) +
  labs(x="Average pH", y="Mean Shell Coloration", fill = "Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = FALSE)+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)
effect_temp

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.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.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=mean, color=pH.normalized), size=3) +
  labs(x="Average Temperature (C)", y="Mean Shell Coloration", color = "Average pH")+
  theme_classic()+
  scale_color_viridis( option = "cividis", breaks = original_values, labels = custom_labels_x, direction = -1)+
  scale_fill_viridis( option= "cividis", guide = FALSE, direction = -1)
effect_temp_swap

juv.arctica.model.plot<-ggarrange(effect_temp,effect_temp_swap,
             labels = c("A", "B"),
              ncol = 2, nrow = 1)
ggsave("02_output/02_plots/color_juv.arctica_linear_both_swap.png", juv.arctica.model.plot, width = 14, height = 4, dpi = 300)

Plot

These plots show how there really doesn’t seem to be much relationship between the color of juvenile Arctica shells and the pH or temperature treatments.

##    temp  pH  N     mean       sd       se        ci
## 1     6 7.4  4 64.92725 8.478716 4.239358  9.976750
## 2     6 7.6  4 58.22325 4.905617 2.452809  5.772350
## 3     6 7.8  4 58.21125 8.986877 4.493438 10.574694
## 4     6   8  4 61.39900 4.614352 2.307176  5.429624
## 5     9 7.4  8 61.37288 3.151784 1.114324  2.111174
## 6     9 7.6  8 62.50863 6.428547 2.272835  4.306064
## 7     9 7.8  8 59.77137 4.823539 1.705379  3.230974
## 8     9   8 11 60.06555 6.548514 1.974451  3.578616
## 9    12 7.4  3 66.82267 3.444560 1.988717  5.807026
## 10   12 7.6  4 51.15500 7.338287 3.669144  8.634829
## 11   12 7.8  4 55.24025 3.043270 1.521635  3.580960
## 12   12   8  4 57.74125 2.830496 1.415248  3.330593