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 66 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. From Diana there should be 66 deaths across these groups.

## [1] 462  15

## [1] 396  15

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

## [1] 35

look at data

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

strange specimens

based on our conversation with Brittany, it seems safe to assume that measurements where growth is below 0 are mostly a reflection of sampling error and not actual loss in shell max height over the experiment. There I would like to make the assumption that any - values here are really a reflection of 0 growth. In some sense we are flooring percent change in height to 0. At this point all specimens marked as dead or lost have been removed so this shouldn’t effect those and all other specimens with missing data are NA for change in height so they wont be included in this either.

Another approach to this would be to add the value of the most negative number to the whole column so that the lowest % change is 0 % change, but this makes less sense to me in the context of this data set.

height.all.w.negatives<-height.all
height.all$`percentage change height` <- replace(height.all$`percentage change height`, which(height.all$`percentage change height` < 0), 0)

height.all$prop.change.height<-height.all$`percentage change height`/100
height.all %>% 
  ggplot(aes(prop.change.height, group = species, fill=species)) +
    geom_boxplot()

## Warning: Removed 35 rows containing non-finite values (`stat_boxplot()`).

height.all %>% 
  ggplot(aes(x=temp.average, y=prop.change.height,color = species, fill=species)) +
    geom_point()+
    geom_smooth(method = loess) #sweet no longer any below 0

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

## Warning: Removed 35 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 35 rows containing missing values (`geom_point()`).

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?

In this version it is updated so each species shows geom_smooth with each temp/pH grouping its own line. pH scale is switched so light color = more acidic. Remove SE coloring as plot was too busy.Noting here that it could be nice to have this plot in Altair in the future, so you could hover over points, but thats a rabbit hole for another day.

## 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, majorly skewed many at 0

## 'data.frame':    115 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  "b47" "b39" "b38" "b43" ...
##  $ Max_heigh_mm            : num  11.4 11.4 12.4 11.9 10.2 ...
##  $ Max_height_mm           : num  17.4 17.2 18.3 17.2 15 ...
##  $ percentage change height: num  53.6 49.9 47.5 44.7 46 ...
##  $ 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.536 0.499 0.475 0.447 0.46 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...

## 
##  Shapiro-Wilk normality test
## 
## data:  mercenaria$prop.change.height
## W = 0.85068, p-value = 2.09e-09

##    vars   n mean  sd median trimmed  mad min  max range skew kurtosis   se
## X1    1 115 0.18 0.2   0.13    0.16 0.19   0 0.71  0.71 0.97    -0.08 0.02

## [1] 0

## [1] 115  17

##     ph temp  n
## 1  7.4    6  9
## 2  7.4    9 14
## 3  7.4   12  5
## 4  7.6    6  9
## 5  7.6    9 14
## 6  7.6   12  9
## 7  7.8    6  8
## 8  7.8    9 10
## 9  7.8   12  6
## 10   8    6  9
## 11   8    9 15
## 12   8   12  7

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 includes temp, observed vs expected is off QQ line, residual vs predicted is wacky. We can see if other models look better or worse.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = mercenaria, na.action = na.fail)

simulationOutput1<-simulateResiduals(height.1)

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

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

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = prop.change.height ~ 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    AIC delta weight
## 3 -0.5860          0.08597                 4 123.013 -238.0  0.00  0.909
## 4 -0.6048  0.04371 0.08646                 5 121.657 -233.3  4.71  0.086
## 8 -0.5182 -0.22050 0.07669         0.0302  6 119.730 -227.5 10.56  0.005
## 1  0.1918                                  3 113.432 -220.9 17.16  0.000
## 2  0.2023 -0.03208                         4 112.839 -217.7 20.35  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: prop.change.height ~ temp.average + (1 | tank)
##    Data: mercenaria
## REML criterion at convergence: -246.0258
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.07756 
##  Residual             0.06728 
## Number of obs: 115, groups:  tank, 16
## Fixed Effects:
##  (Intercept)  temp.average  
##     -0.58596       0.08597

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

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

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	-0.59	-0.77 – -0.40	<0.001
temp average	0.09	0.07 – 0.11	<0.001
Random Effects
σ²	0.00
τ₀₀ _tank	0.01
ICC	0.57
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.750 / 0.893

