this script will be the same as 3 except with changes/additions based on Tredennick et al 2021 and their script ‘02-understand-butterflies’. will cast a wide net in this and pair it down later. WILL NEED TO FIX MERCENARIA PLOT IN NEXT.

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.

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(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?

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

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

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

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

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

interaction.v.main<-anova(main_effect.m, interaction.m) #not different, interaction not significant

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

interaction.v.main

## Data: mercenaria
## Models:
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## interaction.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                     
## interaction.m    6 797.02 813.54 -392.51   785.02 0.4926  1     0.4828

main.v.temp<-anova(temp.m, main_effect.m) #ph significant

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

main.v.temp

## Data: mercenaria
## Models:
## temp.m: mean ~ temp.average + (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## temp.m           4 802.38 813.40 -397.19   794.38                        
## main_effect.m    5 795.51 809.28 -392.75   785.51 8.8729  1   0.002894 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

main.v.ph<-anova(ph.m, main_effect.m) #different! so temp term is significant, when it is not there, the model with pH only is significantly different than the one with both.

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

main.v.ph

## Data: mercenaria
## Models:
## ph.m: mean ~ pH.normalized + (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## ph.m             4 805.51 816.53 -398.76   797.51                         
## main_effect.m    5 795.51 809.28 -392.75   785.51 12.004  1  0.0005308 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

null.v.temp<-anova(temp.m, null.m) #different! have temperature alone is signficantly different than a null of just tank.

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

null.v.temp

## Data: mercenaria
## Models:
## null.m: mean ~ (1 | Tank)
## temp.m: mean ~ temp.average + (1 | Tank)
##        npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)   
## null.m    3 809.23 817.49 -401.61   803.23                       
## temp.m    4 802.38 813.40 -397.19   794.38 8.845  1   0.002939 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

null.v.ph<-anova(ph.m, null.m) #different! have ph alone is signficantly different than a null of just tank.

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

null.v.ph

## Data: mercenaria
## Models:
## null.m: mean ~ (1 | Tank)
## ph.m: mean ~ pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## null.m    3 809.23 817.49 -401.61   803.23                       
## ph.m      4 805.51 816.53 -398.76   797.51 5.7138  1    0.01683 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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: 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

summary(main_effect.m)

## Linear mixed model fit by REML ['lmerMod']
## 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 t value
## (Intercept)    102.895      7.165  14.361
## temp.average     2.780      0.721   3.855
## pH.normalized  -19.633      6.334  -3.100
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.936       
## pH.normalzd -0.374  0.092

summary(temp.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ temp.average + (1 | Tank)
##    Data: mercenaria
## 
## REML criterion at convergence: 789.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.14152 -0.56679 -0.03605  0.62725  2.19102 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 47.98    6.926   
##  Residual             41.47    6.439   
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)   94.5478     8.5421  11.068
## temp.average   2.9783     0.9219   3.231
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.976

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

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
σ²	41.54
τ₀₀ _Tank	26.94
ICC	0.39
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.457 / 0.670

tab_model(temp.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	94.55	77.62 – 111.47	<0.001
temp average	2.98	1.15 – 4.80	0.002
Random Effects
σ²	41.47
τ₀₀ _Tank	47.98
ICC	0.54
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.296 / 0.674

tab_model(ph.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	128.71	121.59 – 135.84	<0.001
pH normalized	-22.03	-39.86 – -4.20	0.016
Random Effects
σ²	41.36
τ₀₀ _Tank	60.83
ICC	0.60
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.203 / 0.677

ph and temperature are significant, not interaction. mercenaria color increases (gets lighter) with increasing temperature. mercenaria color decreases (gets darker) with increasing (LESS acidic) ph.

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(interaction.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 = "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_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=interaction.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="Average pH", y="Mean Shell Coloration", fill = "Average 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="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)

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)

#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

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

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

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

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

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

interaction.v.main<-anova(main_effect.m, interaction.m) #not different, interaction not significant

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

interaction.v.main

## Data: mya
## Models:
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
## interaction.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                     
## interaction.m    6 1110.3 1128.3 -549.14   1098.3 2.1824  1     0.1396

main.v.temp<-anova(temp.m, main_effect.m) #ph significant

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

main.v.temp

## Data: mya
## Models:
## temp.m: mean ~ temp.average + (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## temp.m           4 1133.0 1145.0 -562.48   1125.0                         
## main_effect.m    5 1110.5 1125.5 -550.23   1100.5 24.494  1  7.454e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

main.v.ph<-anova(ph.m, main_effect.m) #on edge, but not different

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

main.v.ph

## Data: mya
## Models:
## ph.m: mean ~ pH.normalized + (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## ph.m             4 1111.9 1123.9 -551.94   1103.9                       
## main_effect.m    5 1110.5 1125.5 -550.23   1100.5 3.4071  1    0.06492 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

null.v.temp<-anova(temp.m, null.m) #not different

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

null.v.temp

## Data: mya
## Models:
## null.m: mean ~ (1 | Tank)
## temp.m: mean ~ temp.average + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## null.m    3 1131.3 1140.3 -562.64   1125.3                     
## temp.m    4 1133.0 1145.0 -562.48   1125.0 0.3079  1      0.579

null.v.ph<-anova(ph.m, null.m) #different! have ph alone is signficantly different than a null of just tan k.

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

null.v.ph

## Data: mya
## Models:
## null.m: mean ~ (1 | Tank)
## ph.m: mean ~ pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## null.m    3 1131.3 1140.3 -562.64   1125.3                         
## ph.m      4 1111.9 1123.9 -551.94   1103.9 21.395  1  3.738e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

summary(main_effect.m)

## Linear mixed model fit by REML ['lmerMod']
## 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 t value
## (Intercept)   162.2850     8.2941  19.566
## temp.average   -1.4654     0.8346  -1.756
## pH.normalized -49.9375     7.2668  -6.872
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.937       
## pH.normalzd -0.376  0.097

summary(temp.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ temp.average + (1 | Tank)
##    Data: mya
## 
## REML criterion at convergence: 1117.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.93389 -0.56565  0.01268  0.66083  2.75864 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 177.30   13.315  
##  Residual              77.04    8.777  
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  140.7302    15.9480   8.824
## temp.average  -0.8963     1.7193  -0.521
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.977

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

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
σ²	77.08
τ₀₀ _Tank	34.73
ICC	0.31
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.540 / 0.683

tab_model(temp.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	140.73	109.21 – 172.25	<0.001
temp average	-0.90	-4.29 – 2.50	0.603
Random Effects
σ²	77.04
τ₀₀ _Tank	177.30
ICC	0.70
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.012 / 0.701

tab_model(ph.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	148.63	142.50 – 154.77	<0.001
pH normalized	-48.71	-64.07 – -33.35	<0.001
Random Effects
σ²	77.04
τ₀₀ _Tank	41.39
ICC	0.35
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.510 / 0.681

formal plot

plot(predictorEffects(interaction.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

as pH increases, coloration decreases. i think we want one line as only pH significant

formal plot

plot(predictorEffects(interaction.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 consistant.

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

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=interaction.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="Average pH", y="Mean Shell Coloration", fill = "Average 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="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)

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

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

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

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

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

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

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

interaction.v.main<-anova(main_effect.m, interaction.m) #not different, interaction not significant

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

interaction.v.main

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

main.v.temp<-anova(temp.m, main_effect.m) #ph not significant

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

main.v.temp

## Data: juv.arctica
## Models:
## temp.m: mean ~ temp.average + (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## temp.m           4 433.25 442.01 -212.63   425.25                     
## main_effect.m    5 432.70 443.65 -211.35   422.70 2.5526  1     0.1101

main.v.ph<-anova(ph.m, main_effect.m) #on edge, but not different

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

main.v.ph

## Data: juv.arctica
## Models:
## ph.m: mean ~ pH.normalized + (1 | Tank)
## main_effect.m: mean ~ temp.average + pH.normalized + (1 | Tank)
##               npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## ph.m             4 434.05 442.81 -213.03   426.05                       
## main_effect.m    5 432.70 443.65 -211.35   422.70 3.3491  1    0.06724 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

null.v.temp<-anova(temp.m, null.m) #not different

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

null.v.temp

## Data: juv.arctica
## Models:
## null.m: mean ~ (1 | Tank)
## temp.m: mean ~ temp.average + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## null.m    3 433.84 440.41 -213.92   427.84                     
## temp.m    4 433.25 442.01 -212.63   425.25 2.5879  1     0.1077

null.v.ph<-anova(ph.m, null.m) #not different

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

null.v.ph

## Data: juv.arctica
## Models:
## null.m: mean ~ (1 | Tank)
## ph.m: mean ~ pH.normalized + (1 | Tank)
##        npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## null.m    3 433.84 440.41 -213.92   427.84                     
## ph.m      4 434.05 442.81 -213.03   426.05 1.7915  1     0.1807

main.v.null<-anova(main_effect.m, null.m)

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

main.v.null # all is not different than none

## Data: juv.arctica
## 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 433.84 440.41 -213.92   427.84                       
## main_effect.m    5 432.70 443.65 -211.35   422.70 5.1406  2    0.07651 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(main_effect.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ temp.average + pH.normalized + (1 | Tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: 417.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.75912 -0.71044 -0.03192  0.51743  2.24832 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept)  1.507   1.228   
##  Residual             35.824   5.985   
## Number of obs: 66, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    68.2868     4.0922  16.687
## temp.average   -0.7177     0.4137  -1.735
## pH.normalized  -5.3219     3.5274  -1.509
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.935       
## pH.normalzd -0.368  0.079

summary(temp.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ temp.average + (1 | Tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: 423.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.61193 -0.64895 -0.02313  0.54228  2.27898 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept)  2.12    1.456   
##  Residual             36.08    6.007   
## Number of obs: 66, groups:  Tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)   65.9905     3.9251  16.812
## temp.average  -0.6649     0.4253  -1.563
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.978

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

tab_model(main_effect.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	68.29	60.10 – 76.47	<0.001
temp average	-0.72	-1.54 – 0.11	0.088
pH normalized	-5.32	-12.38 – 1.73	0.137
Random Effects
σ²	35.82
τ₀₀ _Tank	1.51
ICC	0.04
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.078 / 0.116

tab_model(temp.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	65.99	58.14 – 73.84	<0.001
temp average	-0.66	-1.52 – 0.19	0.123
Random Effects
σ²	36.08
τ₀₀ _Tank	2.12
ICC	0.06
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.042 / 0.095

tab_model(ph.m)

	mean
Predictors	Estimates	CI	p
(Intercept)	61.68	58.54 – 64.82	<0.001
pH normalized	-4.94	-12.61 – 2.73	0.203
Random Effects
σ²	35.44
τ₀₀ _Tank	3.49
ICC	0.09
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.032 / 0.119

no impact.

formal plot

plot(predictorEffects(interaction.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

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

effect_interact.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=interaction.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="Average pH", y="Mean Shell Coloration", fill = "Average 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="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)

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')

A-Maine-Tank-2022_All_Color_V4

McMahon

2023-08-19

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

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