this script will be the same as V4 except with changes/additions based on Tredennick et al 2021 and their script ‘02-understand-butterflies’. will cast a wide net in this and pair it down later.

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.

check the normality assumptions of full model

looks a bit funky/kurtose?

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

simulationOutput1<-simulateResiduals(height.1)

plot(height.1)

#plot(top_model1)
testZeroInflation(height.1)

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

make models and likelihood ratio test

when talking with Diana we decided that 1) there is evidence that temperature effects growth (see long lived arctica records etc.) but that we were not really testing that here. We were more interested in whether there was a relationship with pH. I will therefore test against a null model of … 1) temperature + tank 2) tank alone are we asking vs full model? or main effects?

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

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

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

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

null.m<-lmer(prop.change.height~(1|tank), data = mercenaria, na.action = na.fail)

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

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

interaction.v.main

## Data: mercenaria
## Models:
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
## interaction.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## main_effect.m    5 -249.72 -235.99 129.86  -259.72                     
## interaction.m    6 -248.26 -231.79 130.13  -260.26 0.5481  1     0.4591

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

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

main.v.temp

## Data: mercenaria
## Models:
## temp.m: prop.change.height ~ temp.average + (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## temp.m           4 -251.43 -240.45 129.72  -259.43                     
## main_effect.m    5 -249.72 -235.99 129.86  -259.72 0.2825  1      0.595

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

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

main.v.ph

## Data: mercenaria
## Models:
## ph.m: prop.change.height ~ pH.normalized + (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## ph.m             4 -223.13 -212.15 115.56  -231.13                         
## main_effect.m    5 -249.72 -235.99 129.86  -259.72 28.591  1  8.941e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

null.v.temp

## Data: mercenaria
## Models:
## null.m: prop.change.height ~ (1 | tank)
## temp.m: prop.change.height ~ temp.average + (1 | tank)
##        npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## null.m    3 -225.10 -216.87 115.55  -231.10                         
## temp.m    4 -251.43 -240.45 129.72  -259.43 28.332  1  1.022e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

main.v.null

## Data: mercenaria
## Models:
## null.m: prop.change.height ~ (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## null.m           3 -225.10 -216.87 115.55  -231.10                         
## main_effect.m    5 -249.72 -235.99 129.86  -259.72 28.614  2  6.116e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(main_effect.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##    Data: mercenaria
## 
## REML criterion at convergence: -243.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7026 -0.2263  0.0282  0.5222  2.7092 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.006423 0.08014 
##  Residual             0.004525 0.06727 
## Number of obs: 115, groups:  tank, 16
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)   -0.60480    0.10563  -5.725
## temp.average   0.08646    0.01061   8.147
## pH.normalized  0.04371    0.09307   0.470
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.937       
## pH.normalzd -0.375  0.094

summary(temp.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: mercenaria
## 
## REML criterion at convergence: -246
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.7086 -0.2398  0.0233  0.5195  2.7099 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.006016 0.07756 
##  Residual             0.004526 0.06728 
## Number of obs: 115, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  -0.58596    0.09504  -6.165
## temp.average  0.08597    0.01025   8.384
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.976

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

tab_model(main_effect.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	-0.60	-0.81 – -0.40	<0.001
temp average	0.09	0.07 – 0.11	<0.001
pH normalized	0.04	-0.14 – 0.23	0.639
Random Effects
σ²	0.00
τ₀₀ _tank	0.01
ICC	0.59
N _tank	16
Observations	115
Marginal R² / Conditional R²	0.744 / 0.894

tab_model(temp.m)

	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

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=main_effect.m)
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)

as expected, temperature was a significant predictor in the full model and a full model of main effects was significantly different than a model which didn’t include temperature. What was unexpected was that pH had no influence, a full model compared to a model which did not include pH (null of temp+tank only) were not significantly different.

another plot, ph

based on Tredennick which bases it off of https://biologyforfun.wordpress.com/2015/06/17/confidence-intervals-for-prediction-in-glmms/, should keep in mind that confidence or prediction intervals might underestimate uncertainty around our estimates.

mercenaria_vec <- seq(min(mercenaria$temp.average),max(mercenaria$temp.average))