testZeroInflation(top_model1)

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

plotResiduals(top_model1, mercenaria$ph)

plotResiduals(top_model1, mercenaria$temp)

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

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

	Effect of Temp on Mercenaria Growth
Predictors	Estimates	CI	p
(Intercept)	-0.59	-0.77 – -0.40	<0.001
Temp	0.09	0.07 – 0.11	<0.001
Random Effects
σ²	0.00
τ₀₀ _tank	0.01
ICC	0.57
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.750 / 0.893

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more

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

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

binomial

Warnings for singular fit. The residuals with binomial look extremely wacky, I assume because so many points begin by being @/near 0? Looks way worse than above model.

## 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, 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.800    6.258   2.154                 4 -13.960 35.9  0.00  0.476
## 8  -6.396 -492.300   0.176          42.03  5 -13.032 36.1  0.15  0.443
## 3 -24.490            1.981                 3 -16.725 39.5  3.53  0.081
## 1  -2.593                                  2 -29.039 62.1 26.16  0.000
## 2  -3.504    2.376                         3 -28.055 62.1 26.19  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
##      AIC      BIC   logLik deviance df.resid 
##  35.9196  46.8993 -13.9598  27.9196      111 
## Random effects:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0       
## Number of obs: 115, groups:  tank, 16
## Fixed Effects:
##   (Intercept)  pH.normalized   temp.average  
##       -28.805          6.258          2.154  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

	prop.change.height
Predictors	Odds Ratios	CI	p
(Intercept)	0.00	0.00 – 1.02	0.050
pH normalized	522.21	1.31 – 207426.45	0.040
temp average	8.62	0.80 – 92.73	0.076
Random Effects
σ²	3.29
τ₀₀ _tank	0.00
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.864 / NA

gamma distribution?

gamma requires all positive values, I have 0s here, mainly from flooring done earlier. here I will add a verrrrry small constant (0.00000000000001) to all measurements.

This looks pretty bad too which oddly means linear was best fit… Here we can see it is zero inflated and has some weird dispersion. Non homogenous variance striking when breaking down by group.

mercenaria$prop.change.height.1<-mercenaria$prop.change.height+0.00000000000001
height.1 <- glmer(prop.change.height.1~temp.average*pH.normalized+(1|tank), data = mercenaria, family = Gamma,  na.action = "na.fail")

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

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

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

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

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

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

dd1

## Global model call: glmer(formula = prop.change.height.1 ~ temp.average * pH.normalized + 
##     (1 | tank), data = mercenaria, family = Gamma, na.action = "na.fail")
## ---
## Model selection table 
##    (Int)  pH.nrm tmp.avr pH.nrm:tmp.avr df  logLik     AIC delta weight
## 3 39.590          -3.160                 4 706.199 -1404.4  0.00  0.659
## 4 39.920 -0.8063  -3.166                 5 706.209 -1402.4  1.98  0.245
## 8 37.480  6.4840  -2.955        -0.6316  6 706.217 -1400.4  3.96  0.091
## 1  5.425                                 3 700.160 -1394.3 10.08  0.004
## 2  5.330  0.2915                         4 700.161 -1392.3 12.08  0.002
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | tank

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

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

