Set up

read in relevant files

I took the dfs Diana sent me and read them in here.

remove dead and lost

We only had 59 removed. This seems a bit high to me and I wonder if the scallops that were removed weren’t marked in coloration? or if something else strange happened. But rolling with this for now.

## [1] 462  15

## [1] 403  15

## # A tibble: 4 × 2
##   species         n
##   <chr>       <int>
## 1 juv arctica    66
## 2 mercenaria    122
## 3 mya           150
## 4 scallop        65

look at data

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

## look for trends Scallop doesn’t have all temp pH groups covered, will not be able to model I don’t think… Also mercenaria 12 temp seems to almost be showing two trends, I wonder if this is explained by tank?

## Loading required package: ggpp

## 
## Attaching package: 'ggpp'

## The following object is masked from 'package:ggplot2':
## 
##     annotate

## Loading required package: viridisLite

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

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

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

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

Not totally sure if this is necessary for temperature too - did not do it here.

## [1] 7.39

# Individual Species Groups

mercenaria

check structure

Non normal, Binomial will not work here because values less than 0 - could consider removing those, but there are 34…?

## Classes 'data.table' and 'data.frame':   122 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ Sample_ID_20220415      : chr  "b38" "b39" "b40" "b41" ...
##  $ Max_heigh_mm            : num  12.4 11.4 10.2 10.8 11.5 ...
##  $ Max_height_mm           : num  18.3 17.2 15 15.1 16.3 ...
##  $ percentage change height: num  47.5 49.9 46 39.3 41.7 ...
##  $ 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 ...
##  $ died                    : 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 ...
##  $ prop.change.height      : num  0.475 0.499 0.46 0.393 0.417 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  mercenaria$prop.change.height
## W = 0.88872, p-value = 4.442e-08

##    vars   n mean  sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 122 0.17 0.2   0.11    0.14 0.18 -0.1 0.71  0.81 0.93     -0.1 0.02

## [1] 34

## [1] 88 17

linear model

Binomial will not work here because values less than 0, unless we filter out the 34 w less than 0. Get warnings for singular fit. The residuals with binomial look extremely wacky, I assume because so many points begin by being @/near 0?

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

## Fixed term is "(Intercept)"

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

## Global model call: glmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = mercenaria.small, family = binomial(link = "logit"), 
##     na.action = "na.fail")
## ---
## Model selection table 
##     (Int)   pH.nrm tmp.avr pH.nrm:tmp.avr df  logLik  AIC delta weight
## 4 -28.450    6.256 2.12500                 4 -13.952 35.9  0.00  0.471
## 8  -4.888 -508.900 0.04662          43.43  5 -13.000 36.0  0.09  0.449
## 3 -24.410          1.97400                 3 -16.724 39.4  3.54  0.080
## 2  -3.442    3.087                         3 -25.249 56.5 20.59  0.000
## 1  -2.303                                  2 -26.808 57.6 21.71  0.000
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | tank

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

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: prop.change.height ~ pH.normalized + temp.average + (1 | tank)
##    Data: mercenaria.small
##      AIC      BIC   logLik deviance df.resid 
##  35.9048  45.8141 -13.9524  27.9048       84 
## Random effects:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0       
## Number of obs: 88, groups:  tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       -28.453          6.256          2.125  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

## qu = 0.25, log(sigma) = -6.108266 : outer Newton did not converge fully.

## qu = 0.25, log(sigma) = -6.240541 : outer Newton did not converge fully.

	prop.change.height
Predictors	Odds Ratios	CI	p
(Intercept)	0.00	0.00 – 1.48	0.053
pH normalized	521.31	1.32 – 205594.76	0.040
temp average	8.37	0.78 – 90.29	0.080
Random Effects
σ²	3.29
τ₀₀ _tank	0.00
N _tank	16
Observations	88
Marginal R² / Conditional R²	0.828 / NA

beta distribution?

