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

Set up

read in relevant files

I took the 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

Looking for trends

## `geom_smooth()` using formula = 'y ~ x'
## `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 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)
plot(simulationOutput1)

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(color.1, mercenaria$temp)

plotResiduals(color.1, mercenaria$pH)

#plotResiduals(top_model1, mercenaria$temp)

make models and likelihood ratio test

when talking with Brittany we decided to compare w and w o interaction and then from there test significance of main effects

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

full.v.main<-anova(full.m, main.m, test="Chisq")
full.v.main # no dif - drops AIC, drop interaction

## Data: mercenaria
## Models:
## main.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## full.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## main.m    5 795.51 809.28 -392.75   785.51                     
## full.m    6 797.02 813.54 -392.51   785.02 0.4926  1     0.4828

summary(main.m)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: mercenaria
## 
##      AIC      BIC   logLik deviance df.resid 
##    795.5    809.3   -392.8    785.5      111 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.26256 -0.65541  0.00644  0.59647  2.40795 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 20.35    4.511   
##  Residual             41.62    6.451   
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   102.8706     6.3897  14.1990  16.099 1.62e-10 ***
## temp.average    2.7868     0.6435  14.4332   4.331 0.000646 ***
## pH.normalized -19.5916     5.6590  14.7095  -3.462 0.003574 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.935       
## pH.normalzd -0.374  0.091

confint(main.m)

## Computing profile confidence intervals ...

##                    2.5 %     97.5 %
## .sig01          2.778928   7.279378
## .sigma          5.648080   7.463544
## (Intercept)    89.489723 116.307519
## temp.average    1.429197   4.127774
## pH.normalized -31.489772  -7.787176

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

## Single term deletions using Satterthwaite's method:
## 
## Model:
## mean ~ temp.average + pH.normalized + (1 | Tank)
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## temp.average  780.47  780.47     1 14.433  18.754 0.0006464 ***
## pH.normalized 498.80  498.80     1 14.710  11.986 0.0035737 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(main.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	102.87	90.21 – 115.53	<0.001
temp average	2.79	1.51 – 4.06	<0.001
pH normalized	-19.59	-30.81 – -8.38	0.001
Random Effects
σ²	41.62
τ₀₀ _Tank	20.35
ICC	0.33
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.482 / 0.652

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(full.m, xlevels=list(pH.normalized=seq(0, 0.6, 0.2))),
lines=list(multiline=TRUE))

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

#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=full.m, xlevels=list(temp.average=c(6, 9, 12), pH.normalized=c(0, 0.2, 0.4, 0.6)))

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

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

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

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

##    temp.average pH.normalized      fit       se     lower    upper
## 1             6           0.0 121.4134 3.985039 113.51759 129.3093
## 2             9           0.0 127.8479 2.234808 123.41993 132.2759
## 3            12           0.0 134.2824 4.025759 126.30588 142.2589
## 4             6           0.2 116.4220 2.673760 111.12431 121.7197
## 5             9           0.2 124.0483 1.458269 121.15888 126.9376
## 6            12           0.2 131.6745 2.548673 126.62461 136.7244
## 7             6           0.4 111.4306 2.347622 106.77911 116.0821
## 8             9           0.4 120.2486 1.339659 117.59423 122.9030
## 9            12           0.4 129.0666 2.473412 124.16581 133.9673
## 10            6           0.6 106.4392 3.311728  99.87746 113.0010
## 11            9           0.6 116.4489 1.999697 112.48679 120.4111
## 12           12           0.6 126.4586 3.882374 118.76621 134.1511

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 = "none", breaks = c(6, 9, 12))+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)+
  #scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica"),
    axis.title.y = element_text(size=15,family = "Helvetica"),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp

