Set up

read in relevant files

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

remove dead and lost

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

## [1] 425  18

## [1] 366  18

look at data

## Looking for trends

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

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

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

## [1] 7.39

## [1] 6.14

Plot

## ------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:ggpubr':
## 
##     mutate

##    temp  pH  N      mean       sd       se        ci
## 1     6 7.4 25 124.31220 30.77729 6.155459 10.531264
## 2     6 7.6 27 111.38259 27.62450 5.316337  9.067639
## 3     6 7.8 24 102.63350 25.12788 5.129208  8.790803
## 4     6   8 27 107.49989 24.26398 4.669605  7.964562
## 5     9 7.4 46 127.29761 34.42511 5.075706  8.524280
## 6     9 7.6 45 114.15244 29.07736 4.334596  7.283119
## 7     9 7.8 42 107.31557 26.61527 4.106825  6.911286
## 8     9   8 51 105.50114 26.45628 3.704620  6.208592
## 9    12 7.4 18 123.82133 29.88144 7.043122 12.252263
## 10   12 7.6 22 122.27141 35.53425 7.575928 13.036225
## 11   12 7.8 21  99.05114 26.22239 5.722195  9.869175
## 12   12   8 18 110.24456 29.29434 6.904743 12.011537

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

Individual Species Groups

mercenaria

check structure

Looks pretty normal to me!

## 'data.frame':    116 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ File path      : chr  NA "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\mercenaria\\4th round\\b41.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\mercenaria\\b38.jpeg" NA ...
##  $ file name      : chr  "b38.jpeg" "b42" "b39" "b47" ...
##  $ shell          : chr  "b38" "b42" "b39" "b47" ...
##  $ N              : num  456633 274746 316471 329233 343592 ...
##  $ mean           : num  146 140 142 126 138 ...
##  $ stdev          : num  24.1 16.8 22.7 20.9 21.1 ...
##  $ mode           : num  139 144 128 134 145 126 133 137 130 127 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 4 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 2 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "mercenaria" "mercenaria" "mercenaria" "mercenaria" ...
##  $ temp.average   : num  12.1 12.1 12.1 12.1 12.1 ...
##  $ temp.stdev     : num  0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.63 ...
##  $ pH.average.YSI : num  7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 8 ...
##  $ pH.stdev.YSI   : num  0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.07 ...
##  $ pH.average     : num  7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 8 ...
##  $ pH.stdev       : num  0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.05 ...
##  $ pH.normalized  : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 3.05 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mercenaria$mean
## W = 0.9899, p-value = 0.5504

##    vars   n   mean    sd median trimmed   mad   min    max range  skew kurtosis
## X1    1 116 121.68 11.03 122.41  121.71 12.29 96.69 146.08 49.39 -0.04    -0.67
##      se
## X1 1.02

linear model

‘normal’ linear mixed effect model

first we will try this and see how it looks as per Brittany’s suggestion. This will not specify a family, most similar to what I did way back when with CCA.

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

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

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

simulationOutput1<-simulateResiduals(color.1)

dd1<-dredge(color.1, rank=BIC) #BIC! more conservative

## Warning in dredge(color.1, rank = BIC): comparing models fitted by REML

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = mercenaria, na.action = na.fail)
## ---
## Model selection table 
##    (Int) pH.nrm tmp.avr pH.nrm:tmp.avr df   logLik   BIC delta weight
## 4 102.90 -19.63   2.779                 5 -388.293 800.4  0.00  0.510
## 8 108.60 -37.04   2.132          1.995  6 -386.010 800.5  0.19  0.465
## 2 128.70 -22.03                         4 -394.073 807.2  6.81  0.017
## 3  94.55          2.978                 4 -394.895 808.8  8.45  0.007
## 1 121.40                                3 -399.851 814.0 13.61  0.001
## Models ranked by BIC(x) 
## Random terms (all models): 
##   1 | Tank

top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ pH.normalized + temp.average + (1 | Tank)
##    Data: mercenaria
## REML criterion at convergence: 776.5863
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 5.190   
##  Residual             6.445   
## Number of obs: 116, groups:  Tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       102.895        -19.633          2.779

simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

#plot(top_model1)
tab_model(top_model1)

	mean
Predictors	Estimates	CI	p
(Intercept)	102.90	88.70 – 117.09	<0.001
pH normalized	-19.63	-32.18 – -7.08	0.002
temp average	2.78	1.35 – 4.21	<0.001
Random Effects
σ²	41.54
τ₀₀ _Tank	26.94
ICC	0.39
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.457 / 0.670

testZeroInflation(top_model1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(top_model1, mercenaria$pH)

plotResiduals(top_model1, mercenaria$temp)

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp & pH on Mercenaria Coloration")

	Effect of Temp & pH on Mercenaria Coloration
Predictors	Estimates	CI	p
(Intercept)	102.90	88.70 – 117.09	<0.001
pH (normalized)	-19.63	-32.18 – -7.08	0.002
Temp	2.78	1.35 – 4.21	<0.001
Random Effects
σ²	41.54
τ₀₀ _Tank	26.94
ICC	0.39
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.457 / 0.670

formal plot

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

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

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=top_model1, at=values_of_pH)

