A-Maine-Tank-2022_Bivalve-Height_V4

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

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

## summary table

summary_by_species <- describeBy(height.all, group = "species")

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to min; returning Inf

## Warning in FUN(newX[, i], ...): no non-missing arguments to max; returning -Inf

summary_by_species

## 
##  Descriptive statistics by group 
## species: juv.arctica
##                          vars  n  mean    sd median trimmed   mad   min    max
## tank*                       1 66  8.35  4.65   8.00    8.31  5.93  1.00  16.00
## Sample_ID_20220415*         2 66 33.50 19.20  33.50   33.50 24.46  1.00  66.00
## Max_heigh_mm                3 63 12.94  2.12  12.66   12.97  2.15  6.64  17.35
## Max_height_mm               4 66 21.60  3.04  21.59   21.68  3.22 14.87  28.67
## percentage change height    5 63 68.27 20.16  65.30   67.28 21.07 30.89 127.86
## ph*                         6 66  2.59  1.14   3.00    2.61  1.48  1.00   4.00
## temp*                       7 66  1.98  0.69   2.00    1.98  0.00  1.00   3.00
## died*                       8  0   NaN    NA     NA     NaN    NA   Inf   -Inf
## species*                    9 66  1.00  0.00   1.00    1.00  0.00  1.00   1.00
## temp.average               10 66  9.01  1.94   9.03    8.99  0.58  6.14  12.06
## temp.stdev                 11 66  0.51  0.16   0.44    0.49  0.16  0.33   0.82
## pH.average.YSI             12 66  7.72  0.23   7.82    7.73  0.28  7.36   8.02
## pH.stdev.YSI               13 66  0.07  0.02   0.07    0.07  0.01  0.04   0.14
## pH.average                 14 66  7.74  0.23   7.81    7.74  0.28  7.39   8.05
## pH.stdev                   15 66  0.07  0.02   0.06    0.06  0.01  0.05   0.14
## prop.change.height         16 63  0.68  0.20   0.65    0.67  0.21  0.31   1.28
##                          range  skew kurtosis   se
## tank*                    15.00  0.06    -1.29 0.57
## Sample_ID_20220415*      65.00  0.00    -1.25 2.36
## Max_heigh_mm             10.71 -0.18    -0.20 0.27
## Max_height_mm            13.80 -0.16    -0.28 0.37
## percentage change height 96.97  0.52    -0.25 2.54
## ph*                       3.00 -0.10    -1.43 0.14
## temp*                     2.00  0.02    -0.93 0.08
## died*                     -Inf    NA       NA   NA
## species*                  0.00   NaN      NaN 0.00
## temp.average              5.92  0.08    -0.93 0.24
## temp.stdev                0.49  0.68    -0.99 0.02
## pH.average.YSI            0.66 -0.23    -1.41 0.03
## pH.stdev.YSI              0.10  2.07     5.18 0.00
## pH.average                0.66 -0.25    -1.35 0.03
## pH.stdev                  0.09  2.22     5.12 0.00
## prop.change.height        0.97  0.52    -0.25 0.03
## ------------------------------------------------------------ 
## species: mercenaria
##                          vars   n  mean    sd median trimmed   mad  min    max
## tank*                       1 115  8.07  4.57   8.00    8.01  5.93 1.00  16.00
## Sample_ID_20220415*         2 115 58.00 33.34  58.00   58.00 43.00 1.00 115.00
## Max_heigh_mm                3 115 12.02  1.20  11.86   11.97  1.10 9.60  15.21
## Max_height_mm               4 115 14.13  2.78  13.59   13.93  2.67 9.69  20.75
## percentage change height    5 115 18.43 19.76  13.02   15.79 19.31 0.00  71.26
## ph*                         6 115  2.50  1.13   2.00    2.51  1.48 1.00   4.00
## temp*                       7 115  1.93  0.73   2.00    1.91  1.48 1.00   3.00
## died*                       8   0   NaN    NA     NA     NaN    NA  Inf   -Inf
## species*                    9 115  1.00  0.00   1.00    1.00  0.00 1.00   1.00
## temp.average               10 115  8.88  2.07   9.03    8.83  3.75 6.14  12.06
## temp.stdev                 11 115  0.53  0.16   0.44    0.52  0.16 0.33   0.82
## pH.average.YSI             12 115  7.71  0.23   7.63    7.71  0.33 7.36   8.02
## pH.stdev.YSI               13 115  0.07  0.02   0.07    0.07  0.01 0.04   0.14
## pH.average                 14 115  7.72  0.23   7.67    7.72  0.39 7.39   8.05
## pH.stdev                   15 115  0.07  0.02   0.06    0.07  0.01 0.05   0.14
## prop.change.height         16 115  0.18  0.20   0.13    0.16  0.19 0.00   0.71
##                           range  skew kurtosis   se
## tank*                     15.00  0.05    -1.19 0.43
## Sample_ID_20220415*      114.00  0.00    -1.23 3.11
## Max_heigh_mm               5.61  0.43    -0.56 0.11
## Max_height_mm             11.06  0.58    -0.45 0.26
## percentage change height  71.26  0.97    -0.08 1.84
## ph*                        3.00  0.04    -1.41 0.11
## temp*                      2.00  0.11    -1.16 0.07
## died*                      -Inf    NA       NA   NA
## species*                   0.00   NaN      NaN 0.00
## temp.average               5.92  0.13    -1.14 0.19
## temp.stdev                 0.49  0.47    -1.22 0.02
## pH.average.YSI             0.66 -0.09    -1.43 0.02
## pH.stdev.YSI               0.10  1.95     4.38 0.00
## pH.average                 0.66 -0.12    -1.38 0.02
## pH.stdev                   0.09  2.05     4.19 0.00
## prop.change.height         0.71  0.97    -0.08 0.02
## ------------------------------------------------------------ 
## species: mya
##                          vars   n  mean    sd median trimmed   mad   min    max
## tank*                       1 150  8.40  4.56   8.00    8.39  5.93  1.00  16.00
## Sample_ID_20220415*         2 150 75.50 43.45  75.50   75.50 55.60  1.00 150.00
## Max_heigh_mm                3 150 15.32  1.22  15.29   15.28  1.14 12.07  20.94
## Max_height_mm               4 150 19.17  2.16  19.27   19.25  1.91 12.13  23.93
## percentage change height    5 150 25.43 12.29  25.87   25.42 14.31  0.00  54.09
## ph*                         6 150  2.50  1.12   2.50    2.50  0.74  1.00   4.00
## temp*                       7 150  1.95  0.69   2.00    1.93  0.00  1.00   3.00
## died*                       8   0   NaN    NA     NA     NaN    NA   Inf   -Inf
## species*                    9 150  1.00  0.00   1.00    1.00  0.00  1.00   1.00
## temp.average               10 150  8.92  1.95   9.03    8.88  0.58  6.14  12.06
## temp.stdev                 11 150  0.52  0.16   0.44    0.51  0.16  0.33   0.82
## pH.average.YSI             12 150  7.71  0.23   7.72    7.71  0.23  7.36   8.02
## pH.stdev.YSI               13 150  0.07  0.02   0.07    0.07  0.01  0.04   0.14
## pH.average                 14 150  7.72  0.23   7.74    7.72  0.22  7.39   8.05
## pH.stdev                   15 150  0.07  0.02   0.06    0.06  0.01  0.05   0.14
## prop.change.height         16 150  0.25  0.12   0.26    0.25  0.14  0.00   0.54
##                           range  skew kurtosis   se
## tank*                     15.00  0.00    -1.23 0.37
## Sample_ID_20220415*      149.00  0.00    -1.22 3.55
## Max_heigh_mm               8.87  0.73     2.51 0.10
## Max_height_mm             11.80 -0.36     0.09 0.18
## percentage change height  54.09  0.01    -0.66 1.00
## ph*                        3.00  0.00    -1.37 0.09
## temp*                      2.00  0.07    -0.93 0.06
## died*                      -Inf    NA       NA   NA
## species*                   0.00   NaN      NaN 0.00
## temp.average               5.92  0.09    -0.93 0.16
## temp.stdev                 0.49  0.59    -1.12 0.01
## pH.average.YSI             0.66 -0.14    -1.39 0.02
## pH.stdev.YSI               0.10  2.01     4.75 0.00
## pH.average                 0.66 -0.15    -1.33 0.02
## pH.stdev                   0.09  2.17     4.73 0.00
## prop.change.height         0.54  0.01    -0.66 0.01
## ------------------------------------------------------------ 
## species: scallop
##                          vars  n  mean    sd median trimmed   mad   min    max
## tank*                       1 65  8.51  4.33   9.00    8.49  5.93  1.00  16.00
## Sample_ID_20220415*         2 65 33.00 18.91  33.00   33.00 23.72  1.00  65.00
## Max_heigh_mm                3 65 32.42  3.33  32.81   32.36  3.68 25.01  40.44
## Max_height_mm               4 33 54.40  7.46  55.57   54.61  7.16 34.81  69.96
## percentage change height    5 33 67.52 18.09  68.72   68.34 17.47 28.26 103.84
## ph*                         6 65  2.46  1.12   2.00    2.45  1.48  1.00   4.00
## temp*                       7 65  1.88  0.72   2.00    1.85  1.48  1.00   3.00
## died*                       8  0   NaN    NA     NA     NaN    NA   Inf   -Inf
## species*                    9 65  1.00  0.00   1.00    1.00  0.00  1.00   1.00
## temp.average               10 65  8.71  1.99   9.00    8.63  3.71  6.14  12.06
## temp.stdev                 11 65  0.54  0.17   0.55    0.54  0.25  0.33   0.82
## pH.average.YSI             12 65  7.70  0.23   7.63    7.71  0.30  7.36   8.02
## pH.stdev.YSI               13 65  0.07  0.02   0.07    0.07  0.01  0.04   0.14
## pH.average                 14 65  7.72  0.23   7.67    7.72  0.33  7.39   8.05
## pH.stdev                   15 65  0.07  0.02   0.06    0.06  0.01  0.05   0.14
## prop.change.height         16 33  0.68  0.18   0.69    0.68  0.17  0.28   1.04
##                          range  skew kurtosis   se
## tank*                    15.00  0.01    -1.24 0.54
## Sample_ID_20220415*      64.00  0.00    -1.26 2.35
## Max_heigh_mm             15.43  0.13    -0.35 0.41
## Max_height_mm            35.15 -0.34     0.03 1.30
## percentage change height 75.58 -0.36    -0.35 3.15
## ph*                       3.00  0.03    -1.39 0.14
## temp*                     2.00  0.18    -1.09 0.09
## died*                     -Inf    NA       NA   NA
## species*                  0.00   NaN      NaN 0.00
## temp.average              5.92  0.19    -1.08 0.25
## temp.stdev                0.49  0.31    -1.51 0.02
## pH.average.YSI            0.66 -0.10    -1.41 0.03
## pH.stdev.YSI              0.10  2.20     5.61 0.00
## pH.average                0.66 -0.13    -1.36 0.03
## pH.stdev                  0.09  2.18     4.96 0.00
## prop.change.height        0.76 -0.36    -0.35 0.03

