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.

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] 389  18
## [1] 332  18

look at data

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 23 123.8135 31.52645 6.573718 11.288024
## 2     6 7.6 23 114.8844 28.20548 5.881249 10.098954
## 3     6 7.8 21 104.8837 25.91503 5.655123  9.753494
## 4     6   8 24 108.1858 23.46086 4.788928  8.207608
## 5     9 7.4 42 127.5220 35.29663 5.446388  9.165607
## 6     9 7.6 41 116.8224 28.96624 4.523767  7.617350
## 7     9 7.8 38 109.2931 26.96463 4.374241  7.379753
## 8     9   8 46 106.3783 27.62087 4.072475  6.839427
## 9    12 7.4 15 125.7987 32.07403 8.281480 14.586254
## 10   12 7.6 22 122.2714 35.53425 7.575928 13.036225
## 11   12 7.8 19 101.1314 26.68662 6.122329 10.616509
## 12   12   8 18 110.2446 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.

check the normality assumptions of full model

looks good

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

simulationOutput1<-simulateResiduals(color.1)
shapiro.test(resid(color.1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(color.1)
## W = 0.99444, p-value = 0.9282
plot(color.1)

#plot(top_model1)
testZeroInflation(color.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(top_model1, mercenaria$pH)
#plotResiduals(top_model1, mercenaria$temp)
make models and likelihood ratio test

when talking with Diana we decided that 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(mean~temp.average*pH.normalized+(1|Tank), data = mercenaria, na.action = na.fail)

null.m<-lmer(mean~(1|Tank), data = mercenaria, na.action = na.fail)

all.v.null<-anova(all.m, null.m) # different, but why
## refitting model(s) with ML (instead of REML)
all.v.null
## Data: mercenaria
## Models:
## null.m: mean ~ (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)    
## null.m    3 809.23 817.49 -401.61   803.23                        
## all.m     6 797.02 813.54 -392.51   785.02 18.21  3   0.000398 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(all.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.average * pH.normalized + (1 | Tank)
##    Data: mercenaria
## 
## REML criterion at convergence: 772
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.21899 -0.63490 -0.00922  0.61285  2.33041 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 28.77    5.363   
##  Residual             41.52    6.444   
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##                            Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                 108.612     11.940  10.841   9.097  2.1e-06 ***
## temp.average                  2.132      1.295  11.063   1.646    0.128    
## pH.normalized               -37.038     29.347  10.875  -1.262    0.233    
## temp.average:pH.normalized    1.995      3.281  11.099   0.608    0.555    
## ---
## 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.811       
## tmp.vrg:pH.  0.787 -0.820 -0.975
summary(null.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ (1 | Tank)
##    Data: mercenaria
## 
## REML criterion at convergence: 799.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12038 -0.60171 -0.01291  0.60430  2.16355 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 82.97    9.109   
##  Residual             41.38    6.432   
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  121.450      2.364  14.939   51.38   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(all.m)
  mean
Predictors Estimates CI p
(Intercept) 108.61 84.95 – 132.27 <0.001
temp average 2.13 -0.43 – 4.70 0.103
pH normalized -37.04 -95.20 – 21.12 0.210
temp average × pH
normalized 		1.99 -4.51 – 8.50 0.545
Random Effects
σ2 41.52
Ï„00 Tank 28.77
ICC 0.41
N Tank 16
Observations 116
Marginal R2 / Conditional R2 0.458 / 0.680 
tab_model(null.m)
  mean
Predictors Estimates CI p
(Intercept) 121.45 116.77 – 126.13 <0.001
Random Effects
σ2 41.38
Ï„00 Tank 82.97
ICC 0.67
N Tank 16
Observations 116
Marginal R2 / Conditional R2 0.000 / 0.667 
#######################

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

main.v.null<-anova(main_effect.m, null.m)
## refitting model(s) with ML (instead of REML)
main.v.null # all main effects is different than none
## Data: mercenaria
## Models:
## null.m: mean ~ (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## null.m           3 809.23 817.49 -401.61   803.23                         
## main_effect.m    5 795.51 809.28 -392.75   785.51 17.718  2  0.0001421 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
main.v.all<-anova(main_effect.m, all.m)
## refitting model(s) with ML (instead of REML)
main.v.all # all main effects is different than all
## Data: mercenaria
## Models:
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## main_effect.m    5 795.51 809.28 -392.75   785.51                     
## all.m            6 797.02 813.54 -392.51   785.02 0.4926  1     0.4828
main.v.all.vnull<-anova(all.m,main_effect.m, null.m)
## refitting model(s) with ML (instead of REML)
main.v.all.vnull 
## Data: mercenaria
## Models:
## null.m: mean ~ (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## null.m           3 809.23 817.49 -401.61   803.23                          
## main_effect.m    5 795.51 809.28 -392.75   785.51 17.7179  2  0.0001421 ***
## all.m            6 797.02 813.54 -392.51   785.02  0.4926  1  0.4827662    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(main_effect.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: mercenaria
## 
## REML criterion at convergence: 776.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22557 -0.59810 -0.00499  0.60233  2.34994 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 26.94    5.190   
##  Residual             41.54    6.445   
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    102.895      7.165  11.640  14.361 9.23e-09 ***
## temp.average     2.780      0.721  11.795   3.855  0.00236 ** 
## pH.normalized  -19.633      6.334  11.987  -3.100  0.00921 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.936       
## pH.normalzd -0.374  0.092
summary(main_effect.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: mercenaria
## 
## REML criterion at convergence: 776.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22557 -0.59810 -0.00499  0.60233  2.34994 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 26.94    5.190   
##  Residual             41.54    6.445   
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    102.895      7.165  11.640  14.361 9.23e-09 ***
## temp.average     2.780      0.721  11.795   3.855  0.00236 ** 
## pH.normalized  -19.633      6.334  11.987  -3.100  0.00921 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.936       
## pH.normalzd -0.374  0.092
tab_model(main_effect.m)
  mean
Predictors Estimates CI p
(Intercept) 102.90 88.70 – 117.09 <0.001
temp average 2.78 1.35 – 4.21 <0.001
pH normalized -19.63 -32.18 – -7.08 0.002
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 
tab_model(null.m,main_effect.m, all.m)
  mean mean mean
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 121.45 116.77 – 126.13 <0.001 102.90 88.70 – 117.09 <0.001 108.61 84.95 – 132.27 <0.001
temp average 2.78 1.35 – 4.21 <0.001 2.13 -0.43 – 4.70 0.103
pH normalized -19.63 -32.18 – -7.08 0.002 -37.04 -95.20 – 21.12 0.210
temp average × pH
normalized 		1.99 -4.51 – 8.50 0.545
Random Effects
σ2 41.38 41.54 41.52
Ï„00 82.97 Tank 26.94 Tank 28.77 Tank
ICC 0.67 0.39 0.41
N 16 Tank 16 Tank 16 Tank
Observations 116 116 116
Marginal R2 / Conditional R2 0.000 / 0.667 0.457 / 0.670 0.458 / 0.680

okay I dont know how to make sense of this… I see that these two are different but none are significant in the summary of the full model… I wonder what could cause this? VS. when compared to main effects.. ** should get dianas feedback ** seems like the estimates of main effects vs full aren’t that different? maybe this is a df problem from interaction term.

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?

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

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

#effect_top_model.both<-as.data.frame(effect_top_model.both)
#effect_top_model.both

#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 = "Temperature (C)", y = "Mean Shell Coloration", color = "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_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=all.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 121.4059 4.675654 112.14171 130.6701
## 2             9           0.0 127.8026 2.608967 122.63330 132.9720
## 3            12           0.0 134.1994 4.684550 124.91752 143.4812
## 4             6           0.2 116.3919 3.136092 110.17813 122.6057
## 5             9           0.2 123.9854 1.702243 120.61262 127.3582
## 6            12           0.2 131.5789 2.968334 125.69753 137.4603
## 7             6           0.4 111.3779 2.752368 105.92440 116.8313
## 8             9           0.4 120.1682 1.563516 117.07026 123.2661
## 9            12           0.4 128.9585 2.885508 123.24120 134.6757
## 10            6           0.6 106.3638 3.883700  98.66878 114.0589
## 11            9           0.6 116.3509 2.334001 111.72641 120.9755
## 12           12           0.6 126.3380 4.526639 117.36907 135.3070
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=mean, color=temp.average), size=3) +
  labs(x="pH", y="Mean Shell Coloration", fill = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = FALSE, breaks = c(6, 9, 12))+
  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_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=mean, color=pH.normalized), size=3) +
   labs(x="Temperature (C)", y="Mean Shell Coloration", color = "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"),
             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/color_mercenaria_linear_both_swap.png", mercenaria.model.plot, width = 14, height = 4, dpi = 300, bg = 'white')

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.

check the normality assumptions of full model

looks good

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

simulationOutput1<-simulateResiduals(color.1)

plot(color.1)

shapiro.test(resid(color.1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(color.1)
## W = 0.98572, p-value = 0.1249
#plot(top_model1)
testZeroInflation(color.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(top_model1, mercenaria$pH)
#plotResiduals(top_model1, mercenaria$temp)
make models and likelihood ratio test

when talking with Diana we decided that 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(mean~temp.average*pH.normalized+(1|Tank), data = mya, na.action = na.fail)

null.m<-lmer(mean~(1|Tank), data = mya, na.action = na.fail)

all.v.null<-anova(all.m, null.m) #not different, interaction not significant
## refitting model(s) with ML (instead of REML)
all.v.null
## Data: mya
## Models:
## null.m: mean ~ (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## null.m    3 1131.3 1140.3 -562.64   1125.3                         
## all.m     6 1110.3 1128.3 -549.14   1098.3 26.984  3  5.932e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(all.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.average * pH.normalized + (1 | Tank)
##    Data: mya
## 
## REML criterion at convergence: 1084.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.99104 -0.57071 -0.00478  0.71883  2.69786 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 32.14    5.669   
##  Residual             77.14    8.783   
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##                             Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                148.81995   13.01655  11.75742  11.433 1.01e-07 ***
## temp.average                 0.05589    1.41166  11.98722   0.040    0.969    
## pH.normalized               -8.84105   31.99990  11.78400  -0.276    0.787    
## temp.average:pH.normalized  -4.71195    3.57909  12.04557  -1.317    0.212    
## ---
## 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.786 -0.819 -0.975
summary(null.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ (1 | Tank)
##    Data: mya
## 
## REML criterion at convergence: 1121.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.93080 -0.56309  0.02491  0.67085  2.76181 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 168.34   12.974  
##  Residual              77.04    8.777  
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  132.608      3.323  14.990   39.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#both these give essentially the same answer for coefficent of temp

tab_model(all.m)
  mean
Predictors Estimates CI p
(Intercept) 148.82 123.09 – 174.55 <0.001
temp average 0.06 -2.73 – 2.85 0.968
pH normalized -8.84 -72.09 – 54.41 0.783
temp average × pH
normalized 		-4.71 -11.79 – 2.36 0.190
Random Effects
σ2 77.14
Ï„00 Tank 32.14
ICC 0.29
N Tank 16
Observations 150
Marginal R2 / Conditional R2 0.547 / 0.680 
tab_model(null.m)
  mean
Predictors Estimates CI p
(Intercept) 132.61 126.04 – 139.17 <0.001
Random Effects
σ2 77.04
Ï„00 Tank 168.34
ICC 0.69
N Tank 16
Observations 150
Marginal R2 / Conditional R2 0.000 / 0.686 
#######################

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

main.v.null<-anova(main_effect.m, null.m)
## refitting model(s) with ML (instead of REML)
main.v.null # all is different than none
## Data: mya
## Models:
## null.m: mean ~ (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## null.m           3 1131.3 1140.3 -562.64   1125.3                         
## main_effect.m    5 1110.5 1125.5 -550.23   1100.5 24.802  2  4.115e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
main.v.all<-anova(main_effect.m, all.m)
## refitting model(s) with ML (instead of REML)
main.v.all # all main effects is different than all
## Data: mya
## Models:
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## main_effect.m    5 1110.5 1125.5 -550.23   1100.5                     
## all.m            6 1110.3 1128.3 -549.14   1098.3 2.1824  1     0.1396
main.v.all.vnull<-anova(all.m,main_effect.m, null.m)
## refitting model(s) with ML (instead of REML)
main.v.all.vnull 
## Data: mya
## Models:
## null.m: mean ~ (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## null.m           3 1131.3 1140.3 -562.64   1125.3                          
## main_effect.m    5 1110.5 1125.5 -550.23   1100.5 24.8019  2  4.115e-06 ***
## all.m            6 1110.3 1128.3 -549.14   1098.3  2.1824  1     0.1396    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(main_effect.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: mya
## 
## REML criterion at convergence: 1090.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.98130 -0.54974  0.02177  0.74692  2.70970 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 34.73    5.893   
##  Residual             77.08    8.780   
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   162.2850     8.2941  12.8865  19.566 5.78e-11 ***
## temp.average   -1.4654     0.8346  13.0584  -1.756    0.103    
## pH.normalized -49.9375     7.2668  12.8403  -6.872 1.21e-05 ***
## ---
## 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.376  0.097
tab_model(main_effect.m)
  mean
Predictors Estimates CI p
(Intercept) 162.28 145.89 – 178.68 <0.001
temp average -1.47 -3.11 – 0.18 0.081
pH normalized -49.94 -64.30 – -35.57 <0.001
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 
tab_model(null.m,main_effect.m, all.m)
  mean mean mean
Predictors Estimates CI p Estimates CI p Estimates CI p
(Intercept) 132.61 126.04 – 139.17 <0.001 162.28 145.89 – 178.68 <0.001 148.82 123.09 – 174.55 <0.001
temp average -1.47 -3.11 – 0.18 0.081 0.06 -2.73 – 2.85 0.968
pH normalized -49.94 -64.30 – -35.57 <0.001 -8.84 -72.09 – 54.41 0.783
temp average × pH
normalized 		-4.71 -11.79 – 2.36 0.190
Random Effects
σ2 77.04 77.08 77.14
Ï„00 168.34 Tank 34.73 Tank 32.14 Tank
ICC 0.69 0.31 0.29
N 16 Tank 16 Tank 16 Tank
Observations 150 150 150
Marginal R2 / Conditional R2 0.000 / 0.686 0.540 / 0.683 0.547 / 0.680

We see the same thing here that we did for mercenaria… ** why **

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

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

formal plot both

In the future I think we just want one line, here I will leave it all in the plot to make sure we are consistent.

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

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=all.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 149.1553 5.089799 139.0961 159.2145
## 2             9           0.0 149.3229 2.819453 143.7507 154.8951
## 3            12           0.0 149.4906 5.085560 139.4398 159.5414
## 4             6           0.2 141.7327 3.412069 134.9893 148.4761
## 5             9           0.2 139.0732 1.833839 135.4489 142.6975
## 6            12           0.2 136.4137 3.220383 130.0491 142.7783
## 7             6           0.4 134.3102 2.989733 128.4014 140.2189
## 8             9           0.4 128.8235 1.681199 125.5009 132.1461
## 9            12           0.4 123.3368 3.140786 117.1295 129.5441
## 10            6           0.6 126.8876 4.218403 118.5506 135.2246
## 11            9           0.6 118.5738 2.517878 113.5976 123.5500
## 12           12           0.6 110.2599 4.933958 100.5087 120.0111
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=mean, color=temp.average), size=3) +
  labs(x="pH", y="Mean Shell Coloration", fill = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = FALSE, breaks = c(6, 9, 12))+
  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_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=mean, color=pH.normalized), size=3) +
   labs(x="Temperature (C)", y="Mean Shell Coloration", color = "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("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/color_juv.arctica_linear_both_swap.png", mya.model.plot, width = 14, height = 4, dpi = 300, bg = 'white')
plot using main effect model

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.

check the normality assumptions of full model

looks good

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

simulationOutput1<-simulateResiduals(color.1)

plot(color.1)

#plot(top_model1)
shapiro.test(resid(color.1))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(color.1)
## W = 0.98992, p-value = 0.8723
testZeroInflation(color.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(top_model1, mercenaria$pH)
#plotResiduals(top_model1, mercenaria$temp)
make models and likelihood ratio test

when talking with Diana we decided that 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(mean~temp.average*pH.normalized+(1|Tank), data = juv.arctica, na.action = na.fail)

null.m<-lmer(mean~(1|Tank), data = juv.arctica, na.action = na.fail)

all.v.null<-anova(all.m, null.m) # different, but why
## refitting model(s) with ML (instead of REML)
all.v.null
## Data: juv.arctica
## Models:
## null.m: mean ~ (1 | Tank)
## all.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## null.m    3 433.84 440.41 -213.92   427.84                     
## all.m     6 434.68 447.81 -211.34   422.68 5.1652  3     0.1601
summary(all.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ temp.average * pH.normalized + (1 | Tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: 414.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.69260 -0.71032 -0.01046  0.51751  2.20635 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept)  2.506   1.583   
##  Residual             35.707   5.976   
## Number of obs: 66, groups:  Tank, 16
## 
## Fixed effects:
##                            Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)                 67.2318     7.0681 11.7258   9.512 7.37e-07 ***
## temp.average                -0.5953     0.7717 12.2194  -0.771    0.455    
## pH.normalized               -2.2687    17.2570 11.4526  -0.131    0.898    
## temp.average:pH.normalized  -0.3574     1.9364 11.8376  -0.185    0.857    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr pH.nrm
## temp.averag -0.977              
## pH.normalzd -0.825  0.819       
## tmp.vrg:pH.  0.796 -0.829 -0.977
summary(null.m)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean ~ (1 | Tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: 426.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.74799 -0.67797 -0.04632  0.52422  2.39231 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept)  3.733   1.932   
##  Residual             35.645   5.970   
## Number of obs: 66, groups:  Tank, 16
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  60.0036     0.8823 13.0885   68.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(all.m)
  mean
Predictors Estimates CI p
(Intercept) 67.23 53.09 – 81.37 <0.001
temp average -0.60 -2.14 – 0.95 0.443
pH normalized -2.27 -36.79 – 32.25 0.896
temp average × pH
normalized 		-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 
tab_model(null.m)
  mean
Predictors Estimates CI p
(Intercept) 60.00 58.24 – 61.77 <0.001
Random Effects
σ2 35.64
Ï„00 Tank 3.73
ICC 0.09
N Tank 16
Observations 66
Marginal R2 / Conditional R2 0.000 / 0.095

no impact.

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

#plot(predictorEffects(null.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
#lines=list(multiline=TRUE))
formal plot both
values_of_temp <- data.frame(temp.average = c(6, 9, 12))

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=all.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 63.65973 2.7344558 58.19363 69.12583
## 2             9           0.0 61.87369 1.5218255 58.83161 64.91578
## 3            12           0.0 60.08766 2.8060427 54.47846 65.69686
## 4             6           0.2 62.77711 1.8291353 59.12072 66.43350
## 5             9           0.2 60.77664 0.9885854 58.80049 62.75279
## 6            12           0.2 58.77617 1.7645888 55.24881 62.30353
## 7             6           0.4 61.89450 1.5818602 58.73241 65.05659
## 8             9           0.4 59.67959 0.8705564 57.93937 61.41981
## 9            12           0.4 57.46468 1.6442557 54.17786 60.75150
## 10            6           0.6 61.01188 2.2241594 56.56585 65.45792
## 11            9           0.6 58.58254 1.2874996 56.00886 61.15621
## 12           12           0.6 56.15319 2.5774571 51.00093 61.30546
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=mean, color=temp.average), size=3) +
  labs(x="pH", y="Mean Shell Coloration", fill = "Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma", guide = FALSE, breaks = c(6, 9, 12))+
  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_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=mean, color=pH.normalized), size=3) +
   labs(x="Temperature (C)", y="Mean Shell Coloration", color = "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("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/color_juv.arctica_linear_both_swap.png", juv.arctica.model.plot, width = 14, height = 4, dpi = 300, bg = 'white')

all 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_color_mod_plot.png", all.model.plot, width = 10, height = 12, dpi = 300, bg = "white")