upper_ph <- quantile(mercenaria$pH.normalized, 0.9) 
lower_ph <- quantile(mercenaria$pH.normalized, 0.1)

newdat <- tibble(temp.average = mercenaria_vec,
               Tank = mercenaria$tank[1],
               pH.normalized = lower_ph,
               ph_level = "low") %>%
  bind_rows(
    tibble(temp.average = mercenaria_vec,
           Tank = mercenaria$tank[1],
           pH.normalized = upper_ph,
           ph_level = "high")
  ) %>%
  as.data.frame()

m <- interaction.m
allout <- tibble()  # empty object
for(do_level in unique(newdat$ph_level)) {
  tmpdat <- newdat %>%
    filter(ph_level == do_level)
  
  mm <- model.matrix(~pH.normalized*temp.average, tmpdat)
  tmpdat$y <- mm%*%fixef(m) 
  # predict(m, newdat, re.form=NA) would give the same results
  pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm))
  
  outdat <- data.frame(tmpdat,
                       lower_ci = tmpdat$y-1.96*sqrt(pvar1),
                       upper_ci = tmpdat$y+1.96*sqrt(pvar1))
  allout <- bind_rows(allout, outdat)
}

mercenaria.int.ci.plot<-ggplot(allout, aes(x = temp.average, y = y,
                   color = ph_level, fill = ph_level)) +
  geom_hline(aes(yintercept = 0), linetype = 3) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), color = NA, alpha = 0.25) +
  geom_line() +
  scale_color_manual(values = rev(c("coral","dodgerblue")), 
                     labels = c("High pH", "Low pH"),
                     name = NULL) +
  scale_fill_manual(values = rev(c("coral","dodgerblue")),
                    labels = c("High pH", "Low pH"),
                    name = NULL) +
  #scale_y_continuous(breaks = c(0,1)) +
  #scale_x_continuous(breaks = c(6,9,12)) +
  labs(x = "Mean Temperature", y = "Proportional Change in Height") +
  theme_classic() +
  theme(legend.position = c(0.8, 0.85))
mercenaria.int.ci.plot

ggsave("02_output/02_plots/height_mercenaria_none_int.ci.plot.png", height=3.2, width=4, units = "in")

okay so this plot shows there is no interaction, but is otherwise pretty unhelpful, I am also a bit concerned that the axis scales don’t make sense. I have a feeling all our plots would look like this/the same bc it was only siginificant in one case, therfore I am not sure how helpful it is.

another plot, temp

mercenaria_vec <- seq(min(mercenaria$pH.normalized),max(mercenaria$pH.normalized), length.out=4)

upper_temp <- quantile(mercenaria$temp.average, 0.9) 
lower_temp <- quantile(mercenaria$temp.average, 0.1)

newdat <- tibble(pH.normalized = mercenaria_vec,
               temp.average = lower_temp,
               temp_level = "low") %>%
  bind_rows(
    tibble(pH.normalized = mercenaria_vec,
           temp.average = upper_temp,
           temp_level = "high")
  ) %>%
  as.data.frame()

m <- main_effect.m
allout <- tibble()  # empty object
for(do_level in unique(newdat$temp_level)) {
  tmpdat <- newdat %>%
    filter(temp_level == do_level)
  
  mm <- model.matrix(~pH.normalized+temp.average, tmpdat)
  tmpdat$y <- mm%*%fixef(m) 
  # predict(m, newdat, re.form=NA) would give the same results
  pvar1 <- diag(mm %*% tcrossprod(vcov(m),mm))
  
  outdat <- data.frame(tmpdat,
                       lower_ci = tmpdat$y-1.96*sqrt(pvar1),
                       upper_ci = tmpdat$y+1.96*sqrt(pvar1))
  allout <- bind_rows(allout, outdat)
}