this appears to look much better but I wish I understood the difference between glmmTMB and glmer more? and how beta differs from binomial in glmer. I think this would work if we wanted, but also should consider the fact that 34 were removed… so this is only the mercenaria that grew?

#install.packages("glmmTMB")
library(glmmTMB)

height.1 <- glmmTMB(prop.change.height~temp.average*pH.normalized+(1|tank), data = mercenaria.small, family = beta_family(link = "logit"),  na.action = "na.fail")

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

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

## Fixed terms are "cond((Int))" and "disp((Int))"

dd1

## Global model call: glmmTMB(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = mercenaria.small, family = beta_family(link = "logit"), 
##     na.action = "na.fail", ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##   cnd((Int)) dsp((Int)) cnd(pH.nrm) cnd(tmp.avr) cnd(pH.nrm:tmp.avr) df logLik
## 4     -7.390          +      0.8576       0.5921                      5 98.647
## 8     -7.387          +      0.8486       0.5918           0.0009293  6 98.647
## 3     -7.135          +                   0.5930                      4 96.599
## 1     -1.717          +                                               3 78.450
## 2     -1.825          +      0.3362                                   4 78.482
##      AIC delta weight
## 4 -187.3  0.00  0.582
## 8 -185.3  2.00  0.214
## 3 -185.2  2.09  0.204
## 1 -150.9 36.39  0.000
## 2 -149.0 38.33  0.000
## Models ranked by AIC(x) 
## Random terms (all models): 
##   cond(1 | tank)

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

## Formula:          
## prop.change.height ~ pH.normalized + temp.average + (1 | tank)
## Data: mercenaria.small
##        AIC        BIC     logLik   df.resid 
## -187.29364 -174.90695   98.64682         83 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0.2012  
## 
## Number of obs: 88 / Conditional model: tank, 16
## 
## Dispersion parameter for beta family ():   16 
## 
## Fixed Effects:
## 
## Conditional model:
##   (Intercept)  pH.normalized   temp.average  
##       -7.3897         0.8576         0.5921

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

## qu = 0.75, log(sigma) = -2.681231 : outer Newton did not converge fully.

testDispersion(simulationOutput.top)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.78183, p-value = 0.192
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.00	0.00 – 0.00	<0.001
pH normalized	2.36	1.09 – 5.08	0.029
temp average	1.81	1.63 – 2.01	<0.001
Random Effects
σ²	0.24
τ₀₀ _tank	0.04
ICC	0.15
N _tank	16
Observations	88
Marginal R² / Conditional R²	0.799 / 0.828

plot_model(top_model1, 
                   axis.labels=c("pH", "Temp"),
                   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)	0.00	0.00 – 0.00	<0.001
pH (normalized)	2.36	1.09 – 5.08	0.029
Temp	1.81	1.63 – 2.01	<0.001
Random Effects
σ²	0.24
τ₀₀ _tank	0.04
ICC	0.15
N _tank	16
Observations	88
Marginal R² / Conditional R²	0.799 / 0.828

formal plot

effect_top_model.temp<-effects::effect(term="temp.average", mod=top_model1)
effect_top_model.temp<-as.data.frame(effect_top_model.temp)

