O. Faveolata larvae performance under sediment treatments
Ana M. Palacio & Xaymara Serrano
2024-03-21
General project set-up
# Load all libraries and sources required to run the script
library(tidyverse)
library(ggthemes)
library(survival)
library(survminer)
#library(gtsummary)
library(ggplot2)
library(lme4)
library(multcomp)
library(multcompView)
library(emmeans)
library(effects)
library(lmerTest)
library(modelsummary)
library(performance)
library(kableExtra)
library(sjPlot)
#library(gridExtra)
library(cowplot)
ggthe_bw<-theme_bw()+
theme(panel.grid= element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)1. Survivorship 24h
Data
# Data
Survival.data1<-read.csv("Data/Survival_24h.csv", header = TRUE)
summary(Survival.data1)## Treatment Replicate Dead Swimming
## Length:47 Min. : 1.000 Min. : 0.000 Min. :31.00
## Class :character 1st Qu.: 3.000 1st Qu.: 2.000 1st Qu.:40.00
## Mode :character Median : 5.000 Median : 5.000 Median :45.00
## Mean : 5.255 Mean : 6.043 Mean :43.96
## 3rd Qu.: 7.500 3rd Qu.:10.000 3rd Qu.:48.00
## Max. :10.000 Max. :19.000 Max. :50.00
## Prop.survived Prop.survived.SQRT.ASIN.
## Min. :0.6200 Min. :0.820
## 1st Qu.:0.8000 1st Qu.:0.960
## Median :0.9000 Median :1.060
## Mean :0.8791 Mean :1.061
## 3rd Qu.:0.9600 3rd Qu.:1.130
## Max. :1.0000 Max. :1.250 Survival.data1$Replicate<-factor(Survival.data1$Replicate, ordered = F)
Survival.data1$Treatment<-factor(Survival.data1$Treatment,
levels=c("Control", "Low Reef","High Reef",
"Low Port","High Port"))
summary(Survival.data1$Treatment)## Control Low Reef High Reef Low Port High Port
## 10 8 10 9 10Possible GML and LMs
Model 1: After comparing with other models (below) this was the best model and used in the publication. Because the model was over dispersed a dispersion parameter was added (fit1.1).
## model 1: generalized mixed model
fit1<-glmer(cbind(Swimming, Dead) ~Treatment + (1|Replicate),
data=Survival.data1,family="binomial")
anova(fit1) summary(fit1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Swimming, Dead) ~ Treatment + (1 | Replicate)
## Data: Survival.data1
##
## AIC BIC logLik deviance df.resid
## 285.4 296.5 -136.7 273.4 41
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.75163 -1.20480 -0.01122 1.27217 2.32416
##
## Random effects:
## Groups Name Variance Std.Dev.
## Replicate (Intercept) 0.3545 0.5954
## Number of obs: 47, groups: Replicate, 10
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.71828 0.26015 10.449 < 2e-16 ***
## TreatmentLow Reef 0.04665 0.27783 0.168 0.8667
## TreatmentHigh Reef -0.88778 0.21637 -4.103 4.08e-05 ***
## TreatmentLow Port -0.51335 0.23426 -2.191 0.0284 *
## TreatmentHigh Port -1.23266 0.20950 -5.884 4.01e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.436
## TrtmntHghRf -0.552 0.516
## TrtmntLwPrt -0.510 0.486 0.610
## TrtmntHghPr -0.572 0.534 0.683 0.630 plot(fit1)
check_overdispersion(fit1)## # Overdispersion test
##
## dispersion ratio = 2.185
## Pearson's Chi-Squared = 89.573
## p-value = < 0.001## unique container/observation ID for dispersion parameter
Survival.data1$obs <- factor(seq_len(nrow(Survival.data1)))
Survival.data1$obs<-factor(Survival.data1$obs, ordered = F)
fit1.1<-glmer(cbind(Swimming, Dead) ~Treatment + (1|obs),
data=Survival.data1,family="binomial")
anova(fit1.1) summary(fit1.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(Swimming, Dead) ~ Treatment + (1 | obs)
## Data: Survival.data1
##
## AIC BIC logLik deviance df.resid
## 271.4 282.5 -129.7 259.4 41
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.82185 -0.35112 -0.05736 0.55957 1.38159
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.6935 0.8328
## Number of obs: 47, groups: obs, 47
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.8688 0.3379 8.491 < 2e-16 ***
## TreatmentLow Reef 0.2757 0.5121 0.538 0.59035
## TreatmentHigh Reef -0.8783 0.4493 -1.955 0.05061 .
## TreatmentLow Port -0.5976 0.4651 -1.285 0.19884
## TreatmentHigh Port -1.3504 0.4436 -3.045 0.00233 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.625
## TrtmntHghRf -0.729 0.477
## TrtmntLwPrt -0.715 0.458 0.528
## TrtmntHghPr -0.755 0.478 0.555 0.542 plot(fit1.1) # residuals get worse with this
check_overdispersion(fit1.1) # but over dispersion is fixed ## # Overdispersion test
##
## dispersion ratio = 0.413
## Pearson's Chi-Squared = 16.932
## p-value = 1## Compare binomial models with random effects
anova(fit1, fit1.1) # Model 1.1 is the best AIC and is not over dispersed
av.surv <- data.frame (anova(fit1.1))Other models tested … they were not as good as the one chosen above
## model 2: generalized model - no random effects
fit2<-glm(cbind(Swimming,Dead)~Treatment,
data=Survival.data1,family=binomial)
anova(fit2) #av.surv <- data.frame (anova(fit2))
summary(fit2)##
## Call:
## glm(formula = cbind(Swimming, Dead) ~ Treatment, family = binomial,
## data = Survival.data1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.343 -1.382 0.000 2.104 4.061
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.5867 0.1753 14.758 < 2e-16 ***
## TreatmentLow Reef 0.1648 0.2739 0.602 0.5473
## TreatmentHigh Reef -0.8677 0.2150 -4.035 5.46e-05 ***
## TreatmentLow Port -0.4845 0.2316 -2.092 0.0364 *
## TreatmentHigh Port -1.2004 0.2079 -5.774 7.74e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 241.35 on 46 degrees of freedom
## Residual deviance: 180.40 on 42 degrees of freedom
## AIC: 325.63
##
## Number of Fisher Scoring iterations: 5 check_overdispersion(fit2) # models os overdispersed## # Overdispersion test
##
## dispersion ratio = 3.668
## Pearson's Chi-Squared = 154.058
## p-value = < 0.001 par(mfrow=c(2,2))
plot(fit2)
#par(mfrow=c(2,1))
#plot_model(fit2, type = "est", show.values=T, show.p = T)
#plot_model(fit2, type = "eff", terms="Treatment")
#par(mfrow=c(1,1))
anova(fit1.1, fit2) #Model 1.1 is better and used for downstream multicomp and plotsPairwise comparisons and plot results
Multiple comparisons
# Pairwise comparisons
# To Control:
contrast(emmeans (fit1.1, ~ Treatment), "trt.vs.ctrl1", ref = 'Control')## contrast estimate SE df z.ratio p.value
## Low Reef - Control 0.276 0.512 Inf 0.538 0.9185
## High Reef - Control -0.878 0.449 Inf -1.955 0.1621
## Low Port - Control -0.598 0.465 Inf -1.285 0.5006
## High Port - Control -1.350 0.444 Inf -3.045 0.0088
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: dunnettx method for 4 tests # Among all treatments
Sw.emmc<-emmeans(fit1.1, ~Treatment)
Sw_groups<-cld(Sw.emmc)
Sw_groupsPlot results
# Plot results
odds1<-plot_model(fit1.1, type = "est", show.values=T, show.p = T,
title = "Survivorship compared to control") +theme_sjplot()
odds1
#ggsave(file="Outputs/odds_survival_R1.svg", plot=odds1, width=4, height=4)
effects1<-plot_model(fit1.1, type = "pred", terms="Treatment", colors = c("black"), title="") + # Plots estimated marginal means (or marginal effects).
ggthe_bw +
scale_y_continuous(limits = c(0.6, 1),
#breaks = seq(0,1,0.2),
#expand = c(0.01, 0.01),
name=("Survivorship proportion"))
effects1
#ggsave(file="Outputs/effects_survival.svg", plot=effects1, width=4, height=4)2. Settlement
Data
# Data
Survival.data2<-read.csv("Data/Survival_Settlement_week.csv", header = TRUE)
summary(Survival.data2)## Treatment Replicate N_Dead Prop_dead
## Length:44 Min. :1.000 Min. : 0.00 Min. :0.0000
## Class :character 1st Qu.:3.000 1st Qu.: 3.00 1st Qu.:0.2000
## Mode :character Median :5.000 Median : 4.00 Median :0.2667
## Mean :4.909 Mean : 4.25 Mean :0.2833
## 3rd Qu.:7.000 3rd Qu.: 6.00 3rd Qu.:0.4000
## Max. :9.000 Max. :10.00 Max. :0.6667
## N.alive Prop_survived Prop.survived.SQRT.ASIN. N_settled
## Min. : 5.00 Min. :0.3300 Min. :0.5800 Min. :0.000
## 1st Qu.: 9.00 1st Qu.:0.6000 1st Qu.:0.8000 1st Qu.:1.000
## Median :11.00 Median :0.7300 Median :0.9100 Median :2.000
## Mean :10.75 Mean :0.7164 Mean :0.8998 Mean :2.182
## 3rd Qu.:12.00 3rd Qu.:0.8000 3rd Qu.:0.9600 3rd Qu.:3.000
## Max. :15.00 Max. :1.0000 Max. :1.2500 Max. :6.000
## Prop.settled Prop.settled.SQRT.ASIN. Prop.settled_alive
## Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:0.0700 1st Qu.:0.2600 1st Qu.:0.1000
## Median :0.1300 Median :0.3700 Median :0.1818
## Mean :0.1455 Mean :0.3505 Mean :0.2115
## 3rd Qu.:0.2000 3rd Qu.:0.4500 3rd Qu.:0.2912
## Max. :0.4000 Max. :0.6400 Max. :0.6000 # Random factors
Survival.data2$Replicate<-factor(Survival.data2$Replicate, ordered = F)
## unique container/observation ID for dispersion parameter
Survival.data2$obs <- factor(seq_len(nrow(Survival.data2)))
Survival.data2$obs<-factor(Survival.data2$obs, ordered = F)
# Fixed factors
Survival.data2$Treatment<-factor(Survival.data2$Treatment,
levels=c("Control", "Low Reef", "High Reef", "Low Port", "High Port"))
summary(Survival.data2$Treatment)## Control Low Reef High Reef Low Port High Port
## 9 8 9 9 9 # Long format
Sur_long<-read.csv("Data/Survival_Settlement_week_long.csv", header = TRUE)
summary(Sur_long)## Treatment Replicate Category Number
## Length:132 Min. :1.000 Length:132 Min. : 0
## Class :character 1st Qu.:3.000 Class :character 1st Qu.: 2
## Mode :character Median :5.000 Mode :character Median : 4
## Mean :4.909 Mean : 5
## 3rd Qu.:7.000 3rd Qu.: 8
## Max. :9.000 Max. :14
## Proportion
## Min. :0.0000
## 1st Qu.:0.1333
## Median :0.2683
## Mean :0.3333
## 3rd Qu.:0.5333
## Max. :0.9333 Sur_long$Treatment<-factor(Sur_long$Treatment,
levels=c("Control", "Low Reef","High Reef",
"Low Port","High Port"))
Sur_long$Category<-factor(Sur_long$Category,
levels=c("Dead", "Un_settled","Settled"))Exploratory plots
1-week all outcomes stacked
# Grouped
Spread_<-ggplot(Sur_long, aes(fill=Category, y=Number, x=Treatment)) +
geom_bar(position="dodge", stat="identity")
# Grouped
Stack<-ggplot(Sur_long, aes(fill=Category, y=Number, x=Treatment)) +
geom_bar(position="fill", stat="identity") + ggthe_bw
Stack
GML and for settlement
## model 1: generalized mixed model with random effects
fit3<-glmer(cbind(N_settled, 15-N_settled) ~Treatment + (1|Replicate),
data=Survival.data2,family="binomial")
check_overdispersion(fit3)## # Overdispersion test
##
## dispersion ratio = 0.805
## Pearson's Chi-Squared = 30.609
## p-value = 0.797 fit4<-glmer(cbind(N_settled, 15-N_settled) ~Treatment + (1|obs),
data=Survival.data2,family="binomial")
check_overdispersion(fit4)## # Overdispersion test
##
## dispersion ratio = 0.805
## Pearson's Chi-Squared = 30.609
## p-value = 0.797 ## model 2: generalized model
fit5<-glm(cbind(N_settled, 15-N_settled) ~Treatment,
data=Survival.data2,family="binomial")
check_overdispersion(fit5)## # Overdispersion test
##
## dispersion ratio = 0.785
## Pearson's Chi-Squared = 30.609
## p-value = 0.829 anova(fit3, fit4, fit5, test = "Chi") # All models performed similarly, random effects were disccarded since glm has smaller AIC and less overdispersion (these differneces were small)
anova(fit5) summary(fit5)##
## Call:
## glm(formula = cbind(N_settled, 15 - N_settled) ~ Treatment, family = "binomial",
## data = Survival.data2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7066 -0.4749 -0.1409 0.4762 1.6182
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.29614 0.20956 -6.185 6.21e-10 ***
## TreatmentLow Reef 0.01095 0.30499 0.036 0.97137
## TreatmentHigh Reef -0.71039 0.33885 -2.096 0.03604 *
## TreatmentLow Port -1.03113 0.36794 -2.802 0.00507 **
## TreatmentHigh Port -0.94293 0.35921 -2.625 0.00866 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 54.527 on 43 degrees of freedom
## Residual deviance: 37.789 on 39 degrees of freedom
## AIC: 145.64
##
## Number of Fisher Scoring iterations: 4 #plot(fit4)
av_sett <- data.frame(anova(fit5))
#confint(fit4)
plot_model(fit5, type = "est", show.values=T, show.p = T)
plot_model(fit5, type = "eff", terms="Treatment") 
Pairwise comparisons and plot results
Multiple comparisons
# Pairwise comparisons
# To Control:
contrast(emmeans (fit5, ~ Treatment), "trt.vs.ctrl1", ref = 'Control')## contrast estimate SE df z.ratio p.value
## Low Reef - Control 0.0109 0.305 Inf 0.036 0.9999
## High Reef - Control -0.7104 0.339 Inf -2.096 0.1193
## Low Port - Control -1.0311 0.368 Inf -2.802 0.0187
## High Port - Control -0.9429 0.359 Inf -2.625 0.0313
##
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: dunnettx method for 4 tests # Among all treatments
Set.emmc<-emmeans(fit5, ~Treatment)
Set_groups<-cld(Set.emmc)
Set_groupsPlot results
# Plot results
odds2<-plot_model(fit5, type = "est", show.values=T, show.p = T,
title = "Settlement compared to control") +theme_sjplot()
odds2
#ggsave(file="Outputs/odds_survival_R1.svg", plot=odds1, width=4, height=4)
effects2<-plot_model(fit5, type = "pred", terms="Treatment", colors = c("black"), title="") + # Plots estimated marginal means (or marginal effects).
ggthe_bw +
scale_y_continuous(limits = c(0, 0.35),
#breaks = seq(0,1,0.2),
#expand = c(0.01, 0.01),
name=("Settlement proportion"))
effects2
#ggsave(file="Outputs/effects_survival_R1.svg", plot=effects1, width=4, height=4)Final models
# Anova models into one data frame/table
drop1(fit1.1,test="Chisq") av.surv %>% kable("html", digits=2,
caption="Larvae 24 survivorship") %>%
kable_styling(bootstrap_options = "striped", full_width = F)| npar | Sum.Sq | Mean.Sq | F.value | |
|---|---|---|---|---|
| Treatment | 4 | 15.93 | 3.98 | 3.98 |
drop1(fit5,test="Chisq") av_sett %>% kable("html", digits=2,
caption="Settlement (1 week)") %>%
kable_styling(bootstrap_options = "striped", full_width = F) %>%
pack_rows(., "Settlement (1 week)", 1, 2) | Df | Deviance | Resid..Df | Resid..Dev | |
|---|---|---|---|---|
| Settlement (1 week) | ||||
| NULL | NA | NA | 43 | 54.53 |
| Treatment | 4 | 16.74 | 39 | 37.79 |
#To get the significance for the overall model
#Survivorship
#paste("Survivorship (24 h) X2 =", 1-pchisq(241.35-180.4, 46-42))
#Settlement
#paste("Settlement (1 week) X2 =", 1-pchisq(54.53-37.79, 43-39))Larvae_models<-list("Survivorship" = fit1.1,
"Settlement" = fit5)
modelsummary(Larvae_models, stars = TRUE,
#statistic = c('std.error', 'p.value', 'conf.int'),
title = 'Larvae performance models outputs')| Survivorship | Settlement | |
|---|---|---|
| (Intercept) | 2.869*** | −1.296*** |
| (0.338) | (0.210) | |
| TreatmentLow Reef | 0.276 | 0.011 |
| (0.512) | (0.305) | |
| TreatmentHigh Reef | −0.878+ | −0.710* |
| (0.449) | (0.339) | |
| TreatmentLow Port | −0.598 | −1.031** |
| (0.465) | (0.368) | |
| TreatmentHigh Port | −1.350** | −0.943** |
| (0.444) | (0.359) | |
| SD (Intercept obs) | 0.833 | |
| Num.Obs. | 47 | 44 |
| AIC | 271.4 | |
| BIC | 282.5 | |
| Log.Lik. | ||
| F | 4.065 | |
| RMSE | 0.02 | 0.08 |
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 |
#modelsummary(Larvae_models, estimate = "p.value")Figure 2 (in the manuscript)
Figure_1<- plot_grid(effects1, effects2, labels = "AUTO")
Figure_1
#ggsave(file="Outputs/Figure1_R1.svg", plot=Figure_1, width=8, height=4)Figure S3 (in the manuscript)
Figure_S1<-plot_grid(odds1, odds2, labels = "AUTO")
Figure_S1
#ggsave(file="Outputs/FigureS1_R1.svg", plot=Figure_S1, width=12, height=4)3 Respiration (non-significant)
Data
# Data
data<-read.csv("Data/Respiration_larvae_v2.csv", header = TRUE)
summary(data)## Species Date Timepoint Treatment.code
## Length:75 Length:75 Length:75 Length:75
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Treatment Replicate.vial Plate.. Well
## Length:75 Min. : 1.000 Min. :1.00 Length:75
## Class :character 1st Qu.: 3.500 1st Qu.:1.00 Class :character
## Mode :character Median : 6.000 Median :2.00 Mode :character
## Mean : 5.813 Mean :2.48
## 3rd Qu.: 8.000 3rd Qu.:3.50
## Max. :10.000 Max. :4.00
##
## nmol.min nmol.larva.min
## Length:75 Min. :0.000218
## Class :character 1st Qu.:0.000522
## Mode :character Median :0.000794
## Mean :0.000765
## 3rd Qu.:0.000949
## Max. :0.001581
## NA's :12 data$Replicate.vial<-factor(data$Replicate.vial, ordered = F)
data$Plate<-factor(data$Plate, ordered = F)
data$Well<-factor(data$Well, ordered = F)
data$Treatment<-factor(data$Treatment,
levels=c("Control", "Low Reef","High Reef",
"Low Port","High Port"))
summary(data$Treatment)## Control Low Reef High Reef Low Port High Port
## 15 15 15 15 15Basic plots - by larva
Figure S4
Larvae<- ggplot(data, aes (Treatment, nmol.larva.min)) +
#geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 , color="black")+
stat_summary(fun=mean, geom="point", size=2) +
#geom_point(shape=21)+
geom_jitter(alpha=0.5, shape=21, width = 0.2)+
theme(legend.position = "bottom")+
scale_y_continuous(#limits = c(0, 1),
# breaks = seq(0,1,0.2),
# expand = c(0.01, 0.01),
name=expression(
O[2]~consumption~per~larva~(nmol~min^-1))) +
ggthe_bw +
facet_grid(~Timepoint)
Larvae
#ggsave(file="Outputs/S4_Oxygen.svg", plot=Larvae, width=8, height=4)
LogLarvae<- ggplot(data, aes (Treatment, log(nmol.larva.min))) +
#geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
stat_summary(fun=mean, geom="point") +
#geom_point(shape=21)+
geom_jitter(alpha=0.5, shape=21)+
theme(legend.position = "bottom")+
scale_y_continuous(#limits = c(0, 1),
# breaks = seq(0,1,0.2),
# expand = c(0.01, 0.01),
name=("Log O2 [nmol/min/larva] corrected")) +
ggthe_bw +
facet_grid(~Timepoint)
LogLarvae
Summary stats
Respiration.summary<-data %>%
group_by(Timepoint) %>%
summarize(mean = mean(nmol.larva.min, na.rm = T),
sd = sd(nmol.larva.min, na.rm = T))
Respiration.summaryLMs
Raw data
# Raw data model
fit6<-lmer(nmol.larva.min ~Treatment*Timepoint +
(1|Replicate.vial) + (1|Plate..) + (1|Well), data=data, REML = TRUE)
anova(fit6) ranova(fit6) #summary(fit6)
plot(fit6)
step(fit6)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 14 363.25 -698.49
## (1 | Replicate.vial) 1 13 363.22 -700.43 0.0594 1 0.807442
## (1 | Plate..) 2 12 362.61 -701.22 1.2071 1 0.271905
## (1 | Well) 0 11 358.53 -695.05 8.1730 1 0.004252 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## Treatment:Timepoint 1 7.0410e-08 1.7600e-08 4 40.415 0.5280
## Treatment 2 7.5270e-08 1.8820e-08 4 51.268 0.5877
## Timepoint 0 1.5067e-06 1.5067e-06 1 46.564 47.5877
## Pr(>F)
## Treatment:Timepoint 0.7157
## Treatment 0.6730
## Timepoint 1.22e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## nmol.larva.min ~ Timepoint + (1 | Well) fit7<-lmer(nmol.larva.min ~ Timepoint + (1 | Well), data=data, REML = TRUE)
anova(fit7) ranova(fit7) summary(fit7)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nmol.larva.min ~ Timepoint + (1 | Well)
## Data: data
##
## REML criterion at convergence: -853.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.12476 -0.64982 -0.03057 0.51697 2.61136
##
## Random effects:
## Groups Name Variance Std.Dev.
## Well (Intercept) 1.914e-08 0.0001384
## Residual 3.166e-08 0.0001779
## Number of obs: 63, groups: Well, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 9.276e-04 4.533e-05 3.384e+01 20.462 < 2e-16 ***
## TimepointRecovery -3.189e-04 4.622e-05 4.656e+01 -6.898 1.22e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TimpntRcvry -0.528 plot(fit7)
#Pairwise comparisons
Sw.emmc<-emmeans(fit7, ~Timepoint)
Sw_groups<-cld(Sw.emmc)
Sw_groupslog transformed data
fit6.1<-lmer(log(nmol.larva.min) ~Treatment*Timepoint +
(1|Replicate.vial) + (1|Plate..) + (1|Well), data=data)
anova(fit6.1) ranova(fit6.1) summary(fit6.1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(nmol.larva.min) ~ Treatment * Timepoint + (1 | Replicate.vial) +
## (1 | Plate..) + (1 | Well)
## Data: data
##
## REML criterion at convergence: 40.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.97294 -0.54722 0.01006 0.50311 1.59479
##
## Random effects:
## Groups Name Variance Std.Dev.
## Well (Intercept) 0.043330 0.2082
## Replicate.vial (Intercept) 0.004329 0.0658
## Plate.. (Intercept) 0.012621 0.1123
## Residual 0.054356 0.2331
## Number of obs: 63, groups: Well, 20; Replicate.vial, 9; Plate.., 4
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) -6.928355 0.142129 10.319327 -48.747
## TreatmentLow Reef -0.146367 0.159184 34.428256 -0.919
## TreatmentHigh Reef -0.079119 0.149612 42.007046 -0.529
## TreatmentLow Port -0.134031 0.159212 44.120919 -0.842
## TreatmentHigh Port -0.191753 0.164605 43.515199 -1.165
## TimepointRecovery -0.464883 0.172729 6.042603 -2.691
## TreatmentLow Reef:TimepointRecovery 0.086974 0.210531 34.043876 0.413
## TreatmentHigh Reef:TimepointRecovery -0.087345 0.188238 33.693976 -0.464
## TreatmentLow Port:TimepointRecovery 0.112428 0.187338 31.977021 0.600
## TreatmentHigh Port:TimepointRecovery -0.009071 0.191971 31.784625 -0.047
## Pr(>|t|)
## (Intercept) 1.54e-13 ***
## TreatmentLow Reef 0.3642
## TreatmentHigh Reef 0.5997
## TreatmentLow Port 0.4044
## TreatmentHigh Port 0.2504
## TimepointRecovery 0.0357 *
## TreatmentLow Reef:TimepointRecovery 0.6821
## TreatmentHigh Reef:TimepointRecovery 0.6456
## TreatmentLow Port:TimepointRecovery 0.5526
## TreatmentHigh Port:TimepointRecovery 0.9626
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP TrtmHP TmpntR TLR:TR THR:TR TLP:TR
## TretmntLwRf -0.396
## TrtmntHghRf -0.566 0.364
## TrtmntLwPrt -0.569 0.348 0.604
## TrtmntHghPr -0.544 0.356 0.566 0.549
## TimpntRcvry -0.620 0.312 0.341 0.323 0.304
## TrtmntLR:TR 0.276 -0.753 -0.252 -0.233 -0.240 -0.459
## TrtmntHR:TR 0.331 -0.290 -0.618 -0.270 -0.282 -0.541 0.447
## TrtmntLP:TR 0.335 -0.282 -0.324 -0.581 -0.243 -0.529 0.421 0.487
## TrtmntHP:TR 0.324 -0.259 -0.317 -0.289 -0.604 -0.509 0.399 0.471 0.459 plot(fit6.1)
step(fit6.1)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 14 -20.348 68.696
## (1 | Replicate.vial) 1 13 -20.508 67.015 0.3185 1 0.572514
## (1 | Plate..) 2 12 -21.675 67.350 2.3354 1 0.126463
## (1 | Well) 0 11 -26.205 74.409 9.0588 1 0.002614 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Treatment:Timepoint 1 0.10893 0.02723 4 40.317 0.4167 0.7956
## Treatment 2 0.15384 0.03846 4 51.051 0.6216 0.6492
## Timepoint 0 2.91151 2.91151 1 46.709 47.1201 1.352e-08
##
## Treatment:Timepoint
## Treatment
## Timepoint ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## log(nmol.larva.min) ~ Timepoint + (1 | Well) fit7.2<-lmer(log(nmol.larva.min) ~ Timepoint + (1 | Well), data=data)
anova(fit7.2) ranova(fit7.2) summary(fit7.2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(nmol.larva.min) ~ Timepoint + (1 | Well)
## Data: data
##
## REML criterion at convergence: 30.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4136 -0.5836 0.0567 0.7333 1.6005
##
## Random effects:
## Groups Name Variance Std.Dev.
## Well (Intercept) 0.03925 0.1981
## Residual 0.06179 0.2486
## Number of obs: 63, groups: Well, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -7.02064 0.06409 33.80872 -109.535 < 2e-16 ***
## TimepointRecovery -0.44351 0.06461 46.70855 -6.864 1.35e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TimpntRcvry -0.522 plot(fit7.2)
Only time point and well are significant factors for respiration
Packages used
# Creates bibliography
#knitr::write_bib(c(.packages()), "packages.bib")
Arel-Bundock, Vincent. 2022a. “modelsummary: Data and Model Summaries in
R.” Journal of Statistical Software 103
(1): 1–23. https://doi.org/10.18637/jss.v103.i01.
———. 2022b. Modelsummary: Summary Tables and Plots for Statistical
Models and Data: Beautiful, Customizable, and Publication-Ready. https://vincentarelbundock.github.io/modelsummary/.
Arnold, Jeffrey B. 2021. Ggthemes: Extra Themes, Scales and Geoms
for Ggplot2. https://github.com/jrnold/ggthemes.
Auguie, Baptiste. 2017. gridExtra: Miscellaneous Functions for
"Grid" Graphics. https://CRAN.R-project.org/package=gridExtra.
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015.
“Fitting Linear Mixed-Effects Models Using lme4.” Journal of Statistical
Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.
Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2022.
Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.
Bates, Douglas, Martin Maechler, and Mikael Jagan. 2022. Matrix:
Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.
Fox, John. 2003. “Effect Displays in R for
Generalised Linear Models.” Journal of Statistical
Software 8 (15): 1–27. https://doi.org/10.18637/jss.v008.i15.
Fox, John, and Jangman Hong. 2009. “Effect Displays in
R for Multinomial and Proportional-Odds Logit Models:
Extensions to the effects Package.”
Journal of Statistical Software 32 (1): 1–24. https://doi.org/10.18637/jss.v032.i01.
Fox, John, and Sanford Weisberg. 2018. “Visualizing Fit and Lack
of Fit in Complex Regression Models with Predictor Effect Plots and
Partial Residuals.” Journal of Statistical Software 87
(9): 1–27. https://doi.org/10.18637/jss.v087.i09.
———. 2019. An r Companion to Applied Regression. 3rd ed.
Thousand Oaks CA: Sage. https://socialsciences.mcmaster.ca/jfox/Books/Companion/index.html.
Fox, John, Sanford Weisberg, and Brad Price. 2022. carData:
Companion to Applied Regression Data Sets. https://CRAN.R-project.org/package=carData.
Fox, John, Sanford Weisberg, Brad Price, Michael Friendly, and Jangman
Hong. 2022. Effects: Effect Displays for Linear, Generalized Linear,
and Other Models. https://CRAN.R-project.org/package=effects.
Genz, Alan, and Frank Bretz. 2009. Computation of Multivariate
Normal and t Probabilities. Lecture Notes in Statistics.
Heidelberg: Springer-Verlag.
Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, and Torsten Hothorn.
2021. Mvtnorm: Multivariate Normal and t Distributions. http://mvtnorm.R-forge.R-project.org.
Graves, Spencer, Hans-Peter Piepho, and Luciano Selzer with help from
Sundar Dorai-Raj. 2019. multcompView: Visualizations of Paired
Comparisons. https://CRAN.R-project.org/package=multcompView.
Henry, Lionel, and Hadley Wickham. 2020. Purrr: Functional
Programming Tools. https://CRAN.R-project.org/package=purrr.
Hothorn, Torsten. 2022. TH.data: TH’s Data Archive. https://CRAN.R-project.org/package=TH.data.
Hothorn, Torsten, Frank Bretz, and Peter Westfall. 2008.
“Simultaneous Inference in General Parametric Models.”
Biometrical Journal 50 (3): 346–63.
———. 2022. Multcomp: Simultaneous Inference in General Parametric
Models. https://CRAN.R-project.org/package=multcomp.
Kassambara, Alboukadel. 2020. Ggpubr: Ggplot2 Based Publication
Ready Plots. https://rpkgs.datanovia.com/ggpubr/.
Kassambara, Alboukadel, Marcin Kosinski, and Przemyslaw Biecek. 2021.
Survminer: Drawing Survival Curves Using Ggplot2. https://rpkgs.datanovia.com/survminer/index.html.
Kuznetsova, Alexandra, Per B. Brockhoff, and Rune H. B. Christensen.
2017. “lmerTest Package: Tests in
Linear Mixed Effects Models.” Journal of Statistical
Software 82 (13): 1–26. https://doi.org/10.18637/jss.v082.i13.
Kuznetsova, Alexandra, Per Bruun Brockhoff, and Rune Haubo Bojesen
Christensen. 2020. lmerTest: Tests in Linear Mixed Effects
Models. https://github.com/runehaubo/lmerTestR.
Lenth, Russell V. 2022. Emmeans: Estimated Marginal Means, Aka
Least-Squares Means. https://github.com/rvlenth/emmeans.
Lüdecke, Daniel. 2022. sjPlot: Data Visualization for Statistics in
Social Science. https://strengejacke.github.io/sjPlot/.
Müller, Kirill, and Hadley Wickham. 2022. Tibble: Simple Data
Frames. https://CRAN.R-project.org/package=tibble.
R Core Team. 2021. R: A Language and Environment for Statistical
Computing. Vienna, Austria: R Foundation for Statistical Computing.
https://www.R-project.org/.
Ripley, Brian. 2022. MASS: Support Functions and Datasets for
Venables and Ripley’s MASS. http://www.stats.ox.ac.uk/pub/MASS4/.
Sjoberg, Daniel D., Michael Curry, Joseph Larmarange, Jessica Lavery,
Karissa Whiting, and Emily C. Zabor. 2022. Gtsummary:
Presentation-Ready Data Summary and Analytic Result Tables. https://CRAN.R-project.org/package=gtsummary.
Sjoberg, Daniel D., Karissa Whiting, Michael Curry, Jessica A. Lavery,
and Joseph Larmarange. 2021. “Reproducible Summary Tables with the
Gtsummary Package.” The R Journal 13:
570–80. https://doi.org/10.32614/RJ-2021-053.
Terry M. Therneau, and Patricia M. Grambsch. 2000. Modeling Survival
Data: Extending the Cox Model. New York: Springer.
Therneau, Terry M. 2021. Survival: Survival Analysis. https://github.com/therneau/survival.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics
with s. Fourth. New York: Springer. https://www.stats.ox.ac.uk/pub/MASS4/.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data
Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
———. 2022a. Forcats: Tools for Working with Categorical Variables
(Factors). https://CRAN.R-project.org/package=forcats.
———. 2022b. Stringr: Simple, Consistent Wrappers for Common String
Operations. https://CRAN.R-project.org/package=stringr.
———. 2022c. Tidyverse: Easily Install and Load the Tidyverse.
https://CRAN.R-project.org/package=tidyverse.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy
D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019.
“Welcome to the tidyverse.”
Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen,
Kohske Takahashi, Claus Wilke, Kara Woo, Hiroaki Yutani, and Dewey
Dunnington. 2022. Ggplot2: Create Elegant Data Visualisations Using
the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2022.
Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, and Maximilian Girlich. 2022. Tidyr: Tidy Messy
Data. https://CRAN.R-project.org/package=tidyr.
Wickham, Hadley, Jim Hester, and Jennifer Bryan. 2022. Readr: Read
Rectangular Text Data. https://CRAN.R-project.org/package=readr.
Wilke, Claus O. 2020. Cowplot: Streamlined Plot Theme and Plot
Annotations for Ggplot2. https://wilkelab.org/cowplot/.
Zhu, Hao. 2021. kableExtra: Construct Complex Table with Kable and
Pipe Syntax. https://CRAN.R-project.org/package=kableExtra.