2022-Bivalve-Growth_Mortality_V2
McMahon
2023-09-12
build the appropriate data frames
read in relevant files
This is the file with the binary mortality data for all species.
mr.lrm.start<-read_csv("00_raw_data/LRM_Mortality_2.csv")Check on structure and mutate accordingly
Reminder that for survival a 1 is survived and a 0 is died. In the column ‘died’ a 0 is survived and a 1 is died.
str(mr.lrm.start)## spc_tbl_ [462 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ tank : chr [1:462] "H1" "H1" "H1" "H1" ...
## $ Sample_ID_20220415: chr [1:462] "r22" "r26" "r18" "r31" ...
## $ ph : chr [1:462] "7.6" "7.6" "7.6" "7.6" ...
## $ temp : num [1:462] 12 12 12 12 9 9 9 9 6 6 ...
## $ survived : num [1:462] NA 1 1 0 1 1 1 1 1 1 ...
## $ Species : chr [1:462] "scallop" "scallop" "scallop" "scallop" ...
## - attr(*, "spec")=
## .. cols(
## .. tank = col_character(),
## .. Sample_ID_20220415 = col_character(),
## .. ph = col_character(),
## .. temp = col_double(),
## .. survived = col_double(),
## .. Species = col_character()
## .. )
## - attr(*, "problems")=<externalptr>mr.lrm.start$ph[mr.lrm.start$ph == "ambient"] <- "8"
mr.lrm.start$ph[mr.lrm.start$ph == "amb"] <- "8"
mr.lrm.start<-subset(mr.lrm.start, mr.lrm.start$Species != "scallop") # no longer need scallops
mr.lrm.start<-mr.lrm.start %>%
mutate(ph=as.factor(ph),
temp=as.factor(temp),
tank=as.factor(tank),
Species=as.factor(Species),
died=1-survived) %>%
filter(!is.na(survived)) #nas are those that were lost
str(mr.lrm.start)#looks good ## tibble [389 × 7] (S3: tbl_df/tbl/data.frame)
## $ tank : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Sample_ID_20220415: chr [1:389] "b38" "b39" "b40" "b41" ...
## $ ph : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ temp : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 3 ...
## $ survived : num [1:389] 1 1 1 1 1 1 1 0 1 1 ...
## $ Species : Factor w/ 3 levels "juv.arctica",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ died : num [1:389] 0 0 0 0 0 0 0 1 0 0 ...#maybe use
#mr.lrm.start$survived<-ifelse(test=mr.lrm.start$survived==0, yes="died", no="survived")check out the data
There seems to be deaths across all species groups, temperatures, pH and tanks. It does seem like death was lowered at low temperature. Only have one death for juvenile arctica could cause issues later on.
xtabs(~survived+Species, data = mr.lrm.start)## Species
## survived juv.arctica mercenaria mya
## 0 1 46 11
## 1 66 115 150xtabs(~survived+ph, data = mr.lrm.start)## ph
## survived 7.4 7.6 7.8 8
## 0 17 10 18 13
## 1 80 86 78 87xtabs(~survived+temp, data = mr.lrm.start)## temp
## survived 6 9 12
## 0 6 30 22
## 1 91 166 74xtabs(~survived+tank, data = mr.lrm.start)## tank
## survived H1 H10 H11 H12 H13 H14 H15 H16 H2 H3 H4 H5 H6 H7 H8 H9
## 0 2 2 1 4 9 4 1 3 4 1 3 5 2 3 8 6
## 1 22 22 24 23 15 20 23 22 20 23 21 19 22 21 16 18add real temp and pH data
temp.ph<-read_csv("00_raw_data/A-Maine-Tank-2022_temp-ph.csv")
mr.lrm.start<-merge(mr.lrm.start, temp.ph, by.x = "tank", by.y = "Tank")
#write_csv(weight.all, "CI-Avg-Temp-pH-NA-Dead-Removed.csv")
colnames(mr.lrm.start)## [1] "tank" "Sample_ID_20220415" "ph"
## [4] "temp" "survived" "Species"
## [7] "died" "temp.average" "temp.stdev"
## [10] "pH.average.YSI" "pH.stdev.YSI" "pH.average"
## [13] "pH.stdev"make pH relative
Brittany suggested doing this as it can make models more ‘interpretable’. Ex. if we fit a model which assumes pH has an intercept at 0 then this is really not interpretable for us (since pH of 0 is not any event we would ever expect). To ‘deal’ with this we will subtract the lowest average tank pH (7.39). Moving forward a pH of 0 represents the lowest or most acidic pH observed.
mr.lrm.start$pH.normalized<-mr.lrm.start$pH.average-(min(mr.lrm.start$pH.average))make data frames based on species
mr.lrm.start.mercenaria<-mr.lrm.start %>%
filter(Species=="mercenaria")
mr.lrm.start.mya<-mr.lrm.start %>%
filter(Species=="mya")
mr.lrm.start.juv.arctica<-mr.lrm.start %>%
filter(Species=="juv.arctica")build models
mercenaria
No significant issues with this model based on residuals. I need to look into what the ’model is nearly unidentifiable means. Maybe this model works but is no better than not using any predictive factors (ex. mortality is random)?
Fitting with AIC or BIC score our best model (fitting with dredge) only includes the random effect of tank! So for mercenaria mortality was a tank effect.
str(mr.lrm.start.mercenaria)## 'data.frame': 161 obs. of 14 variables:
## $ tank : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Sample_ID_20220415: chr "b38" "b39" "b40" "b41" ...
## $ ph : Factor w/ 4 levels "7.4","7.6","7.8",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ temp : Factor w/ 3 levels "6","9","12": 3 3 3 3 3 3 3 3 3 3 ...
## $ survived : num 1 1 1 1 1 1 1 0 1 1 ...
## $ Species : Factor w/ 3 levels "juv.arctica",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ died : num 0 0 0 0 0 0 0 1 0 0 ...
## $ 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.42 ...
## $ pH.average.YSI : num 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 7.57 ...
## $ pH.stdev.YSI : num 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 0.04 ...
## $ pH.average : num 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 7.59 ...
## $ pH.stdev : num 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 0.07 ...
## $ pH.normalized : num 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 ...simple.eda.ggplot(mr.lrm.start.mercenaria$survived) #looks pretty normal to me! one outlier identified
shapiro.test(mr.lrm.start.mercenaria$survived)##
## Shapiro-Wilk normality test
##
## data: mr.lrm.start.mercenaria$survived
## W = 0.56587, p-value < 2.2e-16describe(mr.lrm.start.mercenaria$survived)check the normality assumptions of full model
residuals look weird, not sure if this is a normal thing for binomial data?
glmerA<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mercenaria, family=binomial, na.action = "na.fail")
simulationOutput1<-simulateResiduals(glmerA)
plot(glmerA)
#plot(top_model1)
testZeroInflation(glmerA)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0126, p-value = 0.912
## alternative hypothesis: two.sided#plotResiduals(top_model1, mercenaria$pH)
#plotResiduals(top_model1, mercenaria$temp)LRT
all.m<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mercenaria, family=binomial, na.action = "na.fail")
null.m<-glmer(cbind(survived,died)~(1|tank), data = mr.lrm.start.mercenaria, family=binomial, na.action = "na.fail")
all.v.null<-anova(all.m, null.m) # not different
all.v.nullsummary(all.m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survived, died) ~ temp.average * pH.normalized + (1 | tank)
## Data: mr.lrm.start.mercenaria
##
## AIC BIC logLik deviance df.resid
## 198.1 213.5 -94.1 188.1 156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2849 -1.1906 0.4906 0.6327 0.8713
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.1657 0.4071
## Number of obs: 161, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.85474 1.81258 1.575 0.115
## temp.average -0.21218 0.19018 -1.116 0.265
## pH.normalized -0.70002 4.40340 -0.159 0.874
## temp.average:pH.normalized 0.09491 0.47640 0.199 0.842
##
## Correlation of Fixed Effects:
## (Intr) tmp.vr pH.nrm
## temp.averag -0.979
## pH.normalzd -0.820 0.815
## tmp.vrg:pH. 0.791 -0.820 -0.978summary(null.m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survived, died) ~ (1 | tank)
## Data: mr.lrm.start.mercenaria
##
## AIC BIC logLik deviance df.resid
## 195.0 201.2 -95.5 191.0 159
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8962 -1.0560 0.5274 0.6291 0.9470
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.2632 0.513
## Number of obs: 161, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9702 0.2255 4.302 1.69e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tab_model(all.m)| cbind(survived, died) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 17.37 | 0.50 – 606.27 | 0.115 |
| temp average | 0.81 | 0.56 – 1.17 | 0.265 |
| pH normalized | 0.50 | 0.00 – 2780.76 | 0.874 |
|
temp average × pH normalized |
1.10 | 0.43 – 2.80 | 0.842 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 tank | 0.17 | ||
| ICC | 0.05 | ||
| N tank | 16 | ||
| Observations | 161 | ||
| Marginal R2 / Conditional R2 | 0.037 / 0.084 | ||
tab_model(null.m)| cbind(survived, died) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 2.64 | 1.70 – 4.11 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 tank | 0.26 | ||
| ICC | 0.07 | ||
| N tank | 16 | ||
| Observations | 161 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.074 | ||
table
plot
library(viridis)## Loading required package: viridisLitemr.lrm.died.mercenaria<-mr.lrm.start.mercenaria %>%
filter(died==1)
mr.lrm.died.mercenaria$tank <- factor(mr.lrm.died.mercenaria$tank, levels=c('H1','H2','H3','H4','H5','H6','H7','H8', 'H9', 'H10','H11','H12', 'H13','H14','H15','H16'))
mr.lrm.died.mercenaria$ph <- factor(mr.lrm.died.mercenaria$ph, levels=c('8', '7.8', '7.6','7.4'))
mercenaria.plot.ph<-ggplot(mr.lrm.died.mercenaria, aes(x=tank, fill=ph))+
geom_bar()+
theme_classic()+
ylab("Number of Deaths")+
xlab("Tank ID")+
ggtitle("A")+
scale_fill_viridis(discrete = TRUE)+
guides(fill=guide_legend(title="pH Treatment"))+
theme(
legend.key.size = unit(0.25, "cm"),
legend.position = c(0.8, 0.9),
legend.background = element_rect(fill="NA",
color="black"),
legend.margin = margin(t=3, r=3, b=3, l=3),
axis.title = element_text(size = 16),
axis.text = element_text(size = 12)
)
mercenaria.plot.ph
mercenaria.plot.temp<-ggplot(mr.lrm.died.mercenaria, aes(x=tank, fill=temp))+
geom_bar()+
theme_classic()+
ylab("Number of Deaths")+
xlab("Tank ID")+
ggtitle("B")+
scale_fill_viridis(discrete = TRUE)+
guides(fill=guide_legend(title="Temperature Treatment"))+
theme(
legend.key.size = unit(0.25, "cm"),
legend.position = c(0.8, 0.9),
legend.background = element_rect(fill="NA",
color="black"),
legend.margin = margin(t=3, r=3, b=3, l=3),
axis.title = element_text(size = 16),
axis.text = element_text(size = 12)
)
mercenaria.plot.temp
mercenaria.tank.deaths<-ggarrange(mercenaria.plot.ph,mercenaria.plot.temp)
mercenaria.tank.deaths
ggsave("02_output/02_plots/mercenaria.tank.deaths.png", plot=mercenaria.tank.deaths, height = 4, width = 13)Another plot
plot_model(all.m, type = "pred", terms = c("temp.average","pH.normalized"))+
theme_classic()
plot_model(all.m, type = "pred", terms = c("pH.normalized","temp.average"))+
theme_classic()
During my next work session I would like to: -follow through on cross product/predicted probabilities and plotting issue. -continue with scallops, mya and juvenile artica
mya
No significant issues with this model based on residuals. I need to look into what the ’model is nearly unidentifiable means. Maybe this model works but is no better than not using any predictive factors (ex. mortality is random)?
Fitting with AIC or BIC score our best model (fitting with dredge) only includes the random effect of tank! So for mya mortality appears to be a tank effect? or completely random. I am a bit skeptical now that two models have resulted this way even those tank doesnt seem all that significant (when comparing the full model with lm with no random effect). Maybe the random term should include slope along with intercept?
check the normality assumptions of full model
residuals look weird, not sure if this is a normal thing for binomial data?
glmerA<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mya, family=binomial, na.action = "na.fail")## boundary (singular) fit: see help('isSingular')simulationOutput1<-simulateResiduals(glmerA)
plot(glmerA)
#plot(top_model1)
testZeroInflation(glmerA)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0159, p-value = 1
## alternative hypothesis: two.sided#plotResiduals(top_model1, mercenaria$pH)
#plotResiduals(top_model1, mercenaria$temp)LRT
all.m<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mya, family=binomial, na.action = "na.fail")## boundary (singular) fit: see help('isSingular')null.m<-glmer(cbind(survived,died)~(1|tank), data = mr.lrm.start.mya, family=binomial, na.action = "na.fail")
all.v.null<-anova(all.m, null.m) # not different
all.v.nullsummary(all.m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survived, died) ~ temp.average * pH.normalized + (1 | tank)
## Data: mr.lrm.start.mya
##
## AIC BIC logLik deviance df.resid
## 79.3 94.7 -34.7 69.3 156
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -7.9065 0.1381 0.1744 0.2624 0.5260
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 4e-14 2e-07
## Number of obs: 161, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 12.5327 5.3332 2.350 0.0188 *
## temp.average -0.9472 0.4777 -1.983 0.0474 *
## pH.normalized -10.7436 10.6204 -1.012 0.3117
## temp.average:pH.normalized 0.9635 0.9823 0.981 0.3267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) tmp.vr pH.nrm
## temp.averag -0.993
## pH.normalzd -0.915 0.915
## tmp.vrg:pH. 0.888 -0.903 -0.989
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')summary(null.m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survived, died) ~ (1 | tank)
## Data: mr.lrm.start.mya
##
## AIC BIC logLik deviance df.resid
## 82.5 88.6 -39.2 78.5 159
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.7975 0.1894 0.1894 0.2633 0.4344
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0.944 0.9716
## Number of obs: 161, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.0006 0.5594 5.364 8.16e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tab_model(all.m)| cbind(survived, died) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 277247.25 | 8.00 – 9604261677.07 | 0.019 |
| temp average | 0.39 | 0.15 – 0.99 | 0.047 |
| pH normalized | 0.00 | 0.00 – 23670.83 | 0.312 |
|
temp average × pH normalized |
2.62 | 0.38 – 17.97 | 0.327 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 tank | 0.00 | ||
| N tank | 16 | ||
| Observations | 161 | ||
| Marginal R2 / Conditional R2 | 0.369 / NA | ||
tab_model(null.m)| cbind(survived, died) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 20.10 | 6.71 – 60.16 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 tank | 0.94 | ||
| ICC | 0.22 | ||
| N tank | 16 | ||
| Observations | 161 | ||
| Marginal R2 / Conditional R2 | 0.000 / 0.223 | ||
plot
#install.packages("ggpattern")
library(viridis)
mr.lrm.died.mya<-mr.lrm.start.mya %>%
filter(died==1)
#mr.lrm.died.mya$tank <- factor(mr.lrm.died.mya$tank, levels=c('H7', 'H9', 'H12', 'H13'))
#mr.lrm.died.mya$ph <- factor(mr.lrm.died.mya$ph, levels=c('8', '7.8', '7.4', 'H13'))
mya.plot.ph<-ggplot(mr.lrm.died.mya, aes(x=tank, fill=ph))+
geom_bar()+
theme_classic()+
ylab("Number of Deaths")+
xlab("Tank ID")+
ggtitle("A")+
scale_fill_viridis(discrete = TRUE)
mya.plot.ph
mya.plot.temp<-ggplot(mr.lrm.died.mya, aes(x=tank, fill=temp))+
geom_bar()+
theme_classic()+
ylab("Number of Deaths")+
xlab("Tank ID")+
ggtitle("B")+
scale_fill_viridis(discrete = TRUE)
mya.plot.temp
mya.tank.deaths<-ggarrange(mya.plot.ph,mya.plot.temp)
mya.tank.deaths
ggsave("02_output/02_plots/mya.tank.deaths.png", plot=mya.tank.deaths, height = 6, width = 12)make better plot
#original_values <- c(0, 0.2, 0.4, 0.6)
#custom_labels_x <- c("7.4", "7.6", "7.8", "8.0")
plot_model(all.m, type = "pred", terms = c("temp.average","pH.normalized"))+
theme_classic()
plot_model(all.m, type = "pred", terms = c("pH.normalized","temp.average"))+
theme_classic()
# scale_x_discrete(breaks = original_values, labels = custom_labels_x)
#scale_color_viridis( option = "plasma", guide = FALSE, breaks = c(6, 9, 12))+
#scale_fill_viridis( option= "plasma")juvenile arctica
Need to find another method for analysis - GLMER not appropriate since deaths are not spread across groups. So we will likely just have to present this as a singluar value.
check the normality assumptions of full model
residuals look weird, not sure if this is a normal thing for binomial data?
glmerA<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.juv.arctica, family=binomial, na.action = "na.fail")## boundary (singular) fit: see help('isSingular')simulationOutput1<-simulateResiduals(glmerA)
plot(glmerA)
#plot(top_model1)
testZeroInflation(glmerA)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.1312, p-value = 1
## alternative hypothesis: two.sided#plotResiduals(top_model1, mercenaria$pH)
#plotResiduals(top_model1, mercenaria$temp)LRT
all.m<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.juv.arctica, family=binomial, na.action = "na.fail")## boundary (singular) fit: see help('isSingular')null.m<-glmer(cbind(survived,died)~(1|tank), data = mr.lrm.start.juv.arctica, family=binomial, na.action = "na.fail")## boundary (singular) fit: see help('isSingular')all.v.null<-anova(all.m, null.m) # not different
all.v.nullsummary(all.m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survived, died) ~ temp.average * pH.normalized + (1 | tank)
## Data: mr.lrm.start.juv.arctica
##
## AIC BIC logLik deviance df.resid
## 14.5 25.5 -2.2 4.5 62
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7321 0.0000 0.0000 0.0000 0.5774
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0 0
## Number of obs: 67, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 220.08 2954.07 0.075 0.941
## temp.average -19.52 256.41 -0.076 0.939
## pH.normalized -1054.04 4893.62 -0.215 0.829
## temp.average:pH.normalized 161.58 3223.52 0.050 0.960
##
## Correlation of Fixed Effects:
## (Intr) tmp.vr pH.nrm
## temp.averag -0.992
## pH.normalzd -0.005 0.002
## tmp.vrg:pH. 0.020 -0.144 -0.107
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')summary(null.m)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(survived, died) ~ (1 | tank)
## Data: mr.lrm.start.juv.arctica
##
## AIC BIC logLik deviance df.resid
## 14.4 18.8 -5.2 10.4 65
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -8.1240 0.1231 0.1231 0.1231 0.1231
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 6.453e-16 2.54e-08
## Number of obs: 67, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.190 1.008 4.158 3.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')tab_model(all.m)| cbind(survived, died) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 379356347086374808708384740091506241934840281333523730465360903491748337047670246627220400373760.00 | 0.00 – Inf | 0.941 |
| temp average | 0.00 | 0.00 – 604275902057518893069296868944108314383809508420061346124674228534280654878324632365675948000549749042611545484114594052872603264840085123606483000324152288673817950373001708979613512723600263658014649194905600.00 | 0.939 |
| pH normalized | 0.00 | 0.00 – Inf | 0.829 |
|
temp average × pH normalized |
14885045988902448021031203377813201010592776724677563374192181867708416.00 | 0.00 – Inf | 0.960 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 tank | 0.00 | ||
| N tank | 16 | ||
| Observations | 67 | ||
| Marginal R2 / Conditional R2 | 1.000 / NA | ||
tab_model(null.m)| cbind(survived, died) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 66.00 | 9.16 – 475.52 | <0.001 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 tank | 0.00 | ||
| N tank | 16 | ||
| Observations | 67 | ||
| Marginal R2 / Conditional R2 | 0.000 / NA | ||
Inspect one death
mr.lrm.start.juv.arctica.death<-mr.lrm.start.juv.arctica %>%
filter(died==1)
mr.lrm.start.juv.arctica.deathThe one death was in 12 and 7.4. These were the most extreme temp and pH treatments so it’s hard not to make any conclusions about it.