2022-Bivalve-Growth_Mortality_V3
McMahon
2023-09-12
this script will be the same as 2 except with changes/additions based on discussion with Brittany to be more like Cameron 2022 methods. Would also like finalize plots.
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_ [477 × 6] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ tank : chr [1:477] "H1" "H1" "H1" "H1" ...
## $ Sample_ID_20220415: chr [1:477] "r22" "r26" "r18" "r31" ...
## $ ph : chr [1:477] "7.6" "7.6" "7.6" "7.6" ...
## $ temp : num [1:477] 12 12 12 12 9 9 9 9 6 6 ...
## $ survived : num [1:477] NA 1 1 0 1 1 1 1 1 1 ...
## $ Species : chr [1:477] "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 [475 × 7] (S3: tbl_df/tbl/data.frame)
## $ tank : Factor w/ 16 levels "H1","H10","H11",..: 1 1 1 9 9 9 9 10 10 10 ...
## $ Sample_ID_20220415: chr [1:475] "r26" "r18" "r31" "r14" ...
## $ 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 2 2 2 2 1 1 1 ...
## $ survived : num [1:475] 1 1 0 1 1 1 1 1 1 1 ...
## $ Species : Factor w/ 5 levels "ad.arctica","juv.arctica",..: 5 5 5 5 5 5 5 5 5 5 ...
## $ died : num [1:475] 0 0 1 0 0 0 0 0 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 ad.arctica juv.arctica mercenaria mya scallop
## 0 5 1 46 11 6
## 1 10 66 115 150 65xtabs(~survived+ph, data = mr.lrm.start)## ph
## survived 7.4 7.6 7.8 8
## 0 19 14 19 17
## 1 100 103 98 105xtabs(~survived+temp, data = mr.lrm.start)## temp
## survived 6 9 12
## 0 10 33 26
## 1 113 202 91xtabs(~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 3 3 2 4 10 4 2 3 5 2 4 5 2 4 8 8
## 1 25 26 30 28 21 25 28 26 24 29 25 24 27 26 21 21add 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))plot all deaths
mr.lrm.died<-mr.lrm.start %>%
filter(died==1)
species.labs <- c( `ad.arctica` = "Adult Arctica islandica",
`scallop` = "Placopecten magellanicus",
`juv.arctica` = "Juvenile Arctica islandica",
`mercenaria` = "Mercenaria mercenaria",
`mya` = "Mya arenaria"
)
mr.lrm.start$died<-as.numeric(mr.lrm.start$died)
all.dead.plot<-ggplot(mr.lrm.start,aes(x=ph, y=died, color=temp))+
geom_point(size = 2, position = position_jitterdodge(jitter.height = 0.3, dodge.width = 0.4))+
theme_classic()+
geom_abline(intercept = 0.5, slope = 0)+
#scale_color_viridis_d( option = "plasma", guide = FALSE, breaks = c(6, 9, 12))+
scale_color_viridis_d( option= "plasma")+
facet_grid(.~Species, scales = "free", labeller = as_labeller(species.labs))+
theme(axis.text.x = element_text(angle = 90))+
labs(x="pH", y="Mortality", color = "Temp (C)")+
scale_y_continuous(breaks = c(0,1))+
theme(legend.direction = "horizontal",legend.position = "bottom")
all.dead.plot
all.dead.plot.2<-ggplot(mr.lrm.start,aes(x=temp, y=died, color=ph))+
geom_point(size = 2, position = position_jitterdodge(jitter.height = 0.3, dodge.width = 0.4))+
theme_classic()+
geom_abline(intercept = 0.5, slope = 0)+
#scale_color_viridis_d( option = "plasma", guide = FALSE, breaks = c(6, 9, 12))+
scale_color_viridis_d( option= "cividis", direction = -1)+
facet_grid(.~Species, scales = "free", labeller = as_labeller(species.labs))+
theme(axis.text.x = element_text(angle = 90))+
labs(x="Temp (C)", y="Mortality", color = "pH")+
scale_y_continuous(breaks = c(0,1))+
theme(legend.direction = "horizontal",legend.position = "bottom")
all.dead.plot.2
dead.both<-ggarrange(all.dead.plot,all.dead.plot.2,
labels = c("A","B"),
ncol = 1, nrow = 2)
ggsave("02_output/02_plots/mortality_all_species.png", dead.both, width = 10, height = 8, dpi = 600)make data frames based on species
mr.lrm.start$died<-as.factor(mr.lrm.start$died)
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 "b39" "b47" "b43" "b40" ...
## $ 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/ 5 levels "ad.arctica","juv.arctica",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ died : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 2 1 1 ...
## $ 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? non normal residuals, but the rest look okay (qq lines, variance by group) so will proceed
glmerA<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mercenaria, family=binomial(link = "logit"), na.action = "na.fail")## boundary (singular) fit: see help('isSingular')simulationOutput1<-simulateResiduals(glmerA)
plot(simulationOutput1)
plot(glmerA)
shapiro.test(resid(glmerA))##
## Shapiro-Wilk normality test
##
## data: resid(glmerA)
## W = 0.67979, p-value < 2.2e-16#plot(top_model1)
testZeroInflation(glmerA)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.68387, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(glmerA, mr.lrm.start.mercenaria$temp)
plotResiduals(glmerA, mr.lrm.start.mercenaria$ph)
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?
full.m<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mercenaria, family=binomial(link = "logit"), na.action = "na.fail")## boundary (singular) fit: see help('isSingular')main.m<-glmer(cbind(survived,died)~temp.average+pH.normalized+(1|tank), data = mr.lrm.start.mercenaria, family=binomial(link = "logit"), na.action = "na.fail")## boundary (singular) fit: see help('isSingular')full.v.main<-anova(full.m, main.m, test="Chisq")
full.v.main # no dif - drops AIC, drop interaction summary(main.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
## 266.7 279.0 -129.3 258.7 157
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1815 -0.9464 0.2810 0.4374 0.5976
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0 0
## Number of obs: 161, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.06195 0.58735 0.105 0.916
## temp.average -0.07439 0.05946 -1.251 0.211
## pH.normalized 0.06683 0.51375 0.130 0.897
##
## Correlation of Fixed Effects:
## (Intr) tmp.vr
## temp.averag -0.937
## pH.normalzd -0.375 0.096
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')confint(main.m)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 0.0000000 0.35297837
## (Intercept) -1.0938824 1.21460117
## temp.average -0.1918663 0.04172405
## pH.normalized -0.9395940 1.07815743drop1(main.m,test="Chisq")## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')none are significant
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(main.m, type = "pred", terms = c("temp.average","pH.normalized"))+
theme_classic()
plot_model(main.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? quantile deviations more severe here, and more devergance from qq more deviance difference in groups
glmerA<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mya, family=binomial(link = "logit"), na.action = "na.fail")## boundary (singular) fit: see help('isSingular')simulationOutput1<-simulateResiduals(glmerA)
plot(simulationOutput1)
plot(glmerA)
shapiro.test(resid(glmerA))##
## Shapiro-Wilk normality test
##
## data: resid(glmerA)
## W = 0.45431, p-value < 2.2e-16#plot(top_model1)
testZeroInflation(glmerA)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 0.23466, p-value < 2.2e-16
## alternative hypothesis: two.sidedplotResiduals(glmerA, mr.lrm.start.mya$ph)
plotResiduals(glmerA, mr.lrm.start.mya$temp)
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?
full.m<-glmer(cbind(survived,died)~temp.average*pH.normalized+(1|tank), data = mr.lrm.start.mya, family=binomial(link = "logit"), na.action = "na.fail")## boundary (singular) fit: see help('isSingular')main.m<-glmer(cbind(survived,died)~temp.average+pH.normalized+(1|tank), data = mr.lrm.start.mya, family=binomial(link = "logit"), na.action = "na.fail")## boundary (singular) fit: see help('isSingular')full.v.main<-anova(full.m, main.m, test="Chisq")
full.v.main # no dif - drops AIC, drop interaction summary(main.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
## 243.8 256.1 -117.9 235.8 157
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.42042 -0.00922 0.09227 0.11625 0.23376
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 0 0
## Number of obs: 161, groups: tank, 16
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.44211 0.56573 0.781 0.435
## temp.average -0.06162 0.05671 -1.086 0.277
## pH.normalized -0.06567 0.49431 -0.133 0.894
##
## Correlation of Fixed Effects:
## (Intr) tmp.vr
## temp.averag -0.937
## pH.normalzd -0.388 0.109
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')confint(main.m)## Computing profile confidence intervals ...## 2.5 % 97.5 %
## .sig01 0.0000000 0.24856432
## (Intercept) -0.6657490 1.55747064
## temp.average -0.1735216 0.04924852
## pH.normalized -1.0362652 0.90467768drop1(main.m,test="Chisq")## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')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(full.m, type = "pred", terms = c("temp.average","pH.normalized"))+
theme_classic()
plot_model(full.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.
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.