mercenaria.int.ci.plot<-ggplot(allout, aes(x = pH.normalized, y = y,
                   color = temp_level, fill = temp_level)) +
  geom_hline(aes(yintercept = 0), linetype = 3) +
  geom_ribbon(aes(ymin = lower_ci, ymax = upper_ci), color = NA, alpha = 0.25) +
  geom_line() +
  scale_color_manual(values = rev(c("coral","dodgerblue")), 
                     labels = c("High Temp", "Low Temp"),
                     name = NULL) +
  scale_fill_manual(values = rev(c("coral","dodgerblue")),
                    labels = c("High Temp", "Low Temp"),
                    name = NULL) +
  #scale_y_continuous(breaks = c(0,1)) +
  #scale_x_continuous(breaks = c(6,9,12)) +
  labs(x = "pH Normalized", y = "Proportional Change in Height") +
  theme_classic() +
  theme(legend.position = c(0.8, 0.85))
mercenaria.int.ci.plot

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)

plot(height.1)

#plot(top_model1)
testZeroInflation(height.1)

## 
##  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(height.1, mya$ph)

plotResiduals(height.1, 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()

make models and likelihood ratio test

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

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

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

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

null.m<-lmer(prop.change.height~(1|tank), data = mya, na.action = na.fail)

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

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

interaction.v.main

## Data: mya
## Models:
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
## interaction.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)
## main_effect.m    5 -215.91 -200.86 112.96  -225.91                    
## interaction.m    6 -213.91 -195.85 112.96  -225.91 0.001  1      0.975

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

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

main.v.temp

## Data: mya
## Models:
## temp.m: prop.change.height ~ temp.average + (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## temp.m           4 -217.63 -205.59 112.82  -225.63                     
## main_effect.m    5 -215.91 -200.86 112.96  -225.91 0.2775  1     0.5984

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

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

main.v.ph

## Data: mya
## Models:
## ph.m: prop.change.height ~ pH.normalized + (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## ph.m             4 -208.00 -195.95 108.00  -216.00                        
## main_effect.m    5 -215.91 -200.86 112.96  -225.91 9.9171  1   0.001637 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

null.v.temp

## Data: mya
## Models:
## null.m: prop.change.height ~ (1 | tank)
## temp.m: prop.change.height ~ temp.average + (1 | tank)
##        npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)   
## null.m    3 -209.56 -200.52 107.78  -215.56                       
## temp.m    4 -217.63 -205.59 112.82  -225.63 10.08  1   0.001499 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

main.v.null

## Data: mya
## Models:
## null.m: prop.change.height ~ (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)   
## null.m           3 -209.56 -200.52 107.78  -215.56                        
## main_effect.m    5 -215.91 -200.86 112.96  -225.91 10.357  2   0.005636 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(main_effect.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##    Data: mya
## 
## REML criterion at convergence: -206.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.29429 -0.58316  0.00947  0.67847  2.23781 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.001121 0.03348 
##  Residual             0.012438 0.11152 
## Number of obs: 150, groups:  tank, 16
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    0.072040   0.062757   1.148
## temp.average   0.021370   0.006351   3.365
## pH.normalized -0.026586   0.054867  -0.485
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.937       
## pH.normalzd -0.377  0.097

summary(temp.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: mya
## 
## REML criterion at convergence: -210.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.2953 -0.6103  0.0161  0.6556  2.2966 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  tank     (Intercept) 0.0009807 0.03132 
##  Residual             0.0124407 0.11154 
## Number of obs: 150, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  0.060265   0.056452   1.068
## temp.average 0.021708   0.006141   3.535
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.977

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

tab_model(main_effect.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.07	-0.05 – 0.20	0.253
temp average	0.02	0.01 – 0.03	0.001
pH normalized	-0.03	-0.14 – 0.08	0.629
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.08
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.118 / 0.191

tab_model(temp.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.06	-0.05 – 0.17	0.287
temp average	0.02	0.01 – 0.03	0.001
Random Effects
σ²	0.01
τ₀₀ _tank	0.00
ICC	0.07
N _tank	16
Observations	150
Marginal R² / Conditional R²	0.118 / 0.182

This looks very similar to mercenaria, interaction not significant. pH not significant, temperature significant.

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=main_effect.m)
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=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_temp

#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 = 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_temp,
             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)

plot(height.1)

#plot(top_model1)
testZeroInflation(height.1)

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

make models and likelihood ratio test

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

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

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

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

null.m<-lmer(prop.change.height~(1|tank), data = juv.arctica, na.action = na.fail)

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

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

interaction.v.main

## Data: juv.arctica
## Models:
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
## interaction.m: prop.change.height ~ temp.average * pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## main_effect.m    5 -44.763 -34.047 27.381  -54.763                     
## interaction.m    6 -44.256 -31.397 28.128  -56.256 1.4933  1     0.2217

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

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

main.v.temp

## Data: juv.arctica
## Models:
## temp.m: prop.change.height ~ temp.average + (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)
## temp.m           4 -46.763 -38.190 27.381  -54.763                    
## main_effect.m    5 -44.763 -34.047 27.381  -54.763 3e-04  1     0.9868

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

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

main.v.ph

## Data: juv.arctica
## Models:
## ph.m: prop.change.height ~ pH.normalized + (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## ph.m             4 -26.464 -17.892 17.232  -34.464                         
## main_effect.m    5 -44.763 -34.047 27.381  -54.763 20.299  1  6.625e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

null.v.temp

## Data: juv.arctica
## Models:
## null.m: prop.change.height ~ (1 | tank)
## temp.m: prop.change.height ~ temp.average + (1 | tank)
##        npar     AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)    
## null.m    3 -28.150 -21.72 17.075  -34.150                         
## temp.m    4 -46.763 -38.19 27.381  -54.763 20.613  1  5.622e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

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

main.v.null

## Data: juv.arctica
## Models:
## null.m: prop.change.height ~ (1 | tank)
## main_effect.m: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##               npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## null.m           3 -28.150 -21.720 17.075  -34.150                         
## main_effect.m    5 -44.763 -34.047 27.381  -54.763 20.613  2  3.341e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(main_effect.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: prop.change.height ~ temp.average + pH.normalized + (1 | tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: -38.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8291 -0.6244  0.0138  0.6179  3.7785 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  tank     (Intercept) 0.001295 0.03598 
##  Residual             0.024664 0.15705 
## Number of obs: 63, groups:  tank, 16
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)    0.116780   0.109861   1.063
## temp.average   0.063107   0.011044   5.714
## pH.normalized -0.008728   0.097179  -0.090
## 
## Correlation of Fixed Effects:
##             (Intr) tmp.vr
## temp.averag -0.930       
## pH.normalzd -0.380  0.077

summary(temp.m)

## Linear mixed model fit by REML ['lmerMod']
## Formula: prop.change.height ~ temp.average + (1 | tank)
##    Data: juv.arctica
## 
## REML criterion at convergence: -41.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8062 -0.6487  0.0213  0.6443  3.8253 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  tank     (Intercept) 0.0007874 0.02806 
##  Residual             0.0246356 0.15696 
## Number of obs: 63, groups:  tank, 16
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)   0.11391    0.09806   1.162
## temp.average  0.06312    0.01063   5.937
## 
## Correlation of Fixed Effects:
##             (Intr)
## temp.averag -0.977

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

tab_model(main_effect.m)

	prop.change.height
Predictors	Estimates	CI	p
(Intercept)	0.12	-0.10 – 0.34	0.292
temp average	0.06	0.04 – 0.09	<0.001
pH normalized	-0.01	-0.20 – 0.19	0.929
Random Effects
σ²	0.02
τ₀₀ _tank	0.00
ICC	0.05
N _tank	16
Observations	63
Marginal R² / Conditional R²	0.378 / 0.409

tab_model(temp.m)

	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

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=main_effect.m)
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)

Testing this way shows the same thing, just in a clearer way and without worrying about all that arbitary ‘best fit’ that didn’t make sense in our context. Juvenile arctica show the same pattern, temp is significant, pH and interaction are not.

A-Maine-Tank-2022_Bivalve-Height_V5

McMahon

2023-08-19

Set up

read in relevant files

remove dead and lost

look at data

strange specimens

look for trends

summary table

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

check the normality assumptions of full model

make models and likelihood ratio test

formal plot

another plot, ph

another plot, temp

mya

check structure

linear model

‘normal’ linear mixed effect model

make models and likelihood ratio test

formal plot

juvenile arctica

check structure

linear model

‘normal’ linear mixed effect model

make models and likelihood ratio test

formal plot