top_model1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: prop.change.height.1 ~ temp.average + (1 | tank)
##    Data: mercenaria
##        AIC        BIC     logLik   deviance   df.resid 
## -1404.3977 -1393.4180   706.1989 -1412.3977        111 
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.0000  
##  Residual             0.9138  
## Number of obs: 115, groups:  tank, 16
## Fixed Effects:
##  (Intercept)  temp.average  
##        39.59         -3.16  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: prop.change.height.1 ~ temp.average + (1 | tank)
##    Data: mercenaria
## 
##      AIC      BIC   logLik deviance df.resid 
##  -1404.4  -1393.4    706.2  -1412.4      111 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.0943 -1.0780 -0.1626  0.6515  2.6162 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.000    0.0000  
##  Residual             0.835    0.9138  
## Number of obs: 115, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value Pr(>|z|)   
## (Intercept)    39.591     13.803   2.868  0.00413 **
## temp.average   -3.160      1.171  -2.700  0.00694 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.997
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

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.15531, p-value < 2.2e-16
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	prop.change.height.1
Predictors	Estimates	CI	p
(Intercept)	156325575134083200.00	206663.77 – 118248520559533590254375141376.00	0.005
temp average	0.04	0.00 – 0.43	0.008
Random Effects
σ²	0.83
τ₀₀ _tank	0.00
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.981 / NA

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

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

	Effect of Temp and pH on Mercenaria Growth
Predictors	Estimates	CI	p
(Intercept)	156325575134083200.00	206663.77 – 118248520559533590254375141376.00	0.005
temp.average	0.04	0.00 – 0.43	0.008
Random Effects
σ²	0.83
τ₀₀ _tank	0.00
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.981 / NA

mya

check structure

This distribution looks pretty normal.

## '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  "r39" "r40" "r37" "r38" ...
##  $ Max_heigh_mm            : num  13.8 16.5 17.1 13.8 14.3 ...
##  $ Max_height_mm           : num  19.9 21.6 22.2 15.6 15 ...
##  $ percentage change height: num  44.62 31.04 30.25 12.68 4.76 ...
##  $ 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.4462 0.3104 0.3025 0.1268 0.0476 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.61 ...

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

##    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 0.54  0.54 0.01    -0.66 0.01

##     ph temp  n
## 1  7.4    6 10
## 2  7.4    9 20
## 3  7.4   12  7
## 4  7.6    6 10
## 5  7.6    9 19
## 6  7.6   12  9
## 7  7.8    6  9
## 8  7.8    9 20
## 9  7.8   12  9
## 10   8    6 11
## 11   8    9 19
## 12   8   12  7

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 includes none, which seems to make sense based on visual inspection of the data. residuals etc look pretty good but I will try a couple other families. I assume this will be best.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = mya, na.action = na.fail)

simulationOutput1<-simulateResiduals(height.1)

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

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

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = mya, na.action = na.fail)
## ---
## Model selection table 
##     (Int)   pH.nrm tmp.avr pH.nrm:tmp.avr df  logLik    AIC delta weight
## 1 0.25580                                  3 104.546 -203.1  0.00  0.574
## 3 0.06027          0.02171                 4 105.076 -202.2  0.94  0.359
## 2 0.27070 -0.04520                         4 103.025 -198.0  5.04  0.046
## 4 0.07204 -0.02659 0.02137                 5 103.197 -196.4  6.70  0.020
## 8 0.06910 -0.01642 0.02171      -0.001191  6 100.568 -189.1 13.96  0.001
## 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: prop.change.height ~ (1 | tank)
##    Data: mya
## REML criterion at convergence: -209.0921
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.05255 
##  Residual             0.11148 
## Number of obs: 150, groups:  tank, 16
## Fixed Effects:
## (Intercept)  
##      0.2558

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

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.26	0.22 – 0.29	<0.001
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.18
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.000 / 0.182

testZeroInflation(top_model1)

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

plotResiduals(top_model1, mya$ph)

plotResiduals(top_model1, mya$temp)

#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)"),
                  dv.labels= "Effect of Temp on mya Growth")

	Effect of Temp on mya Growth
Predictors	Estimates	CI	p
(Intercept)	0.26	0.22 – 0.29	<0.001
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.18
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.000 / 0.182

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more

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

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

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

binomial

Warnings for singular fit. The residuals with binomial look extremely wacky, I assume because so many points begin by being @/near 0? Looks way worse than above model.