#summary_by_species <- height.all %>%
 # group_by(species) %>%
 # summarise(
  #  mean_height = mean(prop.change.height),
  #  median_height = median(prop.change.height),
  #  min_height = min(prop.change.height),
  #  max_height = max(prop.change.height),
   # N=n ()
  ##)

#write.csv(summary_by_species,'02_output/01_modified_data/height_summary_by_species.csv', row.names=F)

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

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
σ2	0.00
τ00 tank	0.01
ICC	0.57
N tank	16
Observations	115
Marginal R2 / Conditional R2	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
σ2	0.00
τ00 tank	0.01
ICC	0.57
N tank	16
Observations	115
Marginal R2 / Conditional R2	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. cant make a plot with both since pH isnt in model.

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", color = "Average pH")+
  theme_classic()+
  scale_color_viridis(direction = -1, option = "cividis")
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 = 5, height = 4, dpi = 300)

effect_temp_swap <- 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=pH.average, y=prop.change.height, color=temp.average), size=3) +
  labs(x="Average pH", y="Proportional Change in Height", color = "Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis( option = "plasma")
effect_temp_swap

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

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

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
σ2	0.01
τ00 tank	0.00
ICC	0.18
N tank	16
Observations	150
Marginal R2 / Conditional R2	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
σ2	0.01
τ00 tank	0.00
ICC	0.18
N tank	16
Observations	150
Marginal R2 / Conditional R2	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", color = "Average pH")+
  theme_classic()+
  scale_color_viridis(direction = -1, option = "cividis")
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_none, width = 5, height = 4, dpi = 300)

effect_none_pH <- 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=pH.average, y=prop.change.height, color=temp.average), size=3) +
  labs(x="Average Temperature (C)", y="Proportional Change in Height", color = "Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option = "plasma")
effect_none_pH

mya.model.plot<-ggarrange(effect_none_pH,effect_none,
             labels = c("A", "B"),
              ncol = 2, nrow = 1)
ggsave("02_output/02_plots/height_mya_linear_both_swap.png", mya.model.plot, width = 14, height = 4, dpi = 300)

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

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
σ2	0.02
τ00 tank	0.00
ICC	0.03
N tank	16
Observations	63
Marginal R2 / Conditional R2	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
σ2	0.02
τ00 tank	0.00
ICC	0.03
N tank	16
Observations	63
Marginal R2 / Conditional R2	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", color="Average pH")+
  theme_classic()+
  scale_color_viridis(direction = -1, option = "cividis")
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 = 5, height = 4, dpi = 300)

effect_temp_swap <- 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=pH.average, y=prop.change.height, color=temp.average), size=3) +
  labs(x="Average pH", y="Proportional Change in Height", color="Average Temperature (C)")+
  theme_classic()+
  scale_color_viridis(option = "plasma")
effect_temp_swap

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