effect_temp <- ggplot() + 
  geom_point(data=mercenaria.small, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #4
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  #5
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  labs(x="Average Temperature (C)", y="Proportional Change in Height")+
  theme_classic()+
  scale_color_viridis(direction = -1)
effect_temp

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

effect_ph <- ggplot() + 
  geom_point(data=mercenaria.small, mapping = aes(x=pH.normalized, y=prop.change.height, color=temp.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=pH.normalized, y=fit), color="black") +
  #4
  geom_line(data=effect_top_model.ph, aes(x=pH.normalized, y=fit), color="black") +
  #5
  geom_ribbon(data= effect_top_model.ph, aes(x=pH.normalized, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  labs(x="Average pH Normalized", y="Proportional Change in Height")+
  theme_classic()+
  scale_color_viridis(direction = 1)
effect_ph

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

mya

check structure

This distribution looks pretty normal, although it still looks like there are values below zero (2). Again, maybe Brittany has advice on how this could be addressed but I will likely just remove them here.

## Classes 'data.table' and 'data.frame':   150 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 2 ...
##  $ Sample_ID_20220415      : chr  "r37" "r38" "r39" "r40" ...
##  $ Max_heigh_mm            : num  17.1 13.8 13.8 16.5 15.7 ...
##  $ Max_height_mm           : num  22.2 15.6 19.9 21.6 21.4 ...
##  $ percentage change height: num  30.2 12.7 44.6 31 36.1 ...
##  $ 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 ...
##  $ died                    : 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 ...
##  $ prop.change.height      : num  0.302 0.127 0.446 0.31 0.361 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  mya$prop.change.height
## W = 0.99276, p-value = 0.6518

##    vars   n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 150 0.25 0.12   0.26    0.25 0.14 -0.1 0.54  0.64 -0.09    -0.44 0.01

## [1] 2

## [1] 148  17

linear model

A linear model with a binomial family, still getting singular fit warnings and top model is only including random effect of tank, will try beta as before. All residuals are 1 bc there are no fixed effects?

## Fixed term is "(Intercept)"

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

## Global model call: glmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = mya.small, family = binomial(link = "logit"), 
##     na.action = "na.fail")
## ---
## Model selection table 
##    (Int)   pH.nrm  tmp.avr pH.nrm:tmp.avr df  logLik  AIC delta weight
## 1 -3.584                                   2 -18.389 40.8  0.00  0.456
## 3 -5.742           0.23080                 3 -18.006 42.0  1.23  0.246
## 2 -3.521 -0.18960                          3 -18.386 42.8  1.99  0.168
## 4 -5.754  0.02744  0.23110                 4 -18.006 44.0  3.23  0.090
## 8 -3.474 -7.08200 -0.01399         0.7601  5 -17.817 45.6  4.85  0.040
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | tank

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

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: prop.change.height ~ (1 | tank)
##    Data: mya.small
##      AIC      BIC   logLik deviance df.resid 
##  40.7782  46.7727 -18.3891  36.7782      146 
## Random effects:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0       
## Number of obs: 148, groups:  tank, 16
## Fixed Effects:
## (Intercept)  
##      -3.584  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

	prop.change.height
Predictors	Odds Ratios	CI	p
(Intercept)	0.03	0.01 – 0.08	<0.001
Random Effects
σ²	3.29
τ₀₀ _tank	0.00
N _tank	16
Observations	148
Marginal R² / Conditional R²	0.000 / NA

beta distribution?

I wonder why beta seems to work best? I wonder if I need to interpret this difference? Qunatile deviations detected in top model residuals v predicted, lines are looking like parabolas, not sure what that means but I can see whether this is an issue, dispersion seems okay, QQ line looks good.

#install.packages("glmmTMB")
library(glmmTMB)

height.1 <- glmmTMB(prop.change.height~temp.average*pH.normalized+(1|tank), data = mya.small, family = beta_family(link = "logit"),  na.action = "na.fail")

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

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

## Fixed terms are "cond((Int))" and "disp((Int))"

dd1

## Global model call: glmmTMB(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = mya.small, family = beta_family(link = "logit"), 
##     na.action = "na.fail", ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##   cnd((Int)) dsp((Int)) cnd(pH.nrm) cnd(tmp.avr) cnd(pH.nrm:tmp.avr) df  logLik
## 3     -2.064          +                   0.1092                      4 109.089
## 4     -2.072          +     0.01883       0.1095                      5 109.091
## 8     -2.002          +    -0.19030       0.1016             0.02397  6 109.106
## 1     -1.082          +                                               3 104.340
## 2     -1.053          +    -0.08581                                   4 104.369
##      AIC delta weight
## 3 -210.2  0.00  0.650
## 4 -208.2  2.00  0.240
## 8 -206.2  3.96  0.090
## 1 -202.7  7.50  0.015
## 2 -200.7  9.44  0.006
## Models ranked by AIC(x) 
## Random terms (all models): 
##   cond(1 | tank)

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

## Formula:          prop.change.height ~ temp.average + (1 | tank)
## Data: mya.small
##       AIC       BIC    logLik  df.resid 
## -210.1772 -198.1884  109.0886       144 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0.1266  
## 
## Number of obs: 148 / Conditional model: tank, 16
## 
## Dispersion parameter for beta family ():   12 
## 
## Fixed Effects:
## 
## Conditional model:
##  (Intercept)  temp.average  
##      -2.0636        0.1092

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.84488, p-value = 0.192
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.13	0.07 – 0.22	<0.001
temp average	1.12	1.05 – 1.18	<0.001
Random Effects
σ²	0.14
τ₀₀ _tank	0.02
ICC	0.10
N _tank	16
Observations	148
Marginal R² / Conditional R²	0.225 / 0.305

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

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

	Effect of Temp on Mya Growth
Predictors	Estimates	CI	p
(Intercept)	0.13	0.07 – 0.22	<0.001
Temp	1.12	1.05 – 1.18	<0.001
Random Effects
σ²	0.14
τ₀₀ _tank	0.02
ICC	0.10
N _tank	16
Observations	148
Marginal R² / Conditional R²	0.225 / 0.305

formal plot

I would like to understand beta distribution more, this looks fully linear for the predicted line but the last one seemed to be slighltly linear? Not sure why this is.

effect_top_model.temp<-effects::effect(term="temp.average", mod=top_model1)
effect_top_model.temp<-as.data.frame(effect_top_model.temp)

effect_temp <- ggplot() + 
  geom_point(data=mya.small, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #4
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  #5
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  labs(x="Average Temperature (C)", y="Proportional Change in Height")+
  theme_classic()+
  scale_color_viridis(direction = -1)
effect_temp

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

#effect_ph <- ggplot() + 
  #geom_point(data=mya.small, mapping = aes(x=pH.normalized, y=prop.change.height, #color=temp.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=pH.normalized, y=fit), color="black") +
  #4
  #geom_line(data=effect_top_model.ph, aes(x=pH.normalized, y=fit), color="black") +
  #5
  #geom_ribbon(data= effect_top_model.ph, aes(x=pH.normalized, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  #labs(x="Average pH Normalized", y="Proportional Change in Height")+
  #theme_classic()+
  #scale_color_viridis(direction = 1)
#effect_ph

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

juvenile arctica

check structure

This distribution looks pretty normal, although it looks like there is one outlier above 1. Again, maybe Brittany has advice on how this could be addressed but I will likely just remove them here. G33, G44 and G8 are na? Don’t have max height at beginning? will remove

## Classes 'data.table' and 'data.frame':   66 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 2 2 2 2 3 3 ...
##  $ Sample_ID_20220415      : chr  "g29" "g30" "g31" "g32" ...
##  $ Max_heigh_mm            : num  13.3 14.2 16.2 15 13.1 ...
##  $ Max_height_mm           : num  25.5 23.9 28.7 27.7 22.2 ...
##  $ percentage change height: num  92.3 67.6 77.3 84.4 70.3 ...
##  $ 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 ...
##  $ died                    : 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 ...
##  $ prop.change.height      : num  0.923 0.676 0.773 0.844 0.703 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.61 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  juv.arctica$prop.change.height
## W = 0.96809, p-value = 0.1012

##    vars  n mean  sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 63 0.68 0.2   0.65    0.67 0.21 0.31 1.28  0.97 0.52    -0.25 0.03

## [1] 0

## [1] 4

## [1] 58 17

linear model

best fit includes just temp but resid a bit strange.

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

## Fixed term is "(Intercept)"

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

## Global model call: glmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = juv.arctica.small, family = binomial(link = "logit"), 
##     na.action = "na.fail")
## ---
## Model selection table 
##     (Int)   pH.nrm tmp.avr pH.nrm:tmp.avr df  logLik  AIC delta weight
## 3 -2.7820          0.51120                 3 -24.879 55.8  0.00  0.528
## 4 -2.6370  -0.4190 0.51240                 4 -24.842 57.7  1.93  0.202
## 8  0.7135 -10.9300 0.07619          1.382  5 -23.872 57.7  1.99  0.196
## 1  1.4520                                  2 -28.172 60.3  4.59  0.053
## 2  1.6240  -0.4699                         3 -28.124 62.2  6.49  0.021
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | tank

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

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: juv.arctica.small
##      AIC      BIC   logLik deviance df.resid 
##  55.7584  61.9397 -24.8792  49.7584       55 
## Random effects:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0       
## Number of obs: 58, groups:  tank, 16
## Fixed Effects:
##  (Intercept)  temp.average  
##      -2.7822        0.5112  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

## qu = 0.25, log(sigma) = -4.914393 : outer Newton did not converge fully.

	prop.change.height
Predictors	Odds Ratios	CI	p
(Intercept)	0.06	0.00 – 1.95	0.114
temp average	1.67	1.08 – 2.57	0.020
Random Effects
σ²	3.29
τ₀₀ _tank	0.00
N _tank	16
Observations	58
Marginal R² / Conditional R²	0.228 / NA

beta distribution?

beta includes all in top model…., but only ph and int are significant? no idea how to plot interaction term.

#install.packages("glmmTMB")
library(glmmTMB)

height.1 <- glmmTMB(prop.change.height~temp.average*pH.normalized+(1|tank), data = juv.arctica.small, family = beta_family(link = "logit"),  na.action = "na.fail")

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

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

## Fixed terms are "cond((Int))" and "disp((Int))"

dd1

## Global model call: glmmTMB(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = juv.arctica.small, family = beta_family(link = "logit"), 
##     na.action = "na.fail", ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##   cnd((Int)) dsp((Int)) cnd(pH.nrm) cnd(tmp.avr) cnd(pH.nrm:tmp.avr) df logLik
## 8    -0.3788          +     -3.8720       0.1177              0.4501  6 36.981
## 3    -1.7290          +                   0.2740                      4 34.669
## 4    -1.6710          +     -0.1767       0.2748                      5 34.782
## 1     0.7244          +                                               3 25.363
## 2     0.7792          +     -0.1608                                   4 25.386
##     AIC delta weight
## 8 -62.0  0.00  0.492
## 3 -61.3  0.62  0.360
## 4 -59.6  2.40  0.148
## 1 -44.7 17.24  0.000
## 2 -42.8 19.19  0.000
## Models ranked by AIC(x) 
## Random terms (all models): 
##   cond(1 | tank)

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

## Formula:          
## prop.change.height ~ pH.normalized + temp.average + (1 | tank) +  
##     pH.normalized:temp.average
## Data: juv.arctica.small
##       AIC       BIC    logLik  df.resid 
## -61.96188 -49.59922  36.98094        52 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups Name        Std.Dev. 
##  tank   (Intercept) 1.718e-05
## 
## Number of obs: 58 / Conditional model: tank, 16
## 
## Dispersion parameter for beta family (): 10.6 
## 
## Fixed Effects:
## 
## Conditional model:
##                (Intercept)               pH.normalized  
##                    -0.3788                     -3.8715  
##               temp.average  pH.normalized:temp.average  
##                     0.1177                      0.4501

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.89521, p-value = 0.568
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.68	0.17 – 2.82	0.600
pH normalized	0.02	0.00 – 0.67	0.029
temp average	1.12	0.95 – 1.33	0.162
pH normalized × temp average	1.57	1.04 – 2.37	0.033
Random Effects
σ²	-0.05
τ₀₀ _tank	0.00
N _tank	16
Observations	58
Marginal R² / Conditional R²	1.173 / NA

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

tab_model(top_model1, 
                  show.re.var= TRUE, 
                  pred.labels =c("(Intercept)", "pH Normalized", "Temp:pH Interaction"),
                  dv.labels= "Effect of Temp on Juvenile Arctica Growth")

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

	Effect of Temp on Juvenile Arctica Growth
Predictors	Estimates	CI	p
(Intercept)	0.68	0.17 – 2.82	0.600
pH.normalized	0.02	0.00 – 0.67	0.029
temp.average	1.12	0.95 – 1.33	0.162
pH.normalized:temp.average	1.57	1.04 – 2.37	0.033
Random Effects
σ²	-0.05
τ₀₀ _tank	0.00
N _tank	16
Observations	58
Marginal R² / Conditional R²	1.173 / NA

formal plot

Idk how to plot interaction. CI for pH includes 0, how should temp be interpreted if p>0.05

effect_top_model.temp<-effects::effect(term="temp.average", mod=top_model1)

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

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

effect_temp <- ggplot() + 
  geom_point(data=juv.arctica.small, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #4
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  #5
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  labs(x="Average Temperature (C)", y="Proportional Change in Height")+
  theme_classic()+
  scale_color_viridis(direction = -1)
effect_temp

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

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

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

effect_ph <- ggplot() + 
  geom_point(data=juv.arctica.small, mapping = aes(x=pH.normalized, y=prop.change.height, color=temp.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=pH.normalized, y=fit), color="black") +
  #4
  geom_line(data=effect_top_model.ph, aes(x=pH.normalized, y=fit), color="black") +
  #5
  geom_ribbon(data= effect_top_model.ph, aes(x=pH.normalized, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  labs(x="Average pH Normalized", y="Proportional Change in Height")+
  theme_classic()+
  scale_color_viridis(direction = 1)
effect_ph

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

scallop

check structure

This distribution looks pretty normal, although it looks quite kurtose. At least one above 1 that I will need to remove.

## Classes 'data.table' and 'data.frame':   65 obs. of  17 variables:
##  $ tank                    : Factor w/ 16 levels "H1","H10","H11",..: 1 1 2 2 2 3 3 3 3 3 ...
##  $ Sample_ID_20220415      : chr  "r26" "r18" "r34" "r55" ...
##  $ Max_heigh_mm            : num  29.7 35.5 28.1 39.2 37.7 ...
##  $ Max_height_mm           : num  NA NA 57.4 63.2 63.4 ...
##  $ percentage change height: num  NA NA 103.8 61.2 68.2 ...
##  $ ph                      : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 4 4 4 4 4 4 4 4 ...
##  $ temp                    : Factor w/ 3 levels "6","9","12": 3 3 2 2 2 1 1 1 1 1 ...
##  $ died                    : chr  NA NA NA NA ...
##  $ species                 : chr  "scallop" "scallop" "scallop" "scallop" ...
##  $ temp.average            : num  12.06 12.06 9.19 9.19 9.19 ...
##  $ temp.stdev              : num  0.42 0.42 0.63 0.63 0.63 0.78 0.78 0.78 0.78 0.78 ...
##  $ pH.average.YSI          : num  7.57 7.57 8 8 8 8.02 8.02 8.02 8.02 8.02 ...
##  $ pH.stdev.YSI            : num  0.04 0.04 0.07 0.07 0.07 0.06 0.06 0.06 0.06 0.06 ...
##  $ pH.average              : num  7.59 7.59 8 8 8 8.05 8.05 8.05 8.05 8.05 ...
##  $ pH.stdev                : num  0.07 0.07 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 ...
##  $ prop.change.height      : num  NA NA 1.038 0.612 0.682 ...
##  $ pH.normalized           : num  0.2 0.2 0.61 0.61 0.61 ...
##  - attr(*, ".internal.selfref")=<externalptr>

## 
##  Shapiro-Wilk normality test
## 
## data:  scallop$prop.change.height
## W = 0.9722, p-value = 0.5434

##    vars  n mean   sd median trimmed  mad  min  max range  skew kurtosis   se
## X1    1 33 0.68 0.18   0.69    0.68 0.17 0.28 1.04  0.76 -0.36    -0.35 0.03

## [1] 0

## [1] 1

## [1] 32 17

linear model

Not working “Error in pwrssUpdate(pp, resp, tol = tolPwrss, GQmat = GQmat, compDev = compDev, : Downdated VtV is not positive definite” ?? seems like it could be do to 0 or 1 inflation? maybe its due to low sample size here, will comment out to builder html.

beta distribution?

beta includes temp in top model. no apparent issues with residuals etc

#install.packages("glmmTMB")
library(glmmTMB)

height.1 <- glmmTMB(prop.change.height~temp.average*pH.normalized+(1|tank), data = scallop.small, family = beta_family(link = "logit"),  na.action = "na.fail")

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

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

## Fixed terms are "cond((Int))" and "disp((Int))"

dd1

## Global model call: glmmTMB(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = scallop.small, family = beta_family(link = "logit"), 
##     na.action = "na.fail", ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##   cnd((Int)) dsp((Int)) cnd(pH.nrm) cnd(tmp.avr) cnd(pH.nrm:tmp.avr) df logLik
## 3    -1.8490          +                   0.3096                      4 25.914
## 4    -1.8330          +    -0.02868       0.3088                      5 25.916
## 8    -1.8490          +     0.02640       0.3107           -0.006687  6 25.916
## 1     0.8328          +                                               3 19.220
## 2     1.0010          +    -0.52060                                   4 19.429
##     AIC delta weight
## 3 -43.8  0.00  0.663
## 4 -41.8  2.00  0.244
## 8 -39.8  4.00  0.090
## 1 -32.4 11.39  0.002
## 2 -30.9 12.97  0.001
## Models ranked by AIC(x) 
## Random terms (all models): 
##   cond(1 | tank)

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

## Formula:          prop.change.height ~ temp.average + (1 | tank)
## Data: scallop.small
##       AIC       BIC    logLik  df.resid 
## -43.82837 -37.96542  25.91418        28 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0.2807  
## 
## Number of obs: 32 / Conditional model: tank, 14
## 
## Dispersion parameter for beta family (): 20.5 
## 
## Fixed Effects:
## 
## Conditional model:
##  (Intercept)  temp.average  
##      -1.8487        0.3096

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.94911, p-value = 0.864
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.16	0.05 – 0.47	0.001
temp average	1.36	1.20 – 1.55	<0.001
Random Effects
σ²	-0.03
τ₀₀ _tank	0.08
ICC	1.51
N _tank	14
Observations	32
Marginal R² / Conditional R²	0.873 / 1.065

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

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

	Effect of Temp on Scallop Growth
Predictors	Estimates	CI	p
(Intercept)	0.16	0.05 – 0.47	0.001
Temp	1.36	1.20 – 1.55	<0.001
Random Effects
σ²	-0.03
τ₀₀ _tank	0.08
ICC	1.51
N _tank	14
Observations	32
Marginal R² / Conditional R²	0.873 / 1.065

formal plot

Idk how to plot interaction. CI for pH includes 0, how should temp be interpreted if p>0.05

effect_top_model.temp<-effects::effect(term="temp.average", mod=top_model1)
effect_top_model.temp<-as.data.frame(effect_top_model.temp)

effect_temp <- ggplot() + 
  geom_point(data=scallop.small, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) + 
  #3
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  #4
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  #5
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  #6
  labs(x="Average Temperature (C)", y="Proportional Change in Height")+
  theme_classic()+
  scale_color_viridis(direction = -1)
effect_temp

scallop.model.plot<-ggarrange(effect_temp,
              labels = c(""),
              ncol = 2, nrow = 1)
ggsave("02_output/02_plots/height_scallop.png", scallop.model.plot, width = 10, height = 4, dpi = 300)

A-Maine-Tank-2022_Bivalve-Height_V2

McMahon

2023-06-19

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)

mercenaria

check structure

linear model

beta distribution?

formal plot

mya

check structure

linear model

beta distribution?

formal plot

juvenile arctica

check structure

linear model

beta distribution?

formal plot

scallop

check structure

linear model

beta distribution?

formal plot