## NOTE: temp.averagepH.normalized does not appear in the model

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

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

effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = 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")+
  theme_classic()+
  scale_color_viridis(direction = -1)+
  scale_fill_viridis(direction = -1)
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 = 5, height = 4, dpi = 300)

poisson

Warnings for looks horrible, def not right.

## Fixed term is "(Intercept)"

## Global model call: glmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = mercenaria, family = poisson, na.action = "na.fail")
## ---
## Model selection table 
##   (Int)  pH.nrm tmp.avr pH.nrm:tmp.avr df logLik AIC
## 1 4.796                                 2   -Inf Inf
## 2 4.857 -0.1853                         3   -Inf Inf
## 3 4.572         0.02472                 3   -Inf Inf
## 4 4.643 -0.1660 0.02289                 4   -Inf Inf
## 8 4.701 -0.3420 0.01640        0.02007  5   -Inf Inf
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: mean ~ (1 | Tank)
##    Data: mercenaria
##      AIC      BIC   logLik deviance df.resid 
##      Inf      Inf     -Inf      Inf      114 
## Random effects:
##  Groups Name        Std.Dev.
##  Tank   (Intercept) 1       
## Number of obs: 116, groups:  Tank, 16
## Fixed Effects:
## (Intercept)  
##       4.796  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 6148 optimizer warnings; 1 lme4 warnings

	mean
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	121.01	74.11 – 197.59	<0.001
Random Effects
σ²	0.01
τ₀₀ _Tank	1.00
ICC	0.99
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.000 / 0.992

gamma distribution?

gamma requires all positive values.

This looks pretty bad, note that its likely not right because these are quite normal dist.

color.1 <- glmer(mean~temp.average*pH.normalized+(1|Tank), data = mercenaria, family = Gamma(link = "log"),  na.action = "na.fail")

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00292039 (tol = 0.002, component 1)

simulationOutput1<-simulateResiduals(color.1)
plot(simulationOutput1) #wacky red lines

dd1<-dredge(color.1, rank=AIC)

## Fixed term is "(Intercept)"

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00301242 (tol = 0.002, component 1)

dd1

## Global model call: glmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = mercenaria, family = Gamma(link = "log"), na.action = "na.fail")
## ---
## Model selection table 
##   (Int)  pH.nrm tmp.avr pH.nrm:tmp.avr df   logLik   AIC delta weight
## 4 4.643 -0.1661 0.02262                 5 -386.630 783.3  0.00  0.497
## 8 4.705 -0.3527 0.01582        0.02109  6 -386.434 784.9  1.61  0.222
## 3 4.568         0.02462                 4 -388.869 785.7  2.48  0.144
## 2 4.856 -0.1910                         4 -389.367 786.7  3.47  0.087
## 1 4.792                                 3 -390.934 787.9  4.61  0.050
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean ~ pH.normalized + temp.average + (1 | Tank)
##    Data: mercenaria
##       AIC       BIC    logLik  deviance  df.resid 
##  783.2599  797.0278 -386.6299  773.2599       111 
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 0.03029 
##  Residual             0.05497 
## Number of obs: 116, groups:  Tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       4.64319       -0.16608        0.02262

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean ~ pH.normalized + temp.average + (1 | Tank)
##    Data: mercenaria
## 
##      AIC      BIC   logLik deviance df.resid 
##    783.3    797.0   -386.6    773.3      111 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2591 -0.5859  0.0523  0.6655  2.4971 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  Tank     (Intercept) 0.0009174 0.03029 
##  Residual             0.0030222 0.05497 
## Number of obs: 116, groups:  Tank, 16
## 
## Fixed effects:
##                Estimate Std. Error t value Pr(>|z|)    
## (Intercept)    4.643187   0.078375  59.243  < 2e-16 ***
## pH.normalized -0.166084   0.067905  -2.446  0.01445 *  
## temp.average   0.022624   0.007809   2.897  0.00377 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pH.nrm
## pH.normalzd -0.378       
## temp.averag -0.940  0.107

simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

testDispersion(simulationOutput.top)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.00011035, p-value < 2.2e-16
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	mean
Predictors	Estimates	CI	p
(Intercept)	103.87	88.93 – 121.33	<0.001
pH normalized	0.85	0.74 – 0.97	0.016
temp average	1.02	1.01 – 1.04	0.005
Random Effects
σ²	0.00
τ₀₀ _Tank	0.00
ICC	0.23
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.500 / 0.616