## 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 = mya, 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.597                                   2 -18.444 40.9  0.00  0.453
## 3 -5.795          0.235200                 3 -18.046 42.1  1.20  0.248
## 2 -3.545 -0.16080                          3 -18.441 42.9  1.99  0.167
## 4 -5.810  0.03572 0.235600                 4 -18.045 44.1  3.20  0.091
## 8 -3.655 -6.76300 0.003334         0.7277  5 -17.869 45.7  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
##      AIC      BIC   logLik deviance df.resid 
##  40.8871  46.9084 -18.4435  36.8871      148 
## Random effects:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0       
## Number of obs: 150, groups:  tank, 16
## Fixed Effects:
## (Intercept)  
##      -3.597  
## 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.07	<0.001
Random Effects
σ²	3.29
τ₀₀ _tank	0.00
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.000 / NA

gamma distribution?

gamma requires all positive values, I have 0s here, mainly from flooring done earlier. here I will add a verrrrry small constant (0.00000000000001) to all measurements.

This looks pretty bad too which oddly means linear was best fit… Here we can see it is zero inflated and has some weird dispersion.Temp is on edge of significance here but again model output looks much worse than normal linear model.

mya$prop.change.height.1<-mya$prop.change.height+0.00000000000001
height.1 <- glmer(prop.change.height.1~temp.average*pH.normalized+(1|tank), data = mya, family = Gamma,  na.action = "na.fail")

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

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

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

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

dd1

## Global model call: glmer(formula = prop.change.height.1 ~ temp.average * pH.normalized + 
##     (1 | tank), data = mya, family = Gamma, na.action = "na.fail")
## ---
## Model selection table 
##   (Int) pH.nrm tmp.avr pH.nrm:tmp.avr df logLik    AIC delta weight
## 3 7.051        -0.3370                 4 55.739 -103.5  0.00  0.518
## 4 6.879 0.4066 -0.3326                 5 55.782 -101.6  1.91  0.199
## 1 3.932                                3 53.490 -101.0  2.50  0.148
## 8 6.668 1.0570 -0.3100       -0.07066  6 55.787  -99.6  3.91  0.073
## 2 3.716 0.6723                         4 53.610  -99.2  4.26  0.062
## Models ranked by AIC(x) 
## Random terms (all models): 
##   1 | tank

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

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

top_model1

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: prop.change.height.1 ~ temp.average + (1 | tank)
##    Data: mya
##       AIC       BIC    logLik  deviance  df.resid 
## -103.4788  -91.4362   55.7394 -111.4788       146 
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.0000  
##  Residual             0.4647  
## Number of obs: 150, groups:  tank, 16
## Fixed Effects:
##  (Intercept)  temp.average  
##        7.051        -0.337  
## optimizer (Nelder_Mead) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: prop.change.height.1 ~ temp.average + (1 | tank)
##    Data: mya
## 
##      AIC      BIC   logLik deviance df.resid 
##   -103.5    -91.4     55.7   -111.5      146 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15209 -0.73139 -0.02848  0.61134  2.51250 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.0000   0.0000  
##  Residual             0.2159   0.4647  
## Number of obs: 150, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value Pr(>|z|)    
## (Intercept)    7.0506     1.5367   4.588 4.47e-06 ***
## temp.average  -0.3370     0.1571  -2.146   0.0319 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.980
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

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.094964, p-value < 2.2e-16
## alternative hypothesis: two.sided

#plot(top_model1)
tab_model(top_model1)

	prop.change.height.1
Predictors	Estimates	CI	p
(Intercept)	1153.60	55.34 – 24046.55	<0.001
temp average	0.71	0.52 – 0.97	0.034
Random Effects
σ²	0.22
τ₀₀ _tank	0.00
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.666 / NA

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"),
                   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)", "Temp"),
                  dv.labels= "Effect of Temp and pH on mya Growth")

	Effect of Temp and pH on mya Growth
Predictors	Estimates	CI	p
(Intercept)	1153.60	55.34 – 24046.55	<0.001
Temp	0.71	0.52 – 0.97	0.034
Random Effects
σ²	0.22
τ₀₀ _tank	0.00
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.666 / NA

beta distribution?

