Benthic and Nektonic Macroinvertebrate Microplastic Analyses (Entomology 2026)
Grace E. Millsap
Setup: Packages
The data analysis for this research was conducted using The R Project for Statistical Computing, an open-source statistical computing language and software environment (R Core Team 2023).
Packages and Their Uses
tidyverse: data organization, cleaning, and transformation (Wickham et al. 2019)
MASS: fitting negative-binomial generalized linear models (Venables and Ripley 2002)
emmeans: calculating estimated marginal means and post hoc comparisons among categorical predictors (Lenth and Piaskowski 2026)
janitor: standardizing variable names and facilitate initial data cleaning (Firke 2024)
ggplot2: data visualization and graphical summaries (Wickham 2016)
dplyr: data filtering, variable creation, grouping, summarization, and joining datasets (Wickham et al. 2023)
tidyr: restructuring datasets between long and wide formats as required for analysis and table construction (Wickham et al. 2024)
gt: generation of summary tables (Iannone et al. 2026)
library(tidyverse)
library(MASS)
library(emmeans)
library(janitor)
library(ggplot2)
library(dplyr)
library(tidyr)
library(gt)Setup: Working Directory and CSV Files
data_dir <- "."File Information
pitfall: terrestrial pitfall trap data
benthics: macroinvertebrates caught during benthic sampling periods
cnidaria: gelatinous zooplankton caught during nektonic sampling periods
pitfall <- read.csv("mata-pitfalls.csv") %>%
clean_names() %>%
dplyr::select(where(~ !all(is.na(.))))
benthics <- read.csv("arthropods.csv") %>%
clean_names() %>%
dplyr::select(where(~ !all(is.na(.))))
cnidaria <- read.csv("cnidaria.csv") %>%
clean_names() %>%
dplyr::select(where(~ !all(is.na(.))))1. Terrestrial Pitfall Trap Macroinvertebrates
Import and clean data
pitfall <- read.csv(
file.path(data_dir, "mata-pitfalls.csv")
) %>%
clean_names() %>%
dplyr::select(where(~ !all(is.na(.))))pitfall <- pitfall %>%
mutate(
functional_group = trimws(as.character(functional_group)),
functional_group = if_else(
functional_group == "Opportinist",
"Opportunist",
functional_group
),
functional_group = factor(functional_group),
family = factor(family),
dryweight = as.numeric(dryweight),
count = as.numeric(count),
total_mp = as.numeric(total_mp),
fragments = as.numeric(fragments),
fibers = as.numeric(fibers),
films = as.numeric(films),
microbeads = as.numeric(microbeads),
dryweight_mg = dryweight * 1000
)
head(pitfall)## season site transect family functional_group count dryweight
## 1 F25 PRBMF 1 Gryllidae Opportunist 8 0.03720
## 2 F25 PRBMF 1 Formicidae Opportunist 397 0.05950
## 3 F25 PRBMF 1 Mycetophagidae Grazer 20 0.00520
## 4 F25 PRBMF 1 Carabidae Pursuit Predator 2 0.01440
## 5 F25 PRBMF 1 Scarabidae Decomposer 1 0.04240
## 6 F25 PRBMF 1 Pteromalidae Pollinator 1 0.00005
## total_mp fragments fibers films microbeads counter dryweight_mg
## 1 69 69 0 0 0 KBA 37.20
## 2 280 278 2 0 0 GEM 59.50
## 3 174 142 32 0 0 GEM 5.20
## 4 525 517 7 0 1 GEM 14.40
## 5 24 18 8 0 0 GEM 42.40
## 6 217 212 4 1 0 GEM 0.05str(pitfall)## 'data.frame': 122 obs. of 14 variables:
## $ season : chr "F25" "F25" "F25" "F25" ...
## $ site : chr "PRBMF" "PRBMF" "PRBMF" "PRBMF" ...
## $ transect : int 1 1 1 1 1 1 1 1 1 1 ...
## $ family : Factor w/ 60 levels "Acari","Anisolabididae",..: 35 30 42 10 51 47 50 5 22 56 ...
## $ functional_group: Factor w/ 10 levels "Ambush Predator",..: 6 6 5 9 3 8 9 6 4 6 ...
## $ count : num 8 397 20 2 1 1 23 1 1 2 ...
## $ dryweight : num 0.0372 0.0595 0.0052 0.0144 0.0424 0.00005 0.0119 0.0016 0.0018 0.0003 ...
## $ total_mp : num 69 280 174 525 24 217 581 12 33 35 ...
## $ fragments : num 69 278 142 517 18 212 576 12 30 27 ...
## $ fibers : num 0 2 32 7 8 4 5 0 3 8 ...
## $ films : num 0 0 0 0 0 1 0 0 0 0 ...
## $ microbeads : num 0 0 0 1 0 0 0 0 0 0 ...
## $ counter : chr "KBA" "GEM" "GEM" "GEM" ...
## $ dryweight_mg : num 37.2 59.5 5.2 14.4 42.4 0.05 11.9 1.6 1.8 0.3 ...table(pitfall$functional_group)##
## Ambush Predator Aphagous Decomposer Fluid Feeder
## 6 3 22 14
## Grazer Opportunist Parasitoid Pollinator
## 7 29 1 13
## Pursuit Predator Variable
## 15 12summary(pitfall$dryweight_mg)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05 0.10 1.05 20.68 5.65 547.801.1 Is greater total biomass indicative of higher MP load?
model_biomass_null <- glm.nb(
total_mp ~ 1,
data = pitfall
)
model_biomass <- glm.nb(
total_mp ~ dryweight_mg,
data = pitfall
)
summary(model_biomass)##
## Call:
## glm.nb(formula = total_mp ~ dryweight_mg, data = pitfall, init.theta = 1.158108545,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.747033 0.088461 53.662 < 2e-16 ***
## dryweight_mg 0.007631 0.001257 6.073 1.25e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.1581) family taken to be 1)
##
## Null deviance: 174.08 on 121 degrees of freedom
## Residual deviance: 137.65 on 120 degrees of freedom
## AIC: 1446.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.158
## Std. Err.: 0.134
##
## 2 x log-likelihood: -1440.155anova(
model_biomass_null,
model_biomass,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: total_mp
## Model theta Resid. df 2 x log-lik. Test df LR stat.
## 1 1 0.9368323 121 -1472.859
## 2 dryweight_mg 1.1581085 120 -1440.155 1 vs 2 1 32.70313
## Pr(Chi)
## 1
## 2 1.073643e-08RESULT: Arthropod dry biomass was positively associated with total MP abundance. Negative-binomial regression indicated that MP abundance increased significantly with dry weight (beta = 0.00763, pm 0.00126, SE, z = 6.07, p < 0.001); each 1 mg increase in dry biomass corresponded to an approximately 0.77% increase in expected MP abundance. Inclusion of biomass significantly improved model fit relative to the null model (chi^2 = 32.70, p < 0.001).
Figure 1.1
pitfall_pred <- data.frame(
dryweight_mg = exp(
seq(
log(min(pitfall$dryweight_mg)),
log(max(pitfall$dryweight_mg)),
length.out = 200
)
)
)
pitfall_pred$predicted_mp <- predict(
model_biomass,
newdata = pitfall_pred,
type = "response"
)
figure_6.1 <- ggplot(
pitfall,
aes(
x = dryweight_mg,
y = total_mp
)
) +
geom_point(
alpha = 0.7,
size = 2,
color = "cyan4"
) +
geom_line(
data = pitfall_pred,
aes(
x = dryweight_mg,
y = predicted_mp
),
inherit.aes = FALSE,
color = "red4",
linewidth = 1
) +
scale_x_log10() +
scale_y_log10() +
theme_classic(base_size = 13) +
labs(
x = "Dry weight (mg)",
y = "Total microplastics"
) +
theme(
axis.title = element_text(face = "bold"),
axis.text = element_text(size = 11)
)
figure_6.1
1.2 Do functional groups differ in MP load?
model_guild_mass <- glm.nb(
total_mp ~ dryweight_mg,
data = pitfall
)
model_guild_full <- glm.nb(
total_mp ~ functional_group + dryweight_mg,
data = pitfall
)
anova(
model_guild_mass,
model_guild_full,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: total_mp
## Model theta Resid. df 2 x log-lik. Test
## 1 dryweight_mg 1.158109 120 -1440.155
## 2 functional_group + dryweight_mg 1.364428 111 -1415.691 1 vs 2
## df LR stat. Pr(Chi)
## 1
## 2 9 24.46493 0.003623732summary(model_guild_full)##
## Call:
## glm.nb(formula = total_mp ~ functional_group + dryweight_mg,
## data = pitfall, init.theta = 1.364428169, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.887225 0.354192 13.798 < 2e-16 ***
## functional_groupAphagous -0.571716 0.611691 -0.935 0.3500
## functional_groupDecomposer -0.162476 0.398399 -0.408 0.6834
## functional_groupFluid Feeder -0.816911 0.423032 -1.931 0.0535 .
## functional_groupGrazer -0.518820 0.481301 -1.078 0.2811
## functional_groupOpportunist 0.271917 0.385883 0.705 0.4810
## functional_groupParasitoid -0.653680 0.934259 -0.700 0.4841
## functional_groupPollinator -0.170043 0.424589 -0.400 0.6888
## functional_groupPursuit Predator 0.287064 0.415380 0.691 0.4895
## functional_groupVariable -0.905418 0.432935 -2.091 0.0365 *
## dryweight_mg 0.005621 0.001206 4.660 3.16e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.3644) family taken to be 1)
##
## Null deviance: 204.68 on 121 degrees of freedom
## Residual deviance: 135.25 on 111 degrees of freedom
## AIC: 1439.7
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.364
## Std. Err.: 0.160
##
## 2 x log-likelihood: -1415.691Estimated marginal means and pairwise comparisons
guild.emm <- emmeans(
model_guild_full,
~ functional_group,
type = "response"
)
guild.emm## functional_group response SE df asymp.LCL asymp.UCL
## Ambush Predator 148.9 52.4 Inf 74.7 296.7
## Aphagous 84.1 42.0 Inf 31.6 223.7
## Decomposer 126.6 23.4 Inf 88.1 181.8
## Fluid Feeder 65.8 15.3 Inf 41.7 103.8
## Grazer 88.6 29.0 Inf 46.7 168.3
## Opportunist 195.5 31.2 Inf 142.9 267.3
## Parasitoid 77.5 67.0 Inf 14.2 422.0
## Pollinator 125.6 30.3 Inf 78.4 201.4
## Pursuit Predator 198.4 44.6 Inf 127.8 308.2
## Variable 60.2 15.1 Inf 36.8 98.5
##
## Confidence level used: 0.95
## Intervals are back-transformed from the log scalepairs(
guild.emm,
adjust = "tukey"
)## contrast ratio SE df null z.ratio p.value
## Ambush Predator / Aphagous 1.771 1.0800 Inf 1 0.935 0.9954
## Ambush Predator / Decomposer 1.176 0.4690 Inf 1 0.408 1.0000
## Ambush Predator / Fluid Feeder 2.263 0.9580 Inf 1 1.931 0.6477
## Ambush Predator / Grazer 1.680 0.8090 Inf 1 1.078 0.9869
## Ambush Predator / Opportunist 0.762 0.2940 Inf 1 -0.705 0.9995
## Ambush Predator / Parasitoid 1.923 1.8000 Inf 1 0.700 0.9995
## Ambush Predator / Pollinator 1.185 0.5030 Inf 1 0.400 1.0000
## Ambush Predator / Pursuit Predator 0.750 0.3120 Inf 1 -0.691 0.9996
## Ambush Predator / Variable 2.473 1.0700 Inf 1 2.091 0.5336
## Aphagous / Decomposer 0.664 0.3530 Inf 1 -0.770 0.9990
## Aphagous / Fluid Feeder 1.278 0.7030 Inf 1 0.446 1.0000
## Aphagous / Grazer 0.948 0.5650 Inf 1 -0.089 1.0000
## Aphagous / Opportunist 0.430 0.2260 Inf 1 -1.608 0.8446
## Aphagous / Parasitoid 1.085 1.0800 Inf 1 0.082 1.0000
## Aphagous / Pollinator 0.669 0.3720 Inf 1 -0.723 0.9994
## Aphagous / Pursuit Predator 0.424 0.2330 Inf 1 -1.564 0.8656
## Aphagous / Variable 1.396 0.7790 Inf 1 0.598 0.9999
## Decomposer / Fluid Feeder 1.924 0.5680 Inf 1 2.215 0.4462
## Decomposer / Grazer 1.428 0.5350 Inf 1 0.952 0.9947
## Decomposer / Opportunist 0.648 0.1580 Inf 1 -1.775 0.7510
## Decomposer / Parasitoid 1.634 1.4400 Inf 1 0.556 0.9999
## Decomposer / Pollinator 1.008 0.3080 Inf 1 0.025 1.0000
## Decomposer / Pursuit Predator 0.638 0.1870 Inf 1 -1.533 0.8791
## Decomposer / Variable 2.102 0.6520 Inf 1 2.394 0.3292
## Fluid Feeder / Grazer 0.742 0.2970 Inf 1 -0.745 0.9992
## Fluid Feeder / Opportunist 0.337 0.0952 Inf 1 -3.850 0.0046
## Fluid Feeder / Parasitoid 0.849 0.7600 Inf 1 -0.182 1.0000
## Fluid Feeder / Pollinator 0.524 0.1770 Inf 1 -1.919 0.6562
## Fluid Feeder / Pursuit Predator 0.332 0.1080 Inf 1 -3.386 0.0248
## Fluid Feeder / Variable 1.093 0.3720 Inf 1 0.260 1.0000
## Grazer / Opportunist 0.454 0.1650 Inf 1 -2.170 0.4774
## Grazer / Parasitoid 1.144 1.0600 Inf 1 0.146 1.0000
## Grazer / Pollinator 0.706 0.2880 Inf 1 -0.855 0.9976
## Grazer / Pursuit Predator 0.447 0.1780 Inf 1 -2.021 0.5838
## Grazer / Variable 1.472 0.6050 Inf 1 0.940 0.9951
## Opportunist / Parasitoid 2.523 2.2200 Inf 1 1.052 0.9889
## Opportunist / Pollinator 1.556 0.4480 Inf 1 1.533 0.8789
## Opportunist / Pursuit Predator 0.985 0.2710 Inf 1 -0.055 1.0000
## Opportunist / Variable 3.246 0.9660 Inf 1 3.954 0.0031
## Parasitoid / Pollinator 0.617 0.5540 Inf 1 -0.538 0.9999
## Parasitoid / Pursuit Predator 0.390 0.3490 Inf 1 -1.052 0.9890
## Parasitoid / Variable 1.286 1.1600 Inf 1 0.280 1.0000
## Pollinator / Pursuit Predator 0.633 0.2060 Inf 1 -1.403 0.9268
## Pollinator / Variable 2.086 0.7290 Inf 1 2.105 0.5237
## Pursuit Predator / Variable 3.295 1.1200 Inf 1 3.522 0.0156
##
## P value adjustment: tukey method for comparing a family of 10 estimates
## Tests are performed on the log scaleRESULT: Functional group significantly improved prediction of total microplastic abundance after accounting for arthropod dry biomass (likelihood-ratio test: chi^2_9 = 24.46, p = 0.0036), indicating that MP abundance differed among functional groups independently of biomass.
1.2 Figure
figure_6.2 <- ggplot(
pitfall,
aes(
x = functional_group,
y = total_mp,
fill = functional_group
)
) +
geom_boxplot(
outlier.shape = NA,
color = "black",
alpha = 0.85
) +
geom_jitter(
width = 0.2,
color = "black",
alpha = 0.85
) +
scale_y_log10() +
scale_fill_manual(
values = c(
"Ambush Predator" = "#FD5901",
"Aphagous" = "#F78104",
"Decomposer" = "#FAAB36",
"Fluid Feeder" = "#FFD700",
"Grazer" = "#008062",
"Opportunist" = "#008083",
"Parasitoid" = "#005F60",
"Pollinator" = "#003738",
"Pursuit Predator" = "#001F1F",
"Variable" = "#BDBDBD"
)
) +
theme_classic(base_size = 13) +
labs(
x = "Functional Group",
y = "Total Microplastics"
) +
theme(
axis.title = element_text(face = "bold"),
axis.text.x = element_text(
angle = 45,
hjust = 1
),
legend.position = "none"
)
figure_6.2
1.3 Do predaceous groups (ambush predators and pursuit predators) accumulate more MPs than herbivorous groups (grazers and fluid feeders) while accounting for mass?
pitfall <- pitfall %>%
mutate(
feeding_category = case_when(
functional_group %in%
c("Ambush Predator", "Pursuit Predator") ~ "Predator",
functional_group %in%
c("Fluid Feeder", "Grazer") ~ "Herbivore",
TRUE ~ NA_character_
),
feeding_category = factor(feeding_category)
)
pred_herb <- pitfall %>%
filter(
feeding_category %in%
c("Predator", "Herbivore")
)
table(pred_herb$feeding_category)##
## Herbivore Predator
## 21 21model_predherb_null <- glm.nb(
total_mp ~ 1,
data = pred_herb
)
model_predherb_mass <- glm.nb(
total_mp ~ dryweight_mg,
data = pred_herb
)
model_predherb <- glm.nb(
total_mp ~ feeding_category + dryweight_mg,
data = pred_herb
)
anova(
model_predherb_null,
model_predherb_mass,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: total_mp
## Model theta Resid. df 2 x log-lik. Test df LR stat.
## 1 1 0.8893699 41 -511.5131
## 2 dryweight_mg 1.3540203 40 -489.4265 1 vs 2 1 22.08661
## Pr(Chi)
## 1
## 2 2.606215e-06anova(
model_predherb_mass,
model_predherb,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: total_mp
## Model theta Resid. df 2 x log-lik. Test
## 1 dryweight_mg 1.354020 40 -489.4265
## 2 feeding_category + dryweight_mg 1.631966 39 -480.1399 1 vs 2
## df LR stat. Pr(Chi)
## 1
## 2 1 9.286546 0.00230843summary(model_predherb)##
## Call:
## glm.nb(formula = total_mp ~ feeding_category + dryweight_mg,
## data = pred_herb, init.theta = 1.631966102, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.177277 0.172944 24.154 < 2e-16 ***
## feeding_categoryPredator 0.846160 0.261921 3.231 0.001235 **
## dryweight_mg 0.007131 0.002097 3.401 0.000671 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.632) family taken to be 1)
##
## Null deviance: 89.104 on 41 degrees of freedom
## Residual deviance: 45.822 on 39 degrees of freedom
## AIC: 488.14
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.632
## Std. Err.: 0.332
##
## 2 x log-likelihood: -480.140MPs per unit biomass
pred_herb <- pred_herb %>%
mutate(
mp_per_mg = total_mp / dryweight_mg
)
model_mp_per_mg <- lm(
mp_per_mg ~ feeding_category,
data = pred_herb
)
summary(model_mp_per_mg)##
## Call:
## lm(formula = mp_per_mg ~ feeding_category, data = pred_herb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -202.27 -117.78 -89.11 80.31 420.31
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 219.69 43.59 5.040 1.04e-05 ***
## feeding_categoryPredator -113.30 61.65 -1.838 0.0735 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 199.8 on 40 degrees of freedom
## Multiple R-squared: 0.07786, Adjusted R-squared: 0.05481
## F-statistic: 3.377 on 1 and 40 DF, p-value: 0.07353Note: predators and herbivores differed significantly in total MP abundance after accounting for dry biomass
Follow-up: do predaceous groups contain more MPs per gram of weight than herbivorous groups?
1.3 Figure
model_predherb_offset_null <- glm.nb(
total_mp ~ offset(log(dryweight_mg)),
data = pred_herb
)
model_predherb_offset <- glm.nb(
total_mp ~ feeding_category +
offset(log(dryweight_mg)),
data = pred_herb
)
anova(
model_predherb_offset_null,
model_predherb_offset,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: total_mp
## Model theta Resid. df
## 1 offset(log(dryweight_mg)) 0.5517053 41
## 2 feeding_category + offset(log(dryweight_mg)) 0.5833399 40
## 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 -539.0558
## 2 -535.8685 1 vs 2 1 3.187299 0.07421252RESULT: Feeding guild significantly explained variation in total MP abundance after accounting for dry biomass (chi^2 = 9.29, p = 0.0023). However, predator and herbivore groups did not differ significantly in MP abundance per unit biomass when dry mass was incorporated as an offset (chi^2 = 3.19, p = 0.074)
Pitfall descriptive summaries
total_pitfall_mp <- pitfall %>%
summarise(
total_mp = sum(total_mp, na.rm = TRUE)
)
total_pitfall_mp## total_mp
## 1 18741mp_by_functional_group <- pitfall %>%
group_by(functional_group) %>%
summarise(
total_mp = sum(total_mp, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(desc(total_mp))
mp_by_functional_group## # A tibble: 10 × 2
## functional_group total_mp
## <fct> <dbl>
## 1 Opportunist 6425
## 2 Pursuit Predator 4293
## 3 Decomposer 2609
## 4 Pollinator 1902
## 5 Ambush Predator 1148
## 6 Fluid Feeder 823
## 7 Variable 684
## 8 Grazer 563
## 9 Aphagous 225
## 10 Parasitoid 69n_distinct(pitfall$family)## [1] 602. Benthic Macroinvertebrates
Import and clean data
benthics <- read.csv(
file.path(data_dir, "arthropods.csv")
) %>%
clean_names() %>%
dplyr::select(where(~ !all(is.na(.))))benthics <- benthics %>%
mutate(
sample = factor(sample),
site = as.character(site),
family = as.character(family),
season = factor(season),
dryweight_mg = as.numeric(t_dryweight) * 1000,
count = as.numeric(count),
mp_total = as.numeric(mp_total),
fragment = as.numeric(fragment),
fiber = as.numeric(fiber),
film = as.numeric(film),
microbead = as.numeric(microbead)
)
head(benthics)## season site sample family count t_dryweight mp_total fragment fiber
## 1 S25 EM P Solenoceridae 1 0.0547 542 540 1
## 2 S25 EM P Aristeidae 2 0.0130 632 619 13
## 3 F25 CB B Penaeidae 21 0.4682 426 421 5
## 4 F25 CB B Aristeidae 5 0.0030 478 466 12
## 5 F25 TB B Aristeidae 7 0.4731 220 217 3
## 6 F24 TB B Penaeidae 10 0.5237 5445 4962 147
## microbead film initials dryweight_mg
## 1 1 0 GEM 54.7
## 2 0 0 13.0
## 3 0 0 KBA 468.2
## 4 0 0 GAT 3.0
## 5 0 0 GAT 473.1
## 6 331 5 GEM, FA, KAW 523.7str(benthics)## 'data.frame': 16 obs. of 13 variables:
## $ season : Factor w/ 3 levels "F24","F25","S25": 3 3 2 2 2 1 1 2 2 1 ...
## $ site : chr "EM" "EM" "CB" "CB" ...
## $ sample : Factor w/ 2 levels "B","P": 2 2 1 1 1 1 1 1 1 1 ...
## $ family : chr "Solenoceridae" "Aristeidae" "Penaeidae" "Aristeidae" ...
## $ count : num 1 2 21 5 7 10 2 16 1 1 ...
## $ t_dryweight : num 0.0547 0.013 0.4682 0.003 0.4731 ...
## $ mp_total : num 542 632 426 478 220 ...
## $ fragment : num 540 619 421 466 217 ...
## $ fiber : num 1 13 5 12 3 147 107 3 8 15 ...
## $ microbead : num 1 0 0 0 0 331 0 0 0 0 ...
## $ film : num 0 0 0 0 0 5 0 0 0 0 ...
## $ initials : chr "GEM" "" "KBA" "GAT" ...
## $ dryweight_mg: num 54.7 13 468.2 3 473.1 ...summary(benthics$dryweight_mg)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05 3.00 20.20 122.82 121.53 523.702.1 Does total MP abundance increase with specimen dry mass among benthic taxa?
cor.test(
benthics$dryweight_mg,
benthics$mp_total,
method = "spearman",
exact = FALSE
)##
## Spearman's rank correlation rho
##
## data: benthics$dryweight_mg and benthics$mp_total
## S = 371.77, p-value = 0.07785
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4532746model_benths_null <- glm.nb(
mp_total ~ 1,
data = benthics
)
model_benths_mass <- glm.nb(
mp_total ~ dryweight_mg,
data = benthics
)
summary(model_benths_mass)##
## Call:
## glm.nb(formula = mp_total ~ dryweight_mg, data = benthics, init.theta = 1.68393391,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.351505 0.231828 27.398 <2e-16 ***
## dryweight_mg 0.002664 0.001046 2.547 0.0109 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.6839) family taken to be 1)
##
## Null deviance: 25.559 on 15 degrees of freedom
## Residual deviance: 17.516 on 14 degrees of freedom
## AIC: 249.5
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.684
## Std. Err.: 0.548
##
## 2 x log-likelihood: -243.503anova(
model_benths_null,
model_benths_mass,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: mp_total
## Model theta Resid. df 2 x log-lik. Test df LR stat.
## 1 1 1.191541 15 -250.2263
## 2 dryweight_mg 1.683934 14 -243.5031 1 vs 2 1 6.723258
## Pr(Chi)
## 1
## 2 0.009516366Count sensitivity analyses
model_benths_count <- glm.nb(
mp_total ~ count,
data = benthics
)
model_benths_both <- glm.nb(
mp_total ~ dryweight_mg + count,
data = benthics
)
summary(model_benths_count)##
## Call:
## glm.nb(formula = mp_total ~ count, data = benthics, init.theta = 1.193159027,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.848665 0.252984 27.071 <2e-16 ***
## count -0.001928 0.009772 -0.197 0.844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.1932) family taken to be 1)
##
## Null deviance: 18.123 on 15 degrees of freedom
## Residual deviance: 18.095 on 14 degrees of freedom
## AIC: 256.2
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.193
## Std. Err.: 0.378
##
## 2 x log-likelihood: -250.199summary(model_benths_both)##
## Call:
## glm.nb(formula = mp_total ~ dryweight_mg + count, data = benthics,
## init.theta = 1.686946849, link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.367401 0.247032 25.776 <2e-16 ***
## dryweight_mg 0.002662 0.001046 2.546 0.0109 *
## count -0.001486 0.008230 -0.181 0.8567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.6869) family taken to be 1)
##
## Null deviance: 25.605 on 15 degrees of freedom
## Residual deviance: 17.514 on 13 degrees of freedom
## AIC: 251.47
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.687
## Std. Err.: 0.549
##
## 2 x log-likelihood: -243.469anova(
model_benths_null,
model_benths_count,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: mp_total
## Model theta Resid. df 2 x log-lik. Test df LR stat. Pr(Chi)
## 1 1 1.191541 15 -250.2263
## 2 count 1.193159 14 -250.1992 1 vs 2 1 0.02714271 0.8691404anova(
model_benths_mass,
model_benths_both,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: mp_total
## Model theta Resid. df 2 x log-lik. Test df
## 1 dryweight_mg 1.683934 14 -243.5031
## 2 dryweight_mg + count 1.686947 13 -243.4693 1 vs 2 1
## LR stat. Pr(Chi)
## 1
## 2 0.03376103 0.8542161RESULTS: Total MP abundance increased significantly with benthic taxa dry biomass. A negative binomial model indicated a positive biomass effect (beta=0.00266, p=0.011), and inclusion of biomass significantly improved model fit relative to an intercept-only model (likelihood-ratio test: (LR=6.69, df=1, p=0.0097). This corresponded to an approximately 30% increase in expected MP abundance per 100-mg increase in pooled dry biomass.
2.2 Does MP abundance and type composition differ between taxonomic group?
benthic_composition <- benthics %>%
summarise(
fragments = sum(fragment, na.rm = TRUE),
fibers = sum(fiber, na.rm = TRUE),
films = sum(film, na.rm = TRUE),
microbeads = sum(microbead, na.rm = TRUE),
total_mp = sum(mp_total, na.rm = TRUE)
)
benthic_composition## fragments fibers films microbeads total_mp
## 1 14110 329 5 332 14776benthic_composition_percent <- benthics %>%
summarise(
Fragments = sum(fragment, na.rm = TRUE),
Fibers = sum(fiber, na.rm = TRUE),
Films = sum(film, na.rm = TRUE),
Microbeads = sum(microbead, na.rm = TRUE)
) %>%
pivot_longer(
cols = everything(),
names_to = "mp_type",
values_to = "count"
) %>%
mutate(
percent = 100 * count / sum(count)
)
benthic_composition_percent## # A tibble: 4 × 3
## mp_type count percent
## <chr> <dbl> <dbl>
## 1 Fragments 14110 95.5
## 2 Fibers 329 2.23
## 3 Films 5 0.0338
## 4 Microbeads 332 2.25MP abundance by taxon
benthic_taxon_summary <- benthics %>%
group_by(family) %>%
summarise(
n = n(),
mean_mp = mean(mp_total, na.rm = TRUE),
median_mp = median(mp_total, na.rm = TRUE),
min_mp = min(mp_total, na.rm = TRUE),
max_mp = max(mp_total, na.rm = TRUE),
total_mp = sum(mp_total, na.rm = TRUE),
.groups = "drop"
)
benthic_taxon_summary## # A tibble: 5 × 7
## family n mean_mp median_mp min_mp max_mp total_mp
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Aristeidae 5 555. 543 220 904 2777
## 2 Gammaridae 2 406. 406. 190 621 811
## 3 Littorinidae 1 723 723 723 723 723
## 4 Penaeidae 3 2174. 650 426 5445 6521
## 5 Solenoceridae 5 789. 542 117 2317 39443. Gelatinous Zooplankton
Import and clean data
cnidaria <- read.csv(
file.path(data_dir, "cnidaria.csv")
) %>%
clean_names() %>%
dplyr::select(where(~ !all(is.na(.))))cnidaria <- cnidaria %>%
mutate(
sample_type = factor(sample_type),
season = factor(season),
dryweight_mg = as.numeric(t_dryweight) * 1000,
mp_total = as.numeric(mp_total),
fragment = as.numeric(fragment),
fiber = as.numeric(fiber),
film = as.numeric(film),
microbead = as.numeric(microbead)
)
head(cnidaria)## sample_type season site taxon t_dryweight mp_total fragment fiber
## 1 independent F24 KB Stomolophidae 2.27800 265 263 2
## 2 pelagic F25 EM Ctenophora 0.00005 1462 1460 2
## 3 independent F24 KB Stomolophidae 3.74500 735 734 0
## 4 pelagic F25 LB Ctenophora 3.10070 289 287 2
## 5 independent F24 TB Stomolophidae 1.41390 517 513 3
## 6 independent F24 KB Stomolophidae 4.36990 690 686 4
## film microbead initials dryweight_mg
## 1 0 0 KBA 2278.00
## 2 0 0 KBA 0.05
## 3 0 1 KBA 3745.00
## 4 0 0 KBA 3100.70
## 5 1 0 KBA 1413.90
## 6 0 0 KBA 4369.90str(cnidaria)## 'data.frame': 13 obs. of 12 variables:
## $ sample_type : Factor w/ 3 levels "benthic","independent",..: 2 3 2 3 2 2 3 3 1 1 ...
## $ season : Factor w/ 2 levels "F24","F25": 1 2 1 2 1 1 2 2 2 2 ...
## $ site : chr "KB" "EM" "KB" "LB" ...
## $ taxon : chr "Stomolophidae" "Ctenophora" "Stomolophidae" "Ctenophora" ...
## $ t_dryweight : num 2.278 0.00005 3.745 3.1007 1.4139 ...
## $ mp_total : num 265 1462 735 289 517 ...
## $ fragment : num 263 1460 734 287 513 ...
## $ fiber : num 2 2 0 2 3 4 4 6 0 6 ...
## $ film : num 0 0 0 0 1 0 0 0 1 0 ...
## $ microbead : num 0 0 1 0 0 0 0 0 0 0 ...
## $ initials : chr "KBA" "KBA" "KBA" "KBA" ...
## $ dryweight_mg: num 2278 0.05 3745 3100.7 1413.9 ...summary(cnidaria$dryweight_mg)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.05 0.05 0.10 1382.85 3069.10 4369.90Separate standardized nektonic samples and opportunistic Stomolophus samples
stomolophus <- cnidaria %>%
filter(
taxon == "Stomolophidae",
sample_type == "independent"
)
gel.main <- cnidaria %>%
filter(
!(
taxon == "Stomolophidae" &
sample_type == "independent"
)
)
nrow(cnidaria)## [1] 13nrow(gel.main)## [1] 9nrow(stomolophus)## [1] 4table(gel.main$taxon)##
## Ctenophora Ulmaridae
## 8 1table(gel.main$sample_type)##
## benthic independent pelagic
## 2 0 73.1 Does total MP abundance increase with specimen dry mass in the jelly samples?
cor.test(
gel.main$dryweight_mg,
gel.main$mp_total,
method = "spearman",
exact = FALSE
)##
## Spearman's rank correlation rho
##
## data: gel.main$dryweight_mg and gel.main$mp_total
## S = 124.75, p-value = 0.9194
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.0396059model_gel_null <- glm.nb(
mp_total ~ 1,
data = gel.main
)
model_gel_mass <- glm.nb(
mp_total ~ dryweight_mg,
data = gel.main
)
summary(model_gel_mass)##
## Call:
## glm.nb(formula = mp_total ~ dryweight_mg, data = gel.main, init.theta = 1.836831793,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.900e+00 2.791e-01 24.718 <2e-16 ***
## dryweight_mg 6.172e-05 1.919e-04 0.322 0.748
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8368) family taken to be 1)
##
## Null deviance: 9.8973 on 8 degrees of freedom
## Residual deviance: 9.7908 on 7 degrees of freedom
## AIC: 147.31
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.837
## Std. Err.: 0.802
##
## 2 x log-likelihood: -141.313anova(
model_gel_null,
model_gel_mass,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: mp_total
## Model theta Resid. df 2 x log-lik. Test df LR stat.
## 1 1 1.818436 8 -141.4188
## 2 dryweight_mg 1.836832 7 -141.3128 1 vs 2 1 0.1059457
## Pr(Chi)
## 1
## 2 0.7448079Log-mass sensitivity analysis
gel.main <- gel.main %>%
mutate(
log_dryweight = log10(dryweight_mg)
)
model_gel_logmass <- glm.nb(
mp_total ~ log_dryweight,
data = gel.main
)
summary(model_gel_logmass)##
## Call:
## glm.nb(formula = mp_total ~ log_dryweight, data = gel.main, init.theta = 1.834398384,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.9501 0.2476 28.067 <2e-16 ***
## log_dryweight 0.0370 0.1247 0.297 0.767
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(1.8344) family taken to be 1)
##
## Null deviance: 9.8842 on 8 degrees of freedom
## Residual deviance: 9.7918 on 7 degrees of freedom
## AIC: 147.33
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 1.834
## Std. Err.: 0.800
##
## 2 x log-likelihood: -141.327anova(
model_gel_null,
model_gel_logmass,
test = "Chisq"
)## Likelihood ratio tests of Negative Binomial Models
##
## Response: mp_total
## Model theta Resid. df 2 x log-lik. Test df LR stat.
## 1 1 1.818436 8 -141.4188
## 2 log_dryweight 1.834398 7 -141.3268 1 vs 2 1 0.09199762
## Pr(Chi)
## 1
## 2 0.7616528RESULT: Total microplastic abundance was not associated with specimen dry biomass in gelatinous samples. Neither raw dry weight (beta=-4.72^{-5}, p=0.693) nor log-transformed dry weight (beta=-0.0526, p=0.521) significantly predicted MP abundance, and neither model improved fit relative to an intercept-only negative-binomial model.
Note: Biomass appears to be an informative predictor of MP abundance in benthic taxa, but not in gelatinous taxa, suggesting that dry biomass may not represent MP exposure or retention equally across different taxa.
8.2 Does MP abundance and species composition differ between taxonomic group?
gel_composition <- gel.main %>%
summarise(
fragments = sum(fragment, na.rm = TRUE),
fibers = sum(fiber, na.rm = TRUE),
films = sum(film, na.rm = TRUE),
microbeads = sum(microbead, na.rm = TRUE),
total_mp = sum(mp_total, na.rm = TRUE)
)
gel_composition## fragments fibers films microbeads total_mp
## 1 9305 38 2 0 9345gel_composition_percent <- gel.main %>%
summarise(
Fragments = sum(fragment, na.rm = TRUE),
Fibers = sum(fiber, na.rm = TRUE),
Films = sum(film, na.rm = TRUE),
Microbeads = sum(microbead, na.rm = TRUE)
) %>%
pivot_longer(
cols = everything(),
names_to = "mp_type",
values_to = "count"
) %>%
mutate(
percent = 100 * count / sum(count)
)
gel_composition_percent## # A tibble: 4 × 3
## mp_type count percent
## <chr> <dbl> <dbl>
## 1 Fragments 9305 99.6
## 2 Fibers 38 0.407
## 3 Films 2 0.0214
## 4 Microbeads 0 0Standardized gelatinous taxon summary
gel_taxon_summary <- gel.main %>%
group_by(taxon) %>%
summarise(
n = n(),
total_mp = sum(mp_total, na.rm = TRUE),
mean_mp = mean(mp_total, na.rm = TRUE),
sd_mp = sd(mp_total, na.rm = TRUE),
median_mp = median(mp_total, na.rm = TRUE),
min_mp = min(mp_total, na.rm = TRUE),
max_mp = max(mp_total, na.rm = TRUE),
.groups = "drop"
)
gel_taxon_summary## # A tibble: 2 × 8
## taxon n total_mp mean_mp sd_mp median_mp min_mp max_mp
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Ctenophora 8 7226 903. 677. 687 206 2031
## 2 Ulmaridae 1 2119 2119 NA 2119 2119 2119RESULT: Jelly specimens contained 11,552 MP particles, of which 99.56% were fragments. Fragment dominance was consistent among Ctenophora, Stomolophidae, and Ulmaridae. Ctenophore samples contained a mean of 903 MPs per specimen (median = 687, range = 206–2,031), while Stomolophidae contained a mean of 552 MPs (median = 604, range = 265–735). The single Ulmaridae specimen contained 2,119 MPs.
Note: Four Stomolophidae specimens encountered during field sampling were collected opportunistically and were not associated with standardized transect sampling. These specimens were processed to provide descriptive estimates of microplastic burden and were excluded from inferential comparisons requiring standardized sampling effort. Opportunistically collected Stomolophidae specimens (n=4) contained 265–735 MPs per individual (mean = 552; median = 604), with fragments accounting for approximately 99.5% of recovered particles.
Opportunistic Stomolophus descriptive summary
stomolophus_summary <- stomolophus %>%
summarise(
n = n(),
total_mp = sum(mp_total, na.rm = TRUE),
mean_mp = mean(mp_total, na.rm = TRUE),
sd_mp = sd(mp_total, na.rm = TRUE),
median_mp = median(mp_total, na.rm = TRUE),
min_mp = min(mp_total, na.rm = TRUE),
max_mp = max(mp_total, na.rm = TRUE),
mean_mass_mg = mean(dryweight_mg, na.rm = TRUE),
sd_mass_mg = sd(dryweight_mg, na.rm = TRUE),
min_mass_mg = min(dryweight_mg, na.rm = TRUE),
max_mass_mg = max(dryweight_mg, na.rm = TRUE)
)
stomolophus_summary## n total_mp mean_mp sd_mp median_mp min_mp max_mp mean_mass_mg sd_mass_mg
## 1 4 2207 551.75 213.0155 603.5 265 735 2951.7 1348.989
## min_mass_mg max_mass_mg
## 1 1413.9 4369.94. Combined Benthic and Gelatinous Summary
Only standardized gelatinous samples (gel.main) are
included here.
mp_organisms <- bind_rows(
benthics %>%
transmute(
organism_group = "Benthic Macroinvertebrates",
taxon = family,
mp_total = mp_total
),
gel.main %>%
transmute(
organism_group = "Gelatinous Zooplankton",
taxon = taxon,
mp_total = mp_total
)
)Summary table
mp_table <- mp_organisms %>%
group_by(
organism_group,
taxon
) %>%
summarise(
n = n(),
total_mp = sum(mp_total, na.rm = TRUE),
mean_mp = mean(mp_total, na.rm = TRUE),
sd_mp = sd(mp_total, na.rm = TRUE),
median_mp = median(mp_total, na.rm = TRUE),
min_mp = min(mp_total, na.rm = TRUE),
max_mp = max(mp_total, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
`Mean ± SD` = if_else(
is.na(sd_mp),
sprintf("%.1f ± NA", mean_mp),
sprintf("%.1f ± %.1f", mean_mp, sd_mp)
),
Range = paste0(
min_mp,
"–",
max_mp
)
) %>%
dplyr::select(
organism_group,
taxon,
n,
total_mp,
`Mean ± SD`,
median_mp,
Range
)
mp_table## # A tibble: 7 × 7
## organism_group taxon n total_mp `Mean ± SD` median_mp Range
## <chr> <chr> <int> <dbl> <chr> <dbl> <chr>
## 1 Benthic Macroinvertebrates Aristei… 5 2777 555.4 ± 24… 543 220–…
## 2 Benthic Macroinvertebrates Gammari… 2 811 405.5 ± 30… 406. 190–…
## 3 Benthic Macroinvertebrates Littori… 1 723 723.0 ± NA 723 723–…
## 4 Benthic Macroinvertebrates Penaeid… 3 6521 2173.7 ± 2… 650 426–…
## 5 Benthic Macroinvertebrates Solenoc… 5 3944 788.8 ± 87… 542 117–…
## 6 Gelatinous Zooplankton Ctenoph… 8 7226 903.2 ± 67… 687 206–…
## 7 Gelatinous Zooplankton Ulmarid… 1 2119 2119.0 ± NA 2119 2119…Combined organism figure
mp_plot_data <- mp_organisms %>%
mutate(
organism_group = recode(
organism_group,
"Benthic Macroinvertebrates" = "Benthics",
"Gelatinous Zooplankton" = "Jellies"
)
)
figure_mp_organisms <- ggplot(
mp_plot_data,
aes(
x = mp_total,
y = reorder(taxon, mp_total, FUN = median),
color = organism_group
)
) +
geom_jitter(
height = 0.12,
width = 0,
alpha = 0.55,
size = 3
) +
stat_summary(
fun = mean,
geom = "point",
shape = 18,
size = 5
) +
facet_grid(
organism_group ~ .,
scales = "free_y",
space = "free_y"
) +
scale_color_manual(
values = c(
"Benthics" = "deepskyblue1",
"Jellies" = "hotpink"
)
) +
scale_x_log10() +
labs(
x = "MP Abundance Per Sample",
y = NULL
) +
theme_classic(base_size = 16) +
theme(
strip.text = element_text(
face = "bold",
size = 15
),
axis.title = element_text(face = "bold"),
axis.text.y = element_text(
face = "italic",
size = 13
),
legend.position = "none",
panel.spacing = grid::unit(1, "lines")
)
figure_mp_organisms
Session Information
sessionInfo()## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Chicago
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] gt_1.3.0 janitor_2.2.1 emmeans_2.0.2 MASS_7.3-60
## [5] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [9] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
## [13] ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.4 sass_0.4.9 generics_0.1.3 xml2_1.3.8
## [5] stringi_1.8.4 lattice_0.21-9 hms_1.1.3 digest_0.6.37
## [9] magrittr_2.0.3 evaluate_1.0.3 grid_4.3.2 timechange_0.3.0
## [13] estimability_1.5.1 mvtnorm_1.3-3 fastmap_1.2.0 jsonlite_1.8.9
## [17] scales_1.3.0 jquerylib_0.1.4 cli_3.6.3 rlang_1.1.5
## [21] munsell_0.5.1 withr_3.0.2 cachem_1.1.0 yaml_2.3.10
## [25] tools_4.3.2 tzdb_0.4.0 coda_0.19-4.1 colorspace_2.1-1
## [29] vctrs_0.6.5 R6_2.5.1 lifecycle_1.0.4 snakecase_0.11.1
## [33] fs_1.6.5 pkgconfig_2.0.3 pillar_1.10.1 bslib_0.9.0
## [37] gtable_0.3.6 glue_1.8.0 xfun_0.52 tidyselect_1.2.1
## [41] rstudioapi_0.17.1 knitr_1.49 farver_2.1.2 htmltools_0.5.8.1
## [45] rmarkdown_2.30 compiler_4.3.2Citations
Firke, S. 2024. janitor: Simple Tools for Examining and Cleaning Dirty Data. https://CRAN.R-project.org/package=janitor.
Iannone, R., J. Cheng, B. Schloerke, S. Haughton, E. Hughes, A. Lauer, R. François, J. Seo, K. Brevoort, and O. Roy. 2026. gt: Easily Create Presentation-Ready Display Tables. https://CRAN.R-project.org/package=gt.
Lenth, R. V., and J. Piaskowski. 2026. emmeans: Estimated Marginal Means, aka Least-Squares Means. https://CRAN.R-project.org/package=emmeans.
R Core Team. 2023. R: A Language and Environment for Statistical Computing. R Foundation for Spatistical Computing, Vienna, Austria. https://www.R-project.org/.
Venables, W. N., and B. D. Ripley. 2002. Modern Applied Statistics with S. Fourth. Springer, New York.
Wickham, H. 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York.
Wickham, H., M. Averick, J. Bryan, W. Chang, L. D. McGowan, R. François, G. Grolemund, A. Hayes, L. Henry, J. Hester, M. Kuhn, T. L. Pedersen, E. Miller, S. M. Bache, K. Müller, J. Ooms, D. Robinson, D. P. Seidel, V. Spinu, K. Takahashi, D. Vaughan, C. Wilke, K. Woo, and H. Yutani. 2019. Welcome to the tidyverse. Journal of Open Source Software 4:1686.
Wickham, H., R. François, L. Henry, K. Müller, and D. Vaughan. 2023. dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, H., D. Vaughan, and M. Girlich. 2024. tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.