effect_temp_swap <- ggplot() + 
  #geom_point(data=effect_interact.both, aes(x=temp.average, y=fit), color="black") +
   geom_line(data = effect_interact.both, aes(x = temp.average, y = fit, group = pH.normalized, color = pH.normalized), alpha = 0.7) +
   geom_ribbon(data= effect_interact.both, aes(x=temp.average, ymin=lower, ymax=upper, group = pH.normalized, fill=pH.normalized), alpha= 0.3) +
   geom_point(data=mercenaria, mapping = aes(x=temp.average, y=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 = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
    axis.title.y = element_text(size=15,family = "Helvetica",hjust = 0.5),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp_swap

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

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

ggsave("02_output/02_plots/color_mercenaria_linear.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(simulationOutput1)

shapiro.test(resid(color.1))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(color.1)
## W = 0.98572, p-value = 0.1249

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(color.1, mya$temp)

plotResiduals(color.1, mya$pH)

#plotResiduals(top_model1, mercenaria$temp)

make models and likelihood ratio test

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

full.v.main<-anova(full.m, main.m, test="Chisq")
full.v.main # no dif - drops AIC, drop interaction

## Data: mya
## Models:
## main.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## full.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## main.m    5 1110.5 1125.5 -550.23   1100.5                     
## full.m    6 1110.3 1128.3 -549.14   1098.3 2.1824  1     0.1396

summary(main.m)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: mya
## 
##      AIC      BIC   logLik deviance df.resid 
##   1110.5   1125.5   -550.2   1100.5      145 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.97320 -0.56490  0.03388  0.74850  2.71712 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 26.58    5.156   
##  Residual             77.10    8.781   
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   162.3229     7.4692  15.8188  21.732 3.33e-13 ***
## temp.average   -1.4706     0.7521  16.0762  -1.955   0.0682 .  
## pH.normalized -49.9321     6.5427  15.7431  -7.632 1.13e-06 ***
## ---
## 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

confint(main.m)

## Computing profile confidence intervals ...

##                    2.5 %      97.5 %
## .sig01          3.179168   8.2698850
## .sigma          7.825115   9.9480662
## (Intercept)   146.706365 177.8530462
## temp.average   -3.031654   0.1021498
## pH.normalized -63.582793 -36.2936394

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

## Single term deletions using Satterthwaite's method:
## 
## Model:
## mean ~ temp.average + pH.normalized + (1 | Tank)
##               Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## temp.average   294.7   294.7     1 16.076  3.8228   0.06817 .  
## pH.normalized 4490.7  4490.7     1 15.743 58.2439 1.127e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(main.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	162.32	147.56 – 177.09	<0.001
temp average	-1.47	-2.96 – 0.02	0.052
pH normalized	-49.93	-62.86 – -37.00	<0.001
Random Effects
σ²	77.10
τ₀₀ _Tank	26.58
ICC	0.26
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.559 / 0.672

formal plot

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

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

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=full.m, xlevels=list(temp.average=c(6, 9, 12), pH.normalized=c(0, 0.2, 0.4, 0.6)))

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

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

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

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

##    temp.average pH.normalized      fit       se    lower    upper
## 1             6           0.0 149.1485 4.384107 140.4839 157.8130
## 2             9           0.0 149.3366 2.433678 144.5268 154.1464
## 3            12           0.0 149.5247 4.410052 140.8090 158.2405
## 4             6           0.2 141.7378 2.939691 135.9279 147.5476
## 5             9           0.2 139.0711 1.582096 135.9443 142.1978
## 6            12           0.2 136.4044 2.790119 130.8901 141.9186
## 7             6           0.4 134.3271 2.575520 129.2369 139.4172
## 8             9           0.4 128.8055 1.450958 125.9380 131.6731
## 9            12           0.4 123.2840 2.720211 117.9080 128.6601
## 10            6           0.6 126.9164 3.632485 119.7373 134.0954
## 11            9           0.6 118.5400 2.174752 114.2420 122.8381
## 12           12           0.6 110.1637 4.277022 101.7108 118.6166

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 = "none", breaks = c(6, 9, 12))+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)+
  scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica"),
    axis.title.y = element_text(size=15,family = "Helvetica"),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

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 = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
    axis.title.y = element_text(size=15,family = "Helvetica",hjust = 0.5),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp_swap

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

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

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

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(simulationOutput1)

shapiro.test(resid(color.1))

## 
##  Shapiro-Wilk normality test
## 
## data:  resid(color.1)
## W = 0.98992, p-value = 0.8723

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(color.1, juv.arctica$temp)

plotResiduals(color.1, juv.arctica$pH)

#plotResiduals(top_model1, mercenaria$temp)

make models and likelihood ratio test

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

## boundary (singular) fit: see help('isSingular')

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

## boundary (singular) fit: see help('isSingular')

full.v.main<-anova(full.m, main.m, test="Chisq")
full.v.main # no dif - drops AIC, drop interaction