I wonder why beta seems to work better? 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 a bit strange. Similar to binominal - best fit model includes temp. I still think linear is best fit/ probably most conservative.

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

height.1 <- glmmTMB(prop.change.height.1~temp.average*pH.normalized+(1|tank), data = mya, 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.1 ~ temp.average * pH.normalized + 
##     (1 | tank), data = mya, 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.739          +                   0.1675                      4 103.245
## 4     -2.912          +      0.4024       0.1720                      5 103.431
## 8     -4.220          +      4.3720       0.3191             -0.4531  6 104.407
## 1     -1.223          +                                               3 101.114
## 2     -1.307          +      0.2556                                   4 101.172
##      AIC delta weight
## 3 -198.5  0.00  0.430
## 4 -196.9  1.63  0.191
## 8 -196.8  1.68  0.186
## 1 -196.2  2.26  0.139
## 2 -194.3  4.15  0.054
## 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.1 ~ temp.average + (1 | tank)
## Data: mya
##       AIC       BIC    logLik  df.resid 
## -198.4898 -186.4473  103.2449       146 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups Name        Std.Dev.
##  tank   (Intercept) 0.5413  
## 
## Number of obs: 150 / Conditional model: tank, 16
## 
## Dispersion parameter for beta family ():  7.3 
## 
## Fixed Effects:
## 
## Conditional model:
##  (Intercept)  temp.average  
##      -2.7391        0.1675

summary(top_model1)

##  Family: beta  ( logit )
## Formula:          prop.change.height.1 ~ temp.average + (1 | tank)
## Data: mya
## 
##      AIC      BIC   logLik deviance df.resid 
##   -198.5   -186.4    103.2   -206.5      146 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  tank   (Intercept) 0.293    0.5413  
## Number of obs: 150, groups:  tank, 16
## 
## Dispersion parameter for beta family ():  7.3 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.73915    0.70793  -3.869 0.000109 ***
## temp.average  0.16750    0.07606   2.202 0.027657 *  
## ---
## 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.46725, p-value < 2.2e-16
## alternative hypothesis: two.sided

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(top_model1)
tab_model(top_model1)

	prop.change.height.1
Predictors	Estimates	CI	p
(Intercept)	0.06	0.02 – 0.26	<0.001
temp average	1.18	1.02 – 1.37	0.028
Random Effects
σ²	0.25
τ₀₀ _tank	0.29
ICC	0.54
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.163 / 0.612

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.06	0.02 – 0.26	<0.001
Temp	1.18	1.02 – 1.37	0.028
Random Effects
σ²	0.25
τ₀₀ _tank	0.29
ICC	0.54
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.163 / 0.612

juvenile arctica

check structure

This distribution looks pretty normal, although it looks like there is one outlier above 1 (g69). Could possibly make max growth = 1 but I think I will leave for now.

## '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  "g32" "g29" "g30" "g31" ...
##  $ Max_heigh_mm            : num  15 13.3 14.2 16.2 14.4 ...
##  $ Max_height_mm           : num  27.7 25.5 23.9 28.7 23.3 ...
##  $ percentage change height: num  84.4 92.3 67.6 77.3 61.5 ...
##  $ 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.844 0.923 0.676 0.773 0.615 ...
##  $ pH.normalized           : num  0.2 0.2 0.2 0.2 0.61 ...

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

##     ph temp  n
## 1  7.4    6  4
## 2  7.4    9  6
## 3  7.4   12  3
## 4  7.6    6  4
## 5  7.6    9  8
## 6  7.6   12  4
## 7  7.8    6  4
## 8  7.8    9  7
## 9  7.8   12  4
## 10   8    6  4
## 11   8    9 11
## 12   8   12  4

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 includes temp, which seems to make sense based on visual inspection of the data & the biology of arctica. residuals etc look pretty good but I will try a couple other families. I assume this will be best.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = juv.arctica, na.action = na.fail)

