A-Maine-Tank-2022_Bivalve-Boron
McMahon
2023-10-24
This script is for trying to start the analysis process with Diana to look at any associations between temperature, pH and possibly growth metrics as well.
Set up
read in relevant files
I took the df Diana sent me and read them in here.
look at data
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## look for trends
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.
## `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
describe(boron.all)## 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## vars n mean sd median trimmed mad
## SampleID* 1 37 19.00 10.82 19.00 19.00 13.34
## Mass 2 37 8.16 0.37 8.10 8.14 0.44
## pH* 3 37 2.49 1.12 2.00 2.48 1.48
## temp* 4 37 2.00 0.67 2.00 2.00 0.00
## tank* 5 37 8.46 4.58 9.00 8.45 5.93
## Temp with Arctica 6 0 NaN NA NA NaN NA
## avetemp (Tidbit) 7 37 9.11 1.96 9.12 9.11 0.36
## avepH (YSI) 8 37 7.71 0.23 7.66 7.71 0.30
## avepH 9 37 7.72 0.23 7.67 7.72 0.33
## avesalinity (YSI) 10 37 32.05 0.04 32.04 32.05 0.03
## water δ11B Borate 11 37 14.30 0.99 13.79 14.28 1.18
## δ11B 12 37 12.26 1.03 12.26 12.37 0.77
## run conc [B] ppb 13 37 19.00 4.59 19.00 18.84 4.45
## d11Bc-d11Bborate 14 37 -2.04 1.25 -2.10 -1.98 1.44
## Li/Ca 15 37 4.44 0.28 4.39 4.44 0.27
## B/Ca 16 37 22.35 4.34 21.39 22.04 4.84
## Mg/Ca 17 37 0.30 0.04 0.31 0.30 0.04
## Al/Ca 18 27 6.58 5.86 4.68 5.70 3.01
## Mn/Ca 19 37 3.26 1.33 2.82 3.10 1.07
## Ni/Ca 20 37 8.48 0.98 8.26 8.39 0.91
## Cu/Ca 21 20 1.47 2.70 0.72 0.81 0.20
## Zn/Ca 22 37 12.66 8.84 11.58 11.72 7.31
## Sr/Ca 23 37 1.58 0.08 1.58 1.58 0.06
## Cd/Ca 24 17 0.00 0.00 0.00 0.00 0.00
## Sn/Ca 25 2 75.44 32.71 75.44 75.44 34.29
## Ba/Ca 26 37 1.68 1.38 1.41 1.44 0.15
## Nd/Ca 27 13 0.01 0.02 0.00 0.01 0.00
## U/Ca 28 6 1.37 1.16 0.90 1.37 0.37
## Na/Ca 29 37 15.34 0.77 15.22 15.32 0.72
## Fe/Ca 30 32 9.44 33.13 1.72 2.76 1.36
## Co/Ca 31 9 0.01 0.01 0.00 0.01 0.00
## Max_heigh_mm_start 32 36 13.06 1.93 12.72 13.02 2.22
## Dry_weight_g_start 33 36 0.79 0.34 0.72 0.77 0.37
## Bouyant_weight_g_start 34 36 0.20 0.10 0.17 0.18 0.09
## Max_height_mm_end 35 37 21.84 2.71 21.61 21.82 2.56
## Dry_weight_g_end 36 37 3.61 1.26 3.43 3.55 1.13
## Bouyant_weight_g_end 37 37 1.00 0.35 0.96 0.98 0.33
## max_height_percent_change 38 36 68.62 19.91 62.92 68.52 22.09
## dry_wt_percent_change 39 36 391.46 159.21 333.66 382.06 136.80
## bouyant_wt_percent_change 40 36 472.28 205.31 390.98 453.62 185.79
## weekly_growth_rate 41 36 0.43 0.11 0.43 0.42 0.12
## percent_change_per_day 42 36 0.47 0.14 0.43 0.47 0.15
## d13C (per mil) from end of shell 43 9 1.27 0.22 1.31 1.27 0.16
## d18O (per mil) from end of shell 44 9 1.09 0.42 1.04 1.09 0.15
## omega_arag 45 37 1.01 0.44 0.83 1.01 0.52
## omega_calc 46 37 1.60 0.70 1.31 1.60 0.83
## min max range skew kurtosis se
## SampleID* 1.00 37.00 36.00 0.00 -1.30 1.78
## Mass 7.60 9.00 1.40 0.46 -0.60 0.06
## pH* 1.00 4.00 3.00 0.03 -1.42 0.18
## temp* 1.00 3.00 2.00 0.00 -0.81 0.11
## tank* 1.00 16.00 15.00 -0.03 -1.27 0.75
## Temp with Arctica Inf -Inf -Inf NA NA NA
## avetemp (Tidbit) 5.98 12.20 6.22 0.03 -0.80 0.32
## avepH (YSI) 7.32 8.02 0.70 -0.23 -1.28 0.04
## avepH 7.39 8.05 0.66 -0.14 -1.39 0.04
## avesalinity (YSI) 31.96 32.15 0.19 0.40 0.53 0.01
## water δ11B Borate 12.99 15.83 2.84 0.23 -1.41 0.16
## δ11B 8.21 13.71 5.50 -1.58 4.34 0.17
## run conc [B] ppb 10.00 31.00 21.00 0.44 -0.22 0.76
## d11Bc-d11Bborate -5.53 0.18 5.71 -0.53 -0.08 0.21
## Li/Ca 3.84 4.99 1.15 0.21 -0.65 0.05
## B/Ca 15.92 33.00 17.08 0.65 -0.41 0.71
## Mg/Ca 0.22 0.37 0.14 -0.31 -0.88 0.01
## Al/Ca 1.09 28.12 27.03 2.01 4.30 1.13
## Mn/Ca 1.67 7.19 5.52 1.16 0.71 0.22
## Ni/Ca 7.33 11.04 3.71 0.86 -0.27 0.16
## Cu/Ca 0.58 12.78 12.21 3.65 12.37 0.60
## Zn/Ca 1.30 35.76 34.46 0.95 0.11 1.45
## Sr/Ca 1.38 1.77 0.39 -0.13 0.20 0.01
## Cd/Ca 0.00 0.01 0.01 1.72 1.53 0.00
## Sn/Ca 52.31 98.57 46.26 0.00 -2.75 23.13
## Ba/Ca 1.20 9.77 8.56 5.42 28.79 0.23
## Nd/Ca 0.00 0.09 0.09 2.42 4.86 0.01
## U/Ca 0.53 3.65 3.12 1.17 -0.42 0.47
## Na/Ca 13.71 16.84 3.13 0.28 -0.69 0.13
## Fe/Ca 0.60 188.78 188.19 4.95 23.83 5.86
## Co/Ca 0.00 0.03 0.02 1.15 -0.53 0.00
## Max_heigh_mm_start 9.80 17.35 7.55 0.22 -0.92 0.32
## Dry_weight_g_start 0.28 1.76 1.48 0.68 0.02 0.06
## Bouyant_weight_g_start 0.06 0.57 0.51 1.44 2.70 0.02
## Max_height_mm_end 15.90 27.70 11.80 0.06 -0.30 0.44
## Dry_weight_g_end 1.37 6.93 5.56 0.61 0.31 0.21
## Bouyant_weight_g_end 0.37 1.90 1.53 0.56 0.04 0.06
## max_height_percent_change 30.89 106.37 75.47 0.14 -1.11 3.32
## dry_wt_percent_change 158.93 701.21 542.28 0.51 -1.11 26.53
## bouyant_wt_percent_change 156.39 993.22 836.83 0.71 -0.37 34.22
## weekly_growth_rate 0.26 0.62 0.36 0.23 -1.05 0.02
## percent_change_per_day 0.21 0.73 0.52 0.14 -1.11 0.02
## d13C (per mil) from end of shell 0.83 1.57 0.74 -0.64 -0.63 0.07
## d18O (per mil) from end of shell 0.45 1.78 1.33 0.39 -1.03 0.14
## omega_arag 0.38 1.66 1.28 0.18 -1.35 0.07
## omega_calc 0.60 2.64 2.03 0.18 -1.35 0.12#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
look at boron more closely
## spc_tbl_ [37 × 47] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ SampleID : chr [1:37] "g2" "g4" "g5" "g8" ...
## $ Mass : num [1:37] 8.3 8.2 8.7 7.8 7.8 8.7 8.2 7.9 8.6 7.9 ...
## $ pH : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 3 3 2 2 4 4 4 1 ...
## $ temp : Factor w/ 3 levels "6","9","12": 2 2 2 2 2 2 1 2 2 1 ...
## $ tank : Factor w/ 16 levels "h1","h10","h11",..: 9 9 15 15 11 11 3 2 2 7 ...
## $ Temp with Arctica : logi [1:37] NA NA NA NA NA NA ...
## $ avetemp (Tidbit) : num [1:37] 8.98 8.98 9.1 9.1 9.12 9.12 6.48 9.21 9.21 6 ...
## $ avepH (YSI) : num [1:37] 7.62 7.62 7.82 7.82 7.65 7.65 8.02 8 8 7.32 ...
## $ avepH : num [1:37] 7.64 7.64 7.82 7.82 7.64 7.64 8.05 8 8 7.39 ...
## $ avesalinity (YSI) : num [1:37] 32.1 32.1 32 32 32 ...
## $ water δ11B Borate : num [1:37] 13.8 13.8 14.6 14.6 13.8 ...
## $ δ11B : num [1:37] 12.7 12.4 11.5 12.8 12.9 ...
## $ run conc [B] ppb : num [1:37] 13 10 20 17 14 14 23 17 20 26 ...
## $ d11Bc-d11Bborate : num [1:37] -1.07 -1.37 -3.07 -1.81 -0.86 -0.75 -2.39 -2.32 -3.03 -0.67 ...
## $ Li/Ca : num [1:37] 4.42 4.46 4.71 4.22 4.57 ...
## $ B/Ca : num [1:37] 21.6 21.4 21.4 20.7 18 ...
## $ Mg/Ca : num [1:37] 0.301 0.297 0.301 0.315 0.344 ...
## $ Al/Ca : num [1:37] 14.16 6.79 NA NA NA ...
## $ Mn/Ca : num [1:37] 3.59 2.51 2.52 2.47 3.6 ...
## $ Ni/Ca : num [1:37] 11.04 7.36 7.33 7.65 7.42 ...
## $ Cu/Ca : num [1:37] 0.876 1.245 0.721 0.597 NA ...
## $ Zn/Ca : num [1:37] 29.44 33.22 14.3 6.65 15.69 ...
## $ Sr/Ca : num [1:37] 1.56 1.65 1.54 1.67 1.54 ...
## $ Cd/Ca : num [1:37] NA NA NA NA NA NA NA NA NA NA ...
## $ Sn/Ca : num [1:37] NA NA NA NA NA NA NA NA NA NA ...
## $ Ba/Ca : num [1:37] 1.39 1.41 1.39 1.48 1.34 ...
## $ Nd/Ca : num [1:37] NA NA NA NA NA ...
## $ U/Ca : num [1:37] NA NA NA NA NA ...
## $ Na/Ca : num [1:37] 14.9 16 16.7 14.5 15.1 ...
## $ Fe/Ca : num [1:37] 1.85 2.784 1.1 0.786 NA ...
## $ Co/Ca : num [1:37] NA NA NA NA NA NA NA NA NA NA ...
## $ Max_heigh_mm_start : num [1:37] 14.6 12.8 12.1 NA 13.9 ...
## $ Dry_weight_g_start : num [1:37] 0.977 0.861 0.729 NA 0.973 ...
## $ Bouyant_weight_g_start : num [1:37] 0.244 0.22 0.217 NA 0.154 0.321 0.571 0.209 0.181 0.126 ...
## $ Max_height_mm_end : num [1:37] 23.6 20 22.1 23 23.4 ...
## $ Dry_weight_g_end : num [1:37] 3.99 2.98 3.43 4.03 4.15 ...
## $ Bouyant_weight_g_end : num [1:37] 1.09 0.855 0.98 1.052 1.055 ...
## $ max_height_percent_change : num [1:37] 62 56.4 82.2 NA 68.7 ...
## $ dry_wt_percent_change : num [1:37] 308 246 370 NA 327 ...
## $ bouyant_wt_percent_change : num [1:37] 347 289 352 NA 585 ...
## $ weekly_growth_rate : num [1:37] 0.44 0.352 0.486 NA 0.466 ...
## $ percent_change_per_day : num [1:37] 0.428 0.389 0.567 NA 0.474 ...
## $ d13C (per mil) from end of shell: num [1:37] NA NA NA NA NA NA 1.42 1.37 1.57 NA ...
## $ d18O (per mil) from end of shell: num [1:37] NA NA NA NA NA NA 1.71 1.04 0.94 NA ...
## $ omega_arag : num [1:37] 0.795 0.795 1.103 1.103 0.83 ...
## $ omega_calc : num [1:37] 1.26 1.26 1.75 1.75 1.31 ...
## $ pH.normalized : num [1:37] 0.25 0.25 0.43 0.43 0.25 ...
## - attr(*, "spec")=
## .. cols(
## .. SampleID = col_character(),
## .. Mass = col_double(),
## .. pH = col_double(),
## .. temp = col_double(),
## .. tank = col_character(),
## .. `Temp with Arctica` = col_logical(),
## .. `avetemp (Tidbit)` = col_double(),
## .. `avepH (YSI)` = col_double(),
## .. avepH = col_double(),
## .. `avesalinity (YSI)` = col_double(),
## .. `water δ11B Borate` = col_double(),
## .. δ11B = col_double(),
## .. `run conc [B] ppb` = col_double(),
## .. `d11Bc-d11Bborate` = col_double(),
## .. `Li/Ca` = col_double(),
## .. `B/Ca` = col_double(),
## .. `Mg/Ca` = col_double(),
## .. `Al/Ca` = col_double(),
## .. `Mn/Ca` = col_double(),
## .. `Ni/Ca` = col_double(),
## .. `Cu/Ca` = col_double(),
## .. `Zn/Ca` = col_double(),
## .. `Sr/Ca` = col_double(),
## .. `Cd/Ca` = col_double(),
## .. `Sn/Ca` = col_double(),
## .. `Ba/Ca` = col_double(),
## .. `Nd/Ca` = col_double(),
## .. `U/Ca` = col_double(),
## .. `Na/Ca` = col_double(),
## .. `Fe/Ca` = col_double(),
## .. `Co/Ca` = col_double(),
## .. Max_heigh_mm_start = col_double(),
## .. Dry_weight_g_start = col_double(),
## .. Bouyant_weight_g_start = col_double(),
## .. Max_height_mm_end = col_double(),
## .. Dry_weight_g_end = col_double(),
## .. Bouyant_weight_g_end = col_double(),
## .. max_height_percent_change = col_double(),
## .. dry_wt_percent_change = col_double(),
## .. bouyant_wt_percent_change = col_double(),
## .. weekly_growth_rate = col_double(),
## .. percent_change_per_day = col_double(),
## .. `d13C (per mil) from end of shell` = col_double(),
## .. `d18O (per mil) from end of shell` = col_double(),
## .. omega_arag = col_double(),
## .. omega_calc = col_double()
## .. )
## - attr(*, "problems")=<externalptr>
##
## Shapiro-Wilk normality test
##
## data: boron.all$δ11B
## W = 0.86852, p-value = 0.0004392## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 37 12.26 1.03 12.26 12.37 0.77 8.21 13.71 5.5 -1.58 4.34 0.17linear 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?
boron.1 <- lmer(δ11B~`avetemp (Tidbit)`*pH.normalized+(1|tank), data = boron.all, na.action = na.fail)
#boron.1 <- lm(δ11B~`avetemp (Tidbit)`*pH.normalized, data = boron.all, na.action = na.fail) #for small numbers where mixed effect model is not possible
summary(boron.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: δ11B ~ `avetemp (Tidbit)` * pH.normalized + (1 | tank)
## Data: boron.all
##
## REML criterion at convergence: 103.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2643 -0.3608 0.1033 0.5673 2.0943
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.08446 0.2906
## Residual 0.88353 0.9400
## Number of obs: 37, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.0301 1.3815 15.2877 10.156 3.41e-08
## `avetemp (Tidbit)` -0.2296 0.1481 15.5407 -1.551 0.141
## pH.normalized -0.6559 3.5300 15.4922 -0.186 0.855
## `avetemp (Tidbit)`:pH.normalized 0.1800 0.3892 15.4648 0.462 0.650
##
## (Intercept) ***
## `avetemp (Tidbit)`
## pH.normalized
## `avetemp (Tidbit)`:pH.normalized
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) `(Td)` pH.nrm
## `vtm(Tdbt)` -0.975
## pH.normalzd -0.807 0.799
## `(Tdbt)`:H. 0.775 -0.807 -0.976summary(boron.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: δ11B ~ `avetemp (Tidbit)` * pH.normalized + (1 | tank)
## Data: boron.all
##
## REML criterion at convergence: 103.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2643 -0.3608 0.1033 0.5673 2.0943
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.08446 0.2906
## Residual 0.88353 0.9400
## Number of obs: 37, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.0301 1.3815 15.2877 10.156 3.41e-08
## `avetemp (Tidbit)` -0.2296 0.1481 15.5407 -1.551 0.141
## pH.normalized -0.6559 3.5300 15.4922 -0.186 0.855
## `avetemp (Tidbit)`:pH.normalized 0.1800 0.3892 15.4648 0.462 0.650
##
## (Intercept) ***
## `avetemp (Tidbit)`
## pH.normalized
## `avetemp (Tidbit)`:pH.normalized
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) `(Td)` pH.nrm
## `vtm(Tdbt)` -0.975
## pH.normalzd -0.807 0.799
## `(Tdbt)`:H. 0.775 -0.807 -0.976simulationOutput1<-simulateResiduals(boron.1)
plot(boron.1)
#plot(top_model1)
testZeroInflation(boron.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.sidedshapiro.test(resid(boron.1)) #normal resid##
## Shapiro-Wilk normality test
##
## data: resid(boron.1)
## W = 0.93405, p-value = 0.02995make 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?
all.m<-lmer(δ11B~`avetemp (Tidbit)`*pH.normalized+(1|tank), data = boron.all, na.action = na.fail)
temp.m<-lmer(δ11B~`avetemp (Tidbit)`+(1|tank), data = boron.all, na.action = na.fail)
main.m<-lmer(δ11B~`avetemp (Tidbit)`+pH.normalized+(1|tank), data = boron.all, na.action = na.fail)
ph.m<-lmer(δ11B~+pH.normalized+(1|tank), data = boron.all, na.action = na.fail)
null.m<-lmer(δ11B~+(1|tank), data = boron.all, na.action = na.fail)
all.v.temp<-anova(all.m, temp.m) ## refitting model(s) with ML (instead of REML)all.v.temp## Data: boron.all
## Models:
## temp.m: δ11B ~ `avetemp (Tidbit)` + (1 | tank)
## all.m: δ11B ~ `avetemp (Tidbit)` * pH.normalized + (1 | tank)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## temp.m 4 109.13 115.57 -50.565 101.130
## all.m 6 111.09 120.76 -49.547 99.094 2.0359 2 0.3613all.v.ph<-anova(all.m, ph.m) ## refitting model(s) with ML (instead of REML)all.v.ph## Data: boron.all
## Models:
## ph.m: δ11B ~ +pH.normalized + (1 | tank)
## all.m: δ11B ~ `avetemp (Tidbit)` * pH.normalized + (1 | tank)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## ph.m 4 111.79 118.23 -51.893 103.786
## all.m 6 111.09 120.76 -49.547 99.094 4.692 2 0.09575 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1all.v.null<-anova(all.m, null.m) ## refitting model(s) with ML (instead of REML)all.v.null## Data: boron.all
## Models:
## null.m: δ11B ~ +(1 | tank)
## all.m: δ11B ~ `avetemp (Tidbit)` * pH.normalized + (1 | tank)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## null.m 3 111.68 116.51 -52.838 105.677
## all.m 6 111.09 120.76 -49.547 99.094 6.5823 3 0.08647 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1temp.v.null<-anova(temp.m, null.m) ## refitting model(s) with ML (instead of REML)temp.v.null## Data: boron.all
## Models:
## null.m: δ11B ~ +(1 | tank)
## temp.m: δ11B ~ `avetemp (Tidbit)` + (1 | tank)
## npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
## null.m 3 111.68 116.51 -52.838 105.68
## temp.m 4 109.13 115.57 -50.565 101.13 4.5464 1 0.03299 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(all.m)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: δ11B ~ `avetemp (Tidbit)` * pH.normalized + (1 | tank)
## Data: boron.all
##
## REML criterion at convergence: 103.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2643 -0.3608 0.1033 0.5673 2.0943
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.08446 0.2906
## Residual 0.88353 0.9400
## Number of obs: 37, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.0301 1.3815 15.2877 10.156 3.41e-08
## `avetemp (Tidbit)` -0.2296 0.1481 15.5407 -1.551 0.141
## pH.normalized -0.6559 3.5300 15.4922 -0.186 0.855
## `avetemp (Tidbit)`:pH.normalized 0.1800 0.3892 15.4648 0.462 0.650
##
## (Intercept) ***
## `avetemp (Tidbit)`
## pH.normalized
## `avetemp (Tidbit)`:pH.normalized
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) `(Td)` pH.nrm
## `vtm(Tdbt)` -0.975
## pH.normalzd -0.807 0.799
## `(Tdbt)`:H. 0.775 -0.807 -0.976summary(temp.m)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: δ11B ~ `avetemp (Tidbit)` + (1 | tank)
## Data: boron.all
##
## REML criterion at convergence: 106
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4381 -0.4911 0.0167 0.5739 1.7066
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.06835 0.2614
## Residual 0.88885 0.9428
## Number of obs: 37, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 13.92860 0.80151 17.52500 17.378 1.72e-12 ***
## `avetemp (Tidbit)` -0.18381 0.08606 17.75092 -2.136 0.0469 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## `vtm(Tdbt)` -0.978#both these give essentially the same answer for coefficent of temp
tab_model(all.m)|  | δ11B | ||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 14.03 | 11.21 – 16.85 | <0.001 |
| avetemp (Tidbit) | -0.23 | -0.53 – 0.07 | 0.131 |
| pH normalized | -0.66 | -7.86 – 6.54 | 0.854 |
|
avetemp (Tidbit) × pH normalized |
0.18 | -0.61 – 0.97 | 0.647 |
| Random Effects | |||
| σ2 | 0.88 | ||
| Ï„00 tank | 0.08 | ||
| ICC | 0.09 | ||
| N tank | 16 | ||
| Observations | 37 | ||
| Marginal R2 / Conditional R2 | 0.158 / 0.231 | ||
tab_model(temp.m)|  | δ11B | ||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 13.93 | 12.30 – 15.56 | <0.001 |
| avetemp (Tidbit) | -0.18 | -0.36 – -0.01 | 0.040 |
| Random Effects | |||
| σ2 | 0.89 | ||
| Ï„00 tank | 0.07 | ||
| ICC | 0.07 | ||
| N tank | 16 | ||
| Observations | 37 | ||
| Marginal R2 / Conditional R2 | 0.120 / 0.182 | ||
not really effect of any of the tested variables, maybe temp. lets broaden, what Tredennick would call ‘exploration’ instead of inference. Lets consider the questions, which environmental or biological covariates are associated with δ11B
boron.long<-gather(boron.all, covariate, value, -tank, -δ11B, -SampleID, na.rm = TRUE) #make long, Diana will need to double check this looks good## Warning: attributes are not identical across measure variables; they will be
## droppedboron.long$value<-as.numeric(boron.long$value)
cov_cors <- boron.long %>% # from tredennick
group_by(covariate) %>%
summarise(rho = abs(cor(value, δ11B))) %>%
filter(rho > 0.3) %>%
pull(covariate)
boron.long %>% #plot covariate w rho>0.3
filter(covariate %in% cov_cors) %>%
ggplot(aes(x = value, y = δ11B))+
geom_point(size = 1, shape = 1, alpha = 0.4)+
geom_smooth(method = "lm", se = FALSE, linewidth = 0.5, color = "black")+
#stat_regline_equation(label.y = 17)+
stat_cor(label.y = 11.5)+
facet_wrap(~covariate, scales = "free", nrow = 3)+
labs(x = "Covariate value", y = "δ11B")+
theme(axis.text = element_text(size = 6))+
theme_classic()-> boron.corr.plot
boron.corr.plot ## `geom_smooth()` using formula = 'y ~ x'## Warning: Computation failed in `stat_cor()`.
## Caused by error in `cor.test.default()`:
## ! not enough finite observations
ggsave("02_output/02_plots/boron_corr_plot.png", boron.corr.plot, width = 18, height = 18, dpi = 300, bg="white")## `geom_smooth()` using formula = 'y ~ x'## Warning: Computation failed in `stat_cor()`.
## Caused by error in `cor.test.default()`:
## ! not enough finite observationscould then do modeling to drop parameters, but going to bed.