## Data: juv.arctica
## Models:
## main.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## full.m: mean ~ temp.average * pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## main.m    5 432.70 443.65 -211.35   422.70                     
## full.m    6 434.68 447.81 -211.34   422.68 0.0247  1     0.8752

summary(main.m)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: juv.arctica
## 
##      AIC      BIC   logLik deviance df.resid 
##    432.7    443.6   -211.4    422.7       61 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.94750 -0.70080 -0.03117  0.62721  2.35086 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept)  0.0     0.00    
##  Residual             35.4     5.95    
## Number of obs: 66, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)    68.3168     3.7643 66.0000  18.149   <2e-16 ***
## temp.average   -0.7258     0.3809 66.0000  -1.905   0.0611 .  
## pH.normalized  -5.2104     3.2252 66.0000  -1.616   0.1110    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.935       
## pH.normalzd -0.368  0.078
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

confint(main.m)

## Computing profile confidence intervals ...

##                    2.5 %      97.5 %
## .sig01          0.000000  3.21096076
## .sigma          4.978291  7.12927818
## (Intercept)    60.619899 75.94393333
## temp.average   -1.491953  0.05864277
## pH.normalized -11.978240  1.27929098

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

## Single term deletions using Satterthwaite's method:
## 
## Model:
## mean ~ temp.average + pH.normalized + (1 | Tank)
##                Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## temp.average  128.495 128.495     1    66  3.6299 0.06111 .
## pH.normalized  92.391  92.391     1    66  2.6100 0.11096  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

tab_model(main.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	68.32	60.79 – 75.84	<0.001
temp average	-0.73	-1.49 – 0.04	0.061
pH normalized	-5.21	-11.66 – 1.24	0.111
Random Effects
σ²	35.40
τ₀₀ _Tank	0.00
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.082 / NA

formal plot

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

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

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=full.m, xlevels=list(temp.average=c(6, 9, 12), pH.normalized=c(0, 0.2, 0.4, 0.6)))

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

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

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

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

##    temp.average pH.normalized      fit        se    lower    upper
## 1             6           0.0 63.71382 2.4083082 58.89968 68.52796
## 2             9           0.0 61.80431 1.3428174 59.12006 64.48857
## 3            12           0.0 59.89480 2.4967872 54.90380 64.88581
## 4             6           0.2 62.81662 1.6103072 59.59766 66.03558
## 5             9           0.2 60.74562 0.8723342 59.00185 62.48940
## 6            12           0.2 58.67463 1.5675285 55.54118 61.80807
## 7             6           0.4 61.91942 1.3855884 59.14966 64.68917
## 8             9           0.4 59.68693 0.7569363 58.17384 61.20003
## 9            12           0.4 57.45445 1.4434514 54.56903 60.33987
## 10            6           0.6 61.02221 1.9442855 57.13564 64.90878
## 11            9           0.6 58.62824 1.1131555 56.40308 60.85341
## 12           12           0.6 56.23428 2.2612310 51.71414 60.75441

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 = "none", breaks = c(6, 9, 12))+
  scale_fill_viridis( option= "plasma")+
  scale_x_continuous(breaks = original_values, labels = custom_labels_x)+
  scale_fill_viridis( option= "cividis", guide = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica"),
    axis.title.y = element_text(size=15,family = "Helvetica"),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )

## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.

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 = "none", direction = -1)+
  theme(
    axis.text.x = element_text(size= 10, family = "Helvetica"),
    axis.text.y = element_text(size= 10, family = "Helvetica"),
    axis.title.x = element_text(size=15,family = "Helvetica",hjust = 0.5),
    axis.title.y = element_text(size=15,family = "Helvetica",hjust = 0.5),
    legend.title = element_text(size = 15,family = "Helvetica"),
    legend.text = element_text(size = 10,family = "Helvetica")
  )
effect_temp_swap

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

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

ggsave("02_output/02_plots/color_juv.arctica_linear.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 = 5, height = 11, dpi = 300, bg = "white")

A-Maine-Tank-2022_All_Color_V6

McMahon

2023-09-27

Set up

read in relevant files

remove dead and lost

look at data

Looking for trends

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

Plot

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

formal plot both

Plot

mya

check structure

linear model

‘normal’ linear mixed effect model

check the normality assumptions of full model

make models and likelihood ratio test

formal plot

formal plot both

Plot

juv arctica

check structure

linear model

‘normal’ linear mixed effect model

check the normality assumptions of full model

make models and likelihood ratio test

formal plot

formal plot both

all plots