simulationOutput1<-simulateResiduals(height.1)

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

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

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = prop.change.height ~ 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
## 3 0.1139           0.06312                 4 20.757 -33.5  0.00  0.907
## 4 0.1168 -0.008728 0.06311                 5 19.328 -28.7  4.86  0.080
## 8 0.2720 -0.472700 0.04530        0.05378  6 17.829 -23.7  9.86  0.007
## 1 0.6843                                   3 14.719 -23.4 10.08  0.006
## 2 0.7163 -0.095120                         4 14.031 -20.1 13.45  0.001
## 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: prop.change.height ~ temp.average + (1 | tank)
##    Data: juv.arctica
## REML criterion at convergence: -41.5131
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.02806 
##  Residual             0.15696 
## Number of obs: 63, groups:  tank, 16
## Fixed Effects:
##  (Intercept)  temp.average  
##      0.11391       0.06312

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

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.11	-0.08 – 0.31	0.250
temp average	0.06	0.04 – 0.08	<0.001
Random Effects
σ²	0.02
τ₀₀ _tank	0.00
ICC	0.03
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.383 / 0.402

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"),
                   show.values=TRUE, show.p=TRUE,
                   title="Effect of Temp on Juvenile Arctica Growth")+
  theme_classic()

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

	Effect of Temp on Juvenile Arctica Growth
Predictors	Estimates	CI	p
(Intercept)	0.11	-0.08 – 0.31	0.250
Temp	0.06	0.04 – 0.08	<0.001
Random Effects
σ²	0.02
τ₀₀ _tank	0.00
ICC	0.03
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.383 / 0.402

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more

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

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

binomial

Wont do this here do to values over 1

gamma distribution?

gamma requires all positive values, I have no 0s here.

‘NAs producedSimulations from your fitted model produce infinite values’??? This looks really similar to linear. DHARMa cant seem to handle this for some reason so hard to fully assess whats going on. I can dig into this more or just assess whether it is a similar output to the ‘normal’ linear model.

main difference here I believe is that the estimate increases at increasing temp? vs strickly linear.

#juv.arctica$prop.change.height.1<-juv.arctica$prop.change.height+0.00000000000001
height.1 <- glmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = juv.arctica, family = Gamma,  na.action = "na.fail")

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

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

## Fixed term is "(Intercept)"

dd1

## Global model call: glmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = juv.arctica, family = Gamma, na.action = "na.fail")
## ---
## Model selection table 
##   (Int) pH.nrm  tmp.avr pH.nrm:tmp.avr df logLik   AIC delta weight
## 3 2.815        -0.14180                 4 30.411 -52.8  0.00  0.578
## 4 2.753 0.1637 -0.14080                 5 30.572 -51.1  1.68  0.250
## 8 2.325 1.4760 -0.09333        -0.1478  6 31.110 -50.2  2.60  0.157
## 1 1.577                                 3 25.291 -44.6  8.24  0.009
## 2 1.444 0.4114                          4 25.616 -43.2  9.59  0.005
## 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  ( inverse )
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: juv.arctica
##      AIC      BIC   logLik deviance df.resid 
## -52.8213 -44.2487  30.4106 -60.8213       59 
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.1296  
##  Residual             0.2211  
## Number of obs: 63, groups:  tank, 16
## Fixed Effects:
##  (Intercept)  temp.average  
##       2.8150       -0.1418

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: juv.arctica
## 
##      AIC      BIC   logLik deviance df.resid 
##    -52.8    -44.2     30.4    -60.8       59 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6972 -0.5803  0.0542  0.5844  3.7644 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.01681  0.1296  
##  Residual             0.04886  0.2211  
## Number of obs: 63, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value Pr(>|z|)    
## (Intercept)    2.8150     0.3098   9.087  < 2e-16 ***
## temp.average  -0.1418     0.0326  -4.350 1.36e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.978

#simulationOutput.top<-simulateResiduals(top_model1)
#plot(simulationOutput.top)
#testDispersion(simulationOutput.top)
#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	16.69	8.98 – 31.03	<0.001
temp average	0.87	0.81 – 0.93	<0.001
Random Effects
σ²	0.05
τ₀₀ _tank	0.02
ICC	0.26
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.548 / 0.664

#testZeroInflation(top_model1)
#plotResiduals(top_model1, juv.arctica$ph)
#plotResiduals(top_model1, juv.arctica$temp)

plot_model(top_model1, 
                   axis.labels=c("Temp"),
                   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)", "Temp"),
                  dv.labels= "Effect of Temp and pH on juv.arctica Growth")

	Effect of Temp and pH on juv.arctica Growth
Predictors	Estimates	CI	p
(Intercept)	16.69	8.98 – 31.03	<0.001
Temp	0.87	0.81 – 0.93	<0.001
Random Effects
σ²	0.05
τ₀₀ _tank	0.02
ICC	0.26
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.548 / 0.664

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more

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

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

beta distribution?

skipping since values over 1.

scallop

check structure

This distribution looks pretty normal, although it looks quite kurtose. I think all groups are not covered which makes this a no go for modeling. No 7.6 - 12 and no 8 - 12… I will still try but should double check with Brittany.

## '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 ...

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

##     ph temp n
## 1  7.4    6 2
## 2  7.4    9 4
## 3  7.4   12 3
## 4  7.6    6 4
## 5  7.6    9 3
## 6  7.8    6 3
## 7  7.8    9 4
## 8  7.8   12 2
## 9    8    6 3
## 10   8    9 5

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 includes temp, which seems to make sense based on visual inspection of the data. residuals etc look pretty good (QQ good, some patterning in resid vs predicted). I assume this will be best.

height.1 <- lmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = scallop, na.action = na.fail)

simulationOutput1<-simulateResiduals(height.1)

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

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

## Fixed term is "(Intercept)"

dd1

## Global model call: lmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = scallop, na.action = na.fail)
## ---
## Model selection table 
##     (Int)   pH.nrm tmp.avr pH.nrm:tmp.avr df logLik   AIC delta weight
## 3 0.11030          0.06713                 4 15.554 -23.1  0.00  0.817
## 4 0.06707  0.07478 0.06929                 5 14.502 -19.0  4.10  0.105
## 1 0.69220                                  3 11.964 -17.9  5.18  0.061
## 2 0.70540 -0.03961                         4 11.195 -14.4  8.72  0.010
## 8 0.10890 -0.06172 0.06465        0.01594  6 12.664 -13.3  9.78  0.006
## 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: prop.change.height ~ temp.average + (1 | tank)
##    Data: scallop
## REML criterion at convergence: -31.1076
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.05249 
##  Residual             0.11954 
## Number of obs: 33, groups:  tank, 14
## Fixed Effects:
##  (Intercept)  temp.average  
##      0.11030       0.06713

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

#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.11	-0.13 – 0.36	0.365
temp average	0.07	0.04 – 0.10	<0.001
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.16
N _tank	14
Observations	33
Marginal R² / Conditional R²	0.490 / 0.573

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, scallop$ph)

plotResiduals(top_model1, scallop$temp)

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.11	-0.13 – 0.36	0.365
Temp	0.07	0.04 – 0.10	<0.001
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.16
N _tank	14
Observations	33
Marginal R² / Conditional R²	0.490 / 0.573

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, mapping = aes(x=temp.average, y=prop.change.height, color=pH.average), size=3) + 
  #geom_point(data=effect_top_model.ph, aes(x=temp.average, y=fit), color="black") +
  geom_line(data=effect_top_model.temp, aes(x=temp.average, y=fit), color="black") +
  geom_ribbon(data= effect_top_model.temp, aes(x=temp.average, ymin=lower, ymax=upper), alpha= 0.3, fill="gray") +
  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_temp_linear.png", effect_temp, width = 10, height = 4, dpi = 300)

beta distribution?

cant because one over 1

gamma distribution?

gamma requires all positive values, I have no 0s here.

main difference here I believe is that the estimate increases at increasing temp? vs strictly linear.

seem very similar estimate to Juv arctica

#juv.arctica$prop.change.height.1<-juv.arctica$prop.change.height+0.00000000000001
height.1 <- glmer(prop.change.height~temp.average*pH.normalized+(1|tank), data = scallop, family = Gamma,  na.action = "na.fail")

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

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

## Fixed term is "(Intercept)"

dd1

## Global model call: glmer(formula = prop.change.height ~ temp.average * pH.normalized + 
##     (1 | tank), data = scallop, family = Gamma, na.action = "na.fail")
## ---
## Model selection table 
##   (Int)   pH.nrm tmp.avr pH.nrm:tmp.avr df logLik   AIC delta weight
## 3 2.985          -0.1701                 4 18.112 -28.2  0.00  0.596
## 4 3.111 -0.22220 -0.1762                 5 18.300 -26.6  1.62  0.265
## 8 3.225 -0.57920 -0.1891        0.04155  6 18.319 -24.6  3.58  0.099
## 1 1.549                                  3 14.095 -22.2  6.03  0.029
## 2 1.562 -0.04078                         4 14.097 -20.2  8.03  0.011
## 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  ( inverse )
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: scallop
##      AIC      BIC   logLik deviance df.resid 
## -28.2233 -22.2372  18.1116 -36.2233       29 
## Random effects:
##  Groups   Name        Std.Dev.
##  tank     (Intercept) 0.1594  
##  Residual             0.1916  
## Number of obs: 33, groups:  tank, 14
## Fixed Effects:
##  (Intercept)  temp.average  
##       2.9849       -0.1701

summary(top_model1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: scallop
## 
##      AIC      BIC   logLik deviance df.resid 
##    -28.2    -22.2     18.1    -36.2       29 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.36052 -0.42403 -0.00861  0.69406  2.08657 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.02542  0.1594  
##  Residual             0.03669  0.1916  
## Number of obs: 33, groups:  tank, 14
## 
## Fixed effects:
##              Estimate Std. Error t value Pr(>|z|)    
## (Intercept)   2.98487    0.44727   6.674  2.5e-11 ***
## temp.average -0.17007    0.05017  -3.390 0.000699 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.982

#simulationOutput.top<-simulateResiduals(top_model1)
#plot(simulationOutput.top)
#testDispersion(simulationOutput.top)
#plot(top_model1)
tab_model(top_model1)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	19.78	7.93 – 49.38	<0.001
temp average	0.84	0.76 – 0.93	0.002
Random Effects
σ²	0.04
τ₀₀ _tank	0.03
ICC	0.41
N _tank	14
Observations	33
Marginal R² / Conditional R²	0.629 / 0.781

#testZeroInflation(top_model1)
#plotResiduals(top_model1, juv.arctica$ph)
#plotResiduals(top_model1, juv.arctica$temp)

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

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

	Effect of Temp and pH on Scallop Growth
Predictors	Estimates	CI	p
(Intercept)	19.78	7.93 – 49.38	<0.001
Temp	0.84	0.76 – 0.93	0.002
Random Effects
σ²	0.04
τ₀₀ _tank	0.03
ICC	0.41
N _tank	14
Observations	33
Marginal R² / Conditional R²	0.629 / 0.781

formal plot

i dont think we want three lines here because we only have one predictor? might need to think through this a bit more

i dont really think we can model scallop with two groups missing.

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

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

A-Maine-Tank-2022_Bivalve-Height_V3

McMahon

2023-07-29

Set up

read in relevant files

remove dead and lost

look at data

strange specimens

look for trends

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

Individual Species Groups

mercenaria

check structure

linear model

‘normal’ linear mixed effect model

formal plot

binomial

gamma distribution?

mya

check structure

linear model

‘normal’ linear mixed effect model

formal plot

binomial

gamma distribution?

beta distribution?

juvenile arctica

check structure

linear model

‘normal’ linear mixed effect model

formal plot

binomial

gamma distribution?

formal plot

beta distribution?

scallop

check structure

linear model

‘normal’ linear mixed effect model

formal plot

beta distribution?

gamma distribution?

formal plot