testZeroInflation(top_model1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(top_model1, mercenaria$ph)

plotResiduals(top_model1, mercenaria$temp)

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp and pH on Mercenaria Growth")

	Effect of Temp and pH on Mercenaria Growth
Predictors	Estimates	CI	p
(Intercept)	103.87	88.93 – 121.33	<0.001
pH (normalized)	0.85	0.74 – 0.97	0.016
Temp	1.02	1.01 – 1.04	0.005
Random Effects
σ²	0.00
τ₀₀ _Tank	0.00
ICC	0.23
N _Tank	16
Observations	116
Marginal R² / Conditional R²	0.500 / 0.616

Plot

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

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

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

## 'data.frame':    150 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r41.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r40.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r37.jpeg" "C:\\Users\\thatc\\OneDrive - Iowa State University\\Tank Experiment_Boron Isotopes Paper\\Mya arenaria pics\\jpg\\r42.jpeg" ...
##  $ file name      : chr  "r41.jpeg" "r40.jpeg" "r37.jpeg" "r42.jpeg" ...
##  $ shell          : chr  "r41" "r40" "r37" "r42" ...
##  $ N              : num  661697 587852 417483 603232 286395 ...
##  $ mean           : num  130 140 158 139 118 ...
##  $ stdev          : num  45.9 46.1 39.6 46.9 51.9 ...
##  $ mode           : num  127 171 174 96 180 172 117 181 170 112 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 4 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 2 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "mya" "mya" "mya" "mya" ...
##  $ temp.average   : num  12.1 12.1 12.1 12.1 12.1 ...
##  $ temp.stdev     : num  0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.42 0.63 ...
##  $ pH.average.YSI : num  7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 8 ...
##  $ pH.stdev.YSI   : num  0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.07 ...
##  $ pH.average     : num  7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 8 ...
##  $ pH.stdev       : num  0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.05 ...
##  $ pH.normalized  : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 5.92 3.05 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mya$mean
## W = 0.98704, p-value = 0.1759

##    vars   n   mean    sd median trimmed  mad   min    max range skew kurtosis
## X1    1 150 132.61 15.32 129.99  132.29 15.1 96.47 172.74 76.27 0.21    -0.38
##      se
## X1 1.25

linear model

‘normal’ linear mixed effect model

first we will try this and see how it looks as per Brittany’s suggestion. This will not specify a family, most similar to what I did way back when with CCA.

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

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

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

simulationOutput1<-simulateResiduals(color.1)

dd1<-dredge(color.1, rank=AIC) #BIC = more conservative

## Warning in dredge(color.1, rank = AIC): comparing models fitted by REML

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = mean ~ temp.average + pH.normalized + (1 | Tank), 
##     data = mya, na.action = na.fail)
## ---
## Model selection table 
##   (Int) pH.nrm tmp.avr df   logLik    AIC delta weight
## 4 162.3 -49.94 -1.4650  5 -545.348 1100.7  0.00  0.769
## 2 148.6 -48.71          4 -547.552 1103.1  2.41  0.231
## 3 140.7        -0.8963  4 -558.944 1125.9 25.19  0.000
## 1 132.6                 3 -560.532 1127.1 26.37  0.000
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1

## Linear mixed model fit by REML ['lmerMod']
## Formula: mean ~ pH.normalized + temp.average + (1 | Tank)
##    Data: mya
## REML criterion at convergence: 1090.695
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 5.893   
##  Residual             8.780   
## Number of obs: 150, groups:  Tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       162.285        -49.937         -1.465

simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

#plot(top_model1)
tab_model(top_model1)

	mean
Predictors	Estimates	CI	p
(Intercept)	162.28	145.89 – 178.68	<0.001
pH normalized	-49.94	-64.30 – -35.57	<0.001
temp average	-1.47	-3.11 – 0.18	0.081
Random Effects
σ²	77.08
τ₀₀ _Tank	34.73
ICC	0.31
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.540 / 0.683

testZeroInflation(top_model1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(top_model1, mya$pH)

plotResiduals(top_model1, mya$temp)

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  #pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp & pH on Mya Coloration")

	Effect of Temp & pH on Mya Coloration
Predictors	Estimates	CI	p
(Intercept)	162.28	145.89 – 178.68	<0.001
pH normalized	-49.94	-64.30 – -35.57	<0.001
temp average	-1.47	-3.11 – 0.18	0.081
Random Effects
σ²	77.08
τ₀₀ _Tank	34.73
ICC	0.31
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.540 / 0.683

formal plot

i think we want one line as only pH significant

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

effect_top_model<-effects::effect(term= c("pH.normalized"), mod=top_model1)
effect_top_model<-as.data.frame(effect_top_model)
effect_top_model

##   pH.normalized      fit       se    lower    upper
## 1           0.0 149.2078 2.909247 143.4585 154.9572
## 2           0.2 139.2203 1.896676 135.4721 142.9686
## 3           0.3 134.2266 1.657836 130.9503 137.5029
## 4           0.5 124.2391 2.054001 120.1799 128.2983
## 5           0.7 114.2516 3.148646 108.0291 120.4741

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

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

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

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized"), mod=top_model1, at=values_of_temp)

## NOTE: temp.averagepH.normalized does not appear in the model

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

##    temp.average pH.normalized      fit       se    lower    upper
## 1           6.1           0.0 153.3458 3.944119 145.5513 161.1403
## 2           7.6           0.0 151.1477 3.225761 144.7728 157.5225
## 3           9.1           0.0 148.9495 2.896430 143.2255 154.6735
## 4          11.0           0.0 146.1652 3.214795 139.8120 152.5184
## 5          12.0           0.0 144.6998 3.657704 137.4713 151.9282
## 6           6.1           0.2 143.3583 3.166702 137.1002 149.6165
## 7           7.6           0.2 141.1602 2.286697 136.6411 145.6792
## 8           9.1           0.2 138.9620 1.887930 135.2310 142.6930
## 9          11.0           0.2 136.1777 2.440316 131.3551 141.0003
## 10         12.0           0.2 134.7123 3.038648 128.7072 140.7173
## 11          6.1           0.3 138.3646 2.974588 132.4861 144.2431
## 12          7.6           0.3 136.1664 2.055460 132.1044 140.2285
## 13          9.1           0.3 133.9683 1.654077 130.6994 137.2371
## 14         11.0           0.3 131.1840 2.312891 126.6131 135.7548
## 15         12.0           0.3 129.7185 2.957174 123.8745 135.5626
## 16          6.1           0.5 128.3771 3.107561 122.2358 134.5184
## 17          7.6           0.5 126.1789 2.320603 121.5929 130.7650
## 18          9.1           0.5 123.9808 2.061014 119.9077 128.0538
## 19         11.0           0.5 121.1965 2.703059 115.8546 126.5383
## 20         12.0           0.5 119.7310 3.307007 113.1956 126.2664
## 21          6.1           0.7 118.3896 3.832787 110.8151 125.9641
## 22          7.6           0.7 116.1914 3.281773 109.7059 122.6770
## 23          9.1           0.7 113.9933 3.159769 107.7488 120.2377
## 24         11.0           0.7 111.2090 3.672622 103.9510 118.4669
## 25         12.0           0.7 109.7435 4.165605 101.5113 117.9757

effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.both, aes(x=pH.normalized, ymin=lower, ymax=upper, group = temp.average, fill=temp.average), alpha= 0.3) +
  geom_point(data=mercenaria, mapping = aes(x=pH.normalized, y=mean, color=temp.average), size=3) +
  labs(x="Average pH (Normalized)", y="Mean Shell Coloration")+
  theme_classic()+
  scale_color_viridis(direction = -1)+
  scale_fill_viridis(direction = -1, guide=FALSE)
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.

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

poisson

Warnings looks horrible, def not right.

## Fixed term is "(Intercept)"

## Global model call: glmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = mya, family = poisson, na.action = "na.fail")
## ---
## Model selection table 
##   (Int)   pH.nrm    tmp.avr pH.nrm:tmp.avr df logLik AIC
## 1 4.883                                     2   -Inf Inf
## 2 5.002 -0.36260                            3   -Inf Inf
## 3 4.953          -0.0077470                 3   -Inf Inf
## 4 5.113 -0.37260 -0.0118600                 4   -Inf Inf
## 8 5.005 -0.04327  0.0002832       -0.03755  5   -Inf Inf
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: mean ~ (1 | Tank)
##    Data: mya
##      AIC      BIC   logLik deviance df.resid 
##      Inf      Inf     -Inf      Inf      148 
## Random effects:
##  Groups Name        Std.Dev.
##  Tank   (Intercept) 1       
## Number of obs: 150, groups:  Tank, 16
## Fixed Effects:
## (Intercept)  
##       4.883  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 7950 optimizer warnings; 1 lme4 warnings

	mean
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	131.98	80.84 – 215.47	<0.001
Random Effects
σ²	0.01
τ₀₀ _Tank	1.00
ICC	0.99
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.000 / 0.993

gamma distribution?

gamma requires all positive values.

This looks pretty bad, note that its likely not right because these are quite normal dist.

color.1 <- glmer(mean~temp.average*pH.normalized+(1|Tank), data = mya, family = Gamma(link = "log"),  na.action = "na.fail")

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0102884 (tol = 0.002, component 1)

simulationOutput1<-simulateResiduals(color.1)
plot(simulationOutput1) #wacky red lines

dd1<-dredge(color.1, rank=AIC)

## Fixed term is "(Intercept)"

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00609365 (tol = 0.002, component 1)

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0102884 (tol = 0.002, component 1)

dd1

## Global model call: glmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = mya, family = Gamma(link = "log"), na.action = "na.fail")
## ---
## Model selection table 
##   (Int)   pH.nrm    tmp.avr pH.nrm:tmp.avr df   logLik    AIC delta weight
## 2 5.002 -0.36460                            4 -542.409 1092.8  0.00  0.422
## 4 5.112 -0.37390 -1.175e-02                 5 -541.588 1093.2  0.36  0.353
## 8 5.006 -0.05006  5.384e-05       -0.03664  6 -541.101 1094.2  1.38  0.211
## 1 4.882                                     3 -547.222 1100.4  7.63  0.009
## 3 4.951          -7.523e-03                 4 -547.126 1102.3  9.43  0.004
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean ~ pH.normalized + (1 | Tank)
##    Data: mya
##       AIC       BIC    logLik  deviance  df.resid 
## 1092.8175 1104.8601 -542.4088 1084.8175       146 
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 0.03621 
##  Residual             0.06695 
## Number of obs: 150, groups:  Tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized  
##        5.0017        -0.3646

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean ~ pH.normalized + (1 | Tank)
##    Data: mya
## 
##      AIC      BIC   logLik deviance df.resid 
##   1092.8   1104.9   -542.4   1084.8      146 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6130 -0.5833  0.0710  0.7334  2.4021 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 0.001311 0.03621 
##  Residual             0.004483 0.06695 
## Number of obs: 150, groups:  Tank, 16
## 
## Fixed effects:
##               Estimate Std. Error t value Pr(>|z|)    
## (Intercept)    5.00172    0.03382 147.906  < 2e-16 ***
## pH.normalized -0.36457    0.08491  -4.294 1.76e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## pH.normalzd -0.823

simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

testDispersion(simulationOutput.top)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.00021708, p-value < 2.2e-16
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	mean
Predictors	Estimates	CI	p
(Intercept)	148.67	139.06 – 158.94	<0.001
pH normalized	0.69	0.59 – 0.82	<0.001
Random Effects
σ²	0.00
τ₀₀ _Tank	0.00
ICC	0.23
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.544 / 0.647

testZeroInflation(top_model1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(top_model1, mya$ph)

plotResiduals(top_model1, mya$temp)

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp and pH on mya Growth")

## Length of `pred.labels` does not equal number of predictors, no labelling applied.

	Effect of Temp and pH on mya Growth
Predictors	Estimates	CI	p
(Intercept)	148.67	139.06 – 158.94	<0.001
pH.normalized	0.69	0.59 – 0.82	<0.001
Random Effects
σ²	0.00
τ₀₀ _Tank	0.00
ICC	0.23
N _Tank	16
Observations	150
Marginal R² / Conditional R²	0.544 / 0.647

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

scallop

check structure

Outliers, white shells - no point really modeling as we know this shouldnt relate to pH/temp

## 'data.frame':    34 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 2 2 2 3 3 3 4 4 5 5 ...
##  $ File path      : chr  NA NA "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\scallops\\r55.jpg" NA ...
##  $ file name      : chr  NA NA "r55.jpg" NA ...
##  $ shell          : chr  "r67" "r34" "r55" "r39" ...
##  $ N              : num  1239933 2087642 2733522 920068 1839238 ...
##  $ mean           : num  96.6 83.1 104.9 82 80.8 ...
##  $ stdev          : num  32.2 22.8 41.1 33.7 23.1 ...
##  $ mode           : num  91 77 117 73 74 143 101 78 132 100 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 4 4 4 4 4 4 4 4 1 1 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 2 2 2 1 1 1 2 2 3 3 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "scallop" "scallop" "scallop" "scallop" ...
##  $ temp.average   : num  9.19 9.19 9.19 6.26 6.26 ...
##  $ temp.stdev     : num  0.63 0.63 0.63 0.78 0.78 0.78 0.33 0.33 0.75 0.75 ...
##  $ pH.average.YSI : num  8 8 8 8.02 8.02 8.02 8 8 7.44 7.44 ...
##  $ pH.stdev.YSI   : num  0.07 0.07 0.07 0.06 0.06 0.06 0.06 0.06 0.07 0.07 ...
##  $ pH.average     : num  8 8 8 8.05 8.05 8.05 8.01 8.01 7.4 7.4 ...
##  $ pH.stdev       : num  0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.06 0.06 ...
##  $ pH.normalized  : num  0.61 0.61 0.61 0.66 0.66 ...
##  $ temp.normalized: num  3.05 3.05 3.05 0.12 0.12 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  scallop$mean
## W = 0.87395, p-value = 0.001008

##    vars  n  mean    sd median trimmed   mad   min   max range skew kurtosis
## X1    1 34 99.42 21.87  97.03   96.93 17.61 70.42 150.9 80.48 1.05     0.25
##      se
## X1 3.75

Plot

These plots show that there may be something going on with pH but we have poor coverage and large error due to the small sample size. This is also only data with the four outliers removed.

##    temp  pH N      mean        sd        se         ci
## 1     6 7.4 2 130.04750 28.215682 19.951500 125.968813
## 2     6 7.6 4  91.24725 12.042746  6.021373  14.170479
## 3     6 7.8 3  86.88233 10.955324  6.325059  18.469082
## 4     6   8 3 102.01233 35.680181 20.599962  60.151592
## 5     9 7.4 4 124.94100 27.231791 13.615895  32.043150
## 6     9 7.6 4  86.78575 10.831082  5.415541  12.744736
## 7     9 7.8 4  88.52925 13.911565  6.955782  16.369484
## 8     9   8 5  97.43140  8.727595  3.903099   8.320809
## 9    12 7.4 3 113.93433 14.577130  8.416110  24.574920
## 10   12 7.8 2  79.28850  8.353052  5.906500  37.292173

juv arctica

check structure

Appears pretty normal. Lower than others.

## 'data.frame':    66 obs. of  20 variables:
##  $ Tank           : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ File path      : chr  "C:\\Users\\thatc\\OneDrive - Iowa State University\\Desktop\\jpg\\juv Arctica\\g29.jpeg" NA NA NA ...
##  $ file name      : chr  "g29.jpeg" NA NA NA ...
##  $ shell          : chr  "g29" "g30" "g31" "g32" ...
##  $ N              : num  381123 339242 764948 602176 970122 ...
##  $ mean           : num  52.5 41 58.5 52.7 68.8 ...
##  $ stdev          : num  9.96 17.64 19.74 16.74 24.64 ...
##  $ mode           : num  47 28 46 51 68 47 46 50 47 58 ...
##  $ pH             : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 4 4 4 4 4 4 ...
##  $ temp           : Factor w/ 3 levels "6","9","12": 3 3 3 3 2 2 2 2 1 1 ...
##  $ dead?          : chr  NA NA NA NA ...
##  $ species        : chr  "juv.arctica" "juv.arctica" "juv.arctica" "juv.arctica" ...
##  $ temp.average   : num  12.06 12.06 12.06 12.06 9.19 ...
##  $ temp.stdev     : num  0.42 0.42 0.42 0.42 0.63 0.63 0.63 0.63 0.78 0.78 ...
##  $ pH.average.YSI : num  7.57 7.57 7.57 7.57 8 8 8 8 8.02 8.02 ...
##  $ pH.stdev.YSI   : num  0.04 0.04 0.04 0.04 0.07 0.07 0.07 0.07 0.06 0.06 ...
##  $ pH.average     : num  7.59 7.59 7.59 7.59 8 8 8 8 8.05 8.05 ...
##  $ pH.stdev       : num  0.07 0.07 0.07 0.07 0.05 0.05 0.05 0.05 0.05 0.05 ...
##  $ pH.normalized  : num  0.2 0.2 0.2 0.2 0.61 ...
##  $ temp.normalized: num  5.92 5.92 5.92 5.92 3.05 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  juv.arctica$mean
## W = 0.9902, p-value = 0.8847

##    vars  n  mean   sd median trimmed  mad   min   max range  skew kurtosis   se
## X1    1 66 59.97 6.25  59.67   59.92 5.99 40.98 75.74 34.75 -0.03     0.43 0.77

linear model

‘normal’ linear mixed effect model

first we will try this and see how it looks as per Brittany’s suggestion. This will not specify a family, most similar to what I did way back when with CCA.

best fit with AIC includes temp and ph and interaction but none are significant similar to what I saw above, but makes sense here as we do not expect a relationship.

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

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

simulationOutput1<-simulateResiduals(color.1)

dd1<-dredge(color.1, rank=AIC) #BIC = more conservative

## Warning in dredge(color.1, rank = AIC): comparing models fitted by REML

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = juv.arctica, na.action = na.fail)
## ---
## Model selection table 
##   (Int) pH.nrm tmp.avr pH.nrm:tmp.avr df   logLik   AIC delta weight
## 8 67.23 -2.269 -0.5953        -0.3574  6 -207.039 426.1  0.00  0.494
## 4 68.29 -5.322 -0.7177                 5 -208.614 427.2  1.15  0.278
## 2 61.68 -4.936                         4 -210.057 428.1  2.04  0.178
## 3 65.99        -0.6649                 4 -211.918 431.8  5.76  0.028
## 1 60.00                                3 -213.146 432.3  6.21  0.022
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## mean ~ pH.normalized + temp.average + (1 | Tank) + pH.normalized:temp.average
##    Data: juv.arctica
## REML criterion at convergence: 414.0776
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 1.583   
##  Residual             5.976   
## Number of obs: 66, groups:  Tank, 16
## Fixed Effects:
##                (Intercept)               pH.normalized  
##                    67.2318                     -2.2687  
##               temp.average  pH.normalized:temp.average  
##                    -0.5953                     -0.3574

simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

#plot(top_model1)
tab_model(top_model1)

	mean
Predictors	Estimates	CI	p
(Intercept)	67.23	53.09 – 81.37	<0.001
pH normalized	-2.27	-36.79 – 32.25	0.896
temp average	-0.60	-2.14 – 0.95	0.443
pH normalized × temp average	-0.36	-4.23 – 3.52	0.854
Random Effects
σ²	35.71
τ₀₀ _Tank	2.51
ICC	0.07
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.077 / 0.138

testZeroInflation(top_model1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(top_model1, juv.arctica$pH)

plotResiduals(top_model1, juv.arctica$temp)

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  #pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp & pH on Juvenile Arctica Coloration")

	Effect of Temp & pH on Juvenile Arctica Coloration
Predictors	Estimates	CI	p
(Intercept)	67.23	53.09 – 81.37	<0.001
pH normalized	-2.27	-36.79 – 32.25	0.896
temp average	-0.60	-2.14 – 0.95	0.443
pH normalized × temp average	-0.36	-4.23 – 3.52	0.854
Random Effects
σ²	35.71
τ₀₀ _Tank	2.51
ICC	0.07
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.077 / 0.138

formal plot

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

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

effect_top_model<-effects::effect(term= c("pH.normalized"), mod=top_model1)

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

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

##   pH.normalized      fit        se    lower    upper
## 1           0.0 61.86720 1.5220857 58.82459 64.90980
## 2           0.2 60.76937 0.9883871 58.79361 62.74512
## 3           0.3 60.22045 0.8528985 58.51553 61.92537
## 4           0.5 59.12262 1.0341075 57.05547 61.18977
## 5           0.7 58.02479 1.5913851 54.84366 61.20592

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

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

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

effect_top_model.both<-effects::effect(term= c("temp.average", "pH.normalized", "pH.normalized:temp.average"), mod=top_model1, at=values_of_temp)
effect_top_model.both<-as.data.frame(effect_top_model.both)
effect_top_model.both

##    temp.average pH.normalized      fit        se    lower    upper
## 1           6.1           0.0 63.60019 2.6706754 58.26159 68.93880
## 2           7.6           0.0 62.70718 1.8413610 59.02635 66.38800
## 3           9.1           0.0 61.81416 1.5259484 58.76383 64.86449
## 4          11.0           0.0 60.68300 2.1977740 56.28971 65.07629
## 5          12.0           0.0 60.08766 2.8060427 54.47846 65.69686
## 6           6.1           0.2 62.71043 1.7872351 59.13780 66.28306
## 7           7.6           0.2 61.71019 1.2336774 59.24411 64.17628
## 8           9.1           0.2 60.70996 0.9878956 58.73518 62.68473
## 9          11.0           0.2 59.44299 1.3787933 56.68682 62.19916
## 10         12.0           0.2 58.77617 1.7645888 55.24881 62.30353
## 11          6.1           0.3 62.26555 1.5486605 59.16982 65.36128
## 12          7.6           0.3 61.21170 1.0653220 59.08216 63.34125
## 13          9.1           0.3 60.15786 0.8527110 58.45331 61.86240
## 14         11.0           0.3 58.82299 1.1995414 56.42514 61.22084
## 15         12.0           0.3 58.12042 1.5378352 55.04634 61.19451
## 16          6.1           0.5 61.37579 1.7757275 57.82616 64.92542
## 17          7.6           0.5 60.21472 1.2063853 57.80319 62.62625
## 18          9.1           0.5 59.05366 1.0413501 56.97203 61.13529
## 19         11.0           0.5 57.58298 1.5872473 54.41012 60.75584
## 20         12.0           0.5 56.80894 2.0321646 52.74670 60.87118
## 21          6.1           0.7 60.48603 2.6527039 55.18335 65.78870
## 22          7.6           0.7 59.21774 1.7986725 55.62225 62.81324
## 23          9.1           0.7 57.94946 1.6067422 54.73763 61.16129
## 24         11.0           0.7 56.34297 2.5045700 51.33641 61.34953
## 25         12.0           0.7 55.49745 3.2006802 49.09938 61.89552

effect_temp <- ggplot() + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data = effect_top_model.both, aes(x = pH.normalized, y = fit, group = temp.average, color = temp.average), alpha = 0.7) +
  geom_ribbon(data= effect_top_model.both, aes(x=pH.normalized, ymin=lower, ymax=upper, group = temp.average, fill=temp.average), alpha= 0.3) +
  geom_point(data=mercenaria, mapping = aes(x=pH.normalized, y=mean, color=temp.average), size=3) +
  labs(x="Average pH (Normalized)", y="Mean Shell Coloration")+
  theme_classic()+
  scale_color_viridis(direction = -1)+
  scale_fill_viridis(direction = -1, guide=FALSE)
effect_temp

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

poisson

Warnings looks horrible, def not right.

## Fixed term is "(Intercept)"

## Global model call: glmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = juv.arctica, family = poisson, na.action = "na.fail")
## ---
## Model selection table 
##   (Int)    pH.nrm   tmp.avr pH.nrm:tmp.avr df logLik AIC
## 1 4.094                                     2   -Inf Inf
## 2 4.122 -0.085950                           3   -Inf Inf
## 3 4.191           -0.010690                 3   -Inf Inf
## 4 4.232 -0.095830 -0.011740                 4   -Inf Inf
## 8 4.202 -0.004978 -0.008393       -0.01036  5   -Inf Inf
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: mean ~ (1 | Tank)
##    Data: juv.arctica
##      AIC      BIC   logLik deviance df.resid 
##      Inf      Inf     -Inf      Inf       64 
## Random effects:
##  Groups Name        Std.Dev.
##  Tank   (Intercept) 1       
## Number of obs: 66, groups:  Tank, 16
## Fixed Effects:
## (Intercept)  
##       4.094  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 3498 optimizer warnings; 1 lme4 warnings

	mean
Predictors	Incidence Rate Ratios	CI	p
(Intercept)	59.99	36.71 – 98.02	<0.001
Random Effects
σ²	0.02
τ₀₀ _Tank	1.00
ICC	0.98
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.000 / 0.984

gamma distribution?

gamma requires all positive values.

This looks pretty bad, note that its likely not right because these are quite normal dist.

color.1 <- glmer(mean~temp.average*pH.normalized+(1|Tank), data = juv.arctica, family = Gamma(link = "log"),  na.action = "na.fail")

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00424067 (tol = 0.002, component 1)

simulationOutput1<-simulateResiduals(color.1)
plot(simulationOutput1) #wacky red lines

dd1<-dredge(color.1, rank=AIC)

## Fixed term is "(Intercept)"

## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00424065 (tol = 0.002, component 1)

dd1

## Global model call: glmer(formula = mean ~ temp.average * pH.normalized + (1 | Tank), 
##     data = juv.arctica, family = Gamma(link = "log"), na.action = "na.fail")
## ---
## Model selection table 
##   (Int)   pH.nrm   tmp.avr pH.nrm:tmp.avr df   logLik   AIC delta weight
## 1 4.093                                    3 -212.614 431.2  0.00  0.350
## 3 4.194          -0.011130                 4 -212.022 432.0  0.82  0.233
## 2 4.121 -0.08359                           4 -212.189 432.4  1.15  0.197
## 4 4.234 -0.09166 -0.012130                 5 -211.401 432.8  1.57  0.160
## 8 4.211 -0.02396 -0.009608      -0.007745  6 -211.383 434.8  3.54  0.060
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | Tank

top_model1<-get.models(dd1, subset = 1)[[1]]
top_model1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean ~ (1 | Tank)
##    Data: juv.arctica
##       AIC       BIC    logLik  deviance  df.resid 
##  431.2278  437.7968 -212.6139  425.2278        63 
## Random effects:
##  Groups   Name        Std.Dev.
##  Tank     (Intercept) 0.04364 
##  Residual             0.09657 
## Number of obs: 66, groups:  Tank, 16
## Fixed Effects:
## (Intercept)  
##       4.093

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: mean ~ (1 | Tank)
##    Data: juv.arctica
## 
##      AIC      BIC   logLik deviance df.resid 
##    431.2    437.8   -212.6    425.2       63 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.76546 -0.61555 -0.01767  0.49247  2.26119 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Tank     (Intercept) 0.001905 0.04364 
##  Residual             0.009326 0.09657 
## Number of obs: 66, groups:  Tank, 16
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept)  4.09350    0.02107   194.3   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

simulationOutput.top<-simulateResiduals(top_model1)
plot(simulationOutput.top)

testDispersion(simulationOutput.top)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.00068268, p-value < 2.2e-16
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	mean
Predictors	Estimates	CI	p
(Intercept)	59.95	57.48 – 62.53	<0.001
Random Effects
σ²	0.01
τ₀₀ _Tank	0.00
ICC	0.17
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.000 / 0.170

testZeroInflation(top_model1)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided

plotResiduals(top_model1, juv.arctica$ph)

plotResiduals(top_model1, juv.arctica$temp)

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH (normalized)", "Temp"),
                  dv.labels= "Effect of Temp and pH on juv.arctica Growth")

## Length of `pred.labels` does not equal number of predictors, no labelling applied.

	Effect of Temp and pH on juv.arctica Growth
Predictors	Estimates	CI	p
(Intercept)	59.95	57.48 – 62.53	<0.001
Random Effects
σ²	0.01
τ₀₀ _Tank	0.00
ICC	0.17
N _Tank	16
Observations	66
Marginal R² / Conditional R²	0.000 / 0.170

Plot

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

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

A-Maine-Tank-2022_All_Color_V2

McMahon

2023-07-31

Set up

read in relevant files

remove dead and lost

look at data

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

formal plot

poisson

gamma distribution?

Plot

linear model

‘normal’ linear mixed effect model

formal plot

poisson

gamma distribution?

Plot

scallop

check structure

Plot

juv arctica

check structure

linear model

‘normal’ linear mixed effect model

formal plot

poisson

gamma distribution?

Plot