Respiration larvae + sediments
Ana M. Palacio & Xaymara Serrano
2022-02-02
General project set-up
# Load all libraries and sources required to run the script
library(tidyverse)
library(ggthemes)
library(ggplot2)
library(lme4)
library(lmerTest)
#library(multcomp)
#library(multcompView)
#library(emmeans)
#library(effects)
ggthe_bw<-theme_bw()+
theme(panel.grid= element_blank(),
legend.box.background = element_rect(),
panel.background =element_rect(fill = NA, color = "black")
)1. Data
# Data
data<-read.csv("Data/Respiration_larvae.csv", header = TRUE)
summary(data)## Species Date Timepoint Treatment.code Treatment
## Ofav:75 08/23/17:38 Exposure:38 NS:15 Control :15
## 08/29/17:37 Recovery:37 PH:15 High Port:15
## PL:15 High Reef:15
## RH:15 Low Port :15
## RL:15 Low Reef :15
##
##
## Replicate.vial Plate.. Well nmol.min X..larvae
## Min. : 1.000 Min. :1.00 A1 : 4 Min. :0.001305 Min. :1.0
## 1st Qu.: 3.500 1st Qu.:1.00 A3 : 4 1st Qu.:0.003019 1st Qu.:1.0
## Median : 6.000 Median :2.00 A4 : 4 Median :0.004709 Median :3.0
## Mean : 5.813 Mean :2.48 A5 : 4 Mean :0.004569 Mean :2.8
## 3rd Qu.: 8.000 3rd Qu.:3.50 A6 : 4 3rd Qu.:0.005861 3rd Qu.:4.0
## Max. :10.000 Max. :4.00 B1 : 4 Max. :0.009487 Max. :7.0
## (Other):51 NA's :2 NA's :10
## nmol.larva.min
## Min. :0.000550
## 1st Qu.:0.001189
## Median :0.001607
## Mean :0.002100
## 3rd Qu.:0.002409
## Max. :0.007717
## 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 152. Exploratory plots
Plot-well
Well<- ggplot(data, aes (Treatment, nmol.min)) +
geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",
geom = "errorbar", width = 0.2, color="blue")+
stat_summary(fun=mean, geom="point", color="blue") +
#geom_point(shape=21)+
geom_jitter(alpha=0.5, shape=21)+
theme(legend.position = "bottom")+
scale_y_continuous(limits = c(0, 0.01),
# breaks = seq(0,1,0.2),
# expand = c(0.01, 0.01),
name=("O2 [nmol/min]")) +
ggthe_bw +
facet_grid(~Timepoint)
Well
LogWell<- ggplot(data, aes (Treatment, log(nmol.min))) +
geom_boxplot ()+
stat_summary(fun.data = "mean_cl_boot",
geom = "errorbar", width = 0.2, color="blue")+
stat_summary(fun=mean, geom="point", color="blue") +
#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]")) +
ggthe_bw +
facet_grid(~Timepoint)
LogWell
Plot-larvae
Larvae<- ggplot(data, aes (Treatment, 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=("O2 [nmol/min/larva]")) +
ggthe_bw +
facet_grid(~Timepoint)
Larvae
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]")) +
ggthe_bw +
facet_grid(~Timepoint)
LogLarvae
3. Possible GML and LMs
## model 1: LMER for exposure well
fit1<-lmerTest::lmer(nmol.min ~Treatment + (1|Replicate.vial),
data=data[data$Timepoint=="Exposure",])
anova(fit1) ranova(fit1) summary(fit1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nmol.min ~ Treatment + (1 | Replicate.vial)
## Data: data[data$Timepoint == "Exposure", ]
##
## REML criterion at convergence: -308.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.24879 -0.78775 0.01597 0.59450 1.70641
##
## Random effects:
## Groups Name Variance Std.Dev.
## Replicate.vial (Intercept) 4.278e-07 0.000654
## Residual 2.442e-06 0.001563
## Number of obs: 37, groups: Replicate.vial, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.0061999 0.0006378 30.4985473 9.721 7.43e-11 ***
## TreatmentLow Reef -0.0009618 0.0008418 26.3373595 -1.143 0.2635
## TreatmentHigh Reef -0.0001514 0.0008115 25.6907409 -0.187 0.8534
## TreatmentLow Port -0.0006956 0.0008418 26.3373595 -0.826 0.4160
## TreatmentHigh Port -0.0016786 0.0008115 25.6907409 -2.069 0.0488 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.659
## TrtmntHghRf -0.683 0.518
## TrtmntLwPrt -0.659 0.508 0.518
## TrtmntHghPr -0.683 0.518 0.537 0.518 plot(fit1)
#step(fit1)
fit1<-lmerTest::lmer(log10(nmol.min) ~Treatment + (1|Replicate.vial),
data=data[data$Timepoint=="Exposure",])
anova(fit1) ranova(fit1) summary(fit1)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(nmol.min) ~ Treatment + (1 | Replicate.vial)
## Data: data[data$Timepoint == "Exposure", ]
##
## REML criterion at convergence: -22.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.67728 -0.69997 0.06863 0.60093 1.60774
##
## Random effects:
## Groups Name Variance Std.Dev.
## Replicate.vial (Intercept) 0.004018 0.06339
## Residual 0.018204 0.13492
## Number of obs: 37, groups: Replicate.vial, 8
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.21542 0.05605 29.80823 -39.528 <2e-16 ***
## TreatmentLow Reef -0.09796 0.07278 26.29701 -1.346 0.1898
## TreatmentHigh Reef -0.01485 0.07010 25.68829 -0.212 0.8339
## TreatmentLow Port -0.06940 0.07278 26.29701 -0.954 0.3490
## TreatmentHigh Port -0.15276 0.07010 25.68829 -2.179 0.0387 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.649
## TrtmntHghRf -0.672 0.519
## TrtmntLwPrt -0.649 0.509 0.519
## TrtmntHghPr -0.672 0.519 0.537 0.519 plot(fit1)
#step(fit1)
## model 2: LMER for exposure larva
fit2<-lmer(nmol.larva.min ~Treatment + (1|Replicate.vial) + (1|Plate..) + (1|Well),
data=data[data$Timepoint=="Exposure",])
fit2<-lmer(nmol.larva.min ~Treatment + (1|Plate..),
data=data[data$Timepoint=="Exposure",])
fit2<-lmer(log(nmol.larva.min) ~Treatment + (1|Plate..),
data=data[data$Timepoint=="Exposure",])
anova(fit2) ranova(fit2) summary(fit2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(nmol.larva.min) ~ Treatment + (1 | Plate..)
## Data: data[data$Timepoint == "Exposure", ]
##
## REML criterion at convergence: 38.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.84658 -0.44652 -0.05528 0.41667 1.99608
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plate.. (Intercept) 0.2472 0.4972
## Residual 0.1624 0.4030
## Number of obs: 31, groups: Plate.., 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -5.8378 0.3833 1.2915 -15.231 0.0205 *
## TreatmentLow Reef -0.1665 0.2539 25.0240 -0.656 0.5179
## TreatmentHigh Reef -0.1747 0.2154 25.0000 -0.811 0.4250
## TreatmentLow Port -0.2865 0.2164 25.0212 -1.324 0.1975
## TreatmentHigh Port -0.4820 0.2246 25.0088 -2.146 0.0418 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.236
## TrtmntHghRf -0.281 0.424
## TrtmntLwPrt -0.282 0.412 0.498
## TrtmntHghPr -0.268 0.413 0.479 0.471 plot(fit2)
step(fit2)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 -19.336 52.672
## (1 | Plate..) 0 6 -25.739 63.478 12.806 1 0.0003455 ***
## ---
## 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 1 0.807 0.20175 4 25.023 1.2425 0.3186
##
## Model found:
## log(nmol.larva.min) ~ (1 | Plate..)## model 3: LMER for recovery well
# fit3<-lmer(nmol.min ~Treatment + (1|Replicate.vial) + (1|Plate..) + (1|Well),
# data=data[data$Timepoint=="Recovery",])
fit3<-lmer(nmol.min ~Treatment + (1|Well),
data=data[data$Timepoint=="Recovery",])
anova(fit3) ranova(fit3) summary(fit3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: nmol.min ~ Treatment + (1 | Well)
## Data: data[data$Timepoint == "Recovery", ]
##
## REML criterion at convergence: -323.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.16815 -0.62384 -0.04064 0.56904 1.33203
##
## Random effects:
## Groups Name Variance Std.Dev.
## Well (Intercept) 9.780e-07 0.0009889
## Residual 5.597e-07 0.0007482
## Number of obs: 36, groups: Well, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.0040769 0.0004331 30.7821139 9.414 1.42e-10 ***
## TreatmentLow Reef -0.0001311 0.0004537 15.1010849 -0.289 0.7765
## TreatmentHigh Reef -0.0010641 0.0005665 26.2639863 -1.878 0.0715 .
## TreatmentLow Port -0.0003695 0.0005648 28.9232131 -0.654 0.5181
## TreatmentHigh Port -0.0009104 0.0006139 29.6098025 -1.483 0.1486
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.615
## TrtmntHghRf -0.696 0.554
## TrtmntLwPrt -0.719 0.530 0.676
## TrtmntHghPr -0.678 0.467 0.599 0.708 plot(fit3)
step(fit3)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 161.59 -309.19
## (1 | Well) 0 6 159.09 -306.17 5.0154 1 0.02512 *
## ---
## 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 1 3.3137e-06 8.2842e-07 4 15.506 1.48 0.2561
##
## Model found:
## nmol.min ~ (1 | Well) fit3<-lmer(log(nmol.min) ~Treatment + (1|Well),
data=data[data$Timepoint=="Recovery",])
anova(fit3) ranova(fit3) summary(fit3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(nmol.min) ~ Treatment + (1 | Well)
## Data: data[data$Timepoint == "Recovery", ]
##
## REML criterion at convergence: 30.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.52057 -0.47973 0.03375 0.55573 1.34214
##
## Random effects:
## Groups Name Variance Std.Dev.
## Well (Intercept) 0.08695 0.2949
## Residual 0.05195 0.2279
## Number of obs: 36, groups: Well, 21
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -5.54535 0.13059 30.78001 -42.463 <2e-16 ***
## TreatmentLow Reef -0.07873 0.13781 15.34690 -0.571 0.576
## TreatmentHigh Reef -0.28176 0.17139 26.53973 -1.644 0.112
## TreatmentLow Port -0.09542 0.17069 29.10724 -0.559 0.580
## TreatmentHigh Port -0.25601 0.18547 29.76387 -1.380 0.178
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.618
## TrtmntHghRf -0.696 0.553
## TrtmntLwPrt -0.720 0.529 0.673
## TrtmntHghPr -0.679 0.467 0.596 0.704 plot(fit3)
step(fit3)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 -15.400 44.801
## (1 | Well) 0 6 -17.855 47.710 4.9093 1 0.02671 *
## ---
## 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 1 0.22792 0.056979 4 15.788 1.0968 0.3921
##
## Model found:
## log(nmol.min) ~ (1 | Well)## model 4: LMER for exposure larva
fit4<-lmer(nmol.larva.min ~Treatment + (1|Replicate.vial) + (1|Plate..) + (1|Well),
data=data[data$Timepoint=="Recovery",])
fit4<-lmer(log(nmol.larva.min) ~Treatment + (1|Plate..),
data=data[data$Timepoint=="Recovery",])
anova(fit4) ranova(fit4) summary(fit4)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(nmol.larva.min) ~ Treatment + (1 | Plate..)
## Data: data[data$Timepoint == "Recovery", ]
##
## REML criterion at convergence: 50.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.73254 -0.50753 0.09466 0.55750 2.25525
##
## Random effects:
## Groups Name Variance Std.Dev.
## Plate.. (Intercept) 0.05322 0.2307
## Residual 0.25172 0.5017
## Number of obs: 32, groups: Plate.., 2
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -6.46320 0.25040 3.11305 -25.811 9.8e-05 ***
## TreatmentLow Reef 0.08070 0.27955 26.04884 0.289 0.775
## TreatmentHigh Reef 0.03975 0.26916 26.11666 0.148 0.884
## TreatmentLow Port -0.39489 0.27955 26.04884 -1.413 0.170
## TreatmentHigh Port -0.12991 0.27955 26.04884 -0.465 0.646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) TrtmLR TrtmHR TrtmLP
## TretmntLwRf -0.511
## TrtmntHghRf -0.537 0.473
## TrtmntLwPrt -0.511 0.463 0.473
## TrtmntHghPr -0.511 0.463 0.473 0.463 plot(fit4)
#step(fit4)
# Pairwise comparisons
# Day specific comparisons
# Sw.emmc<-emmeans(fit4, ~Treatment)
# Sw_groups<-cld(Sw.emmc)
# Sw_groupsPackages used
# Creates bibliography
#knitr::write_bib(c(.packages()), "packages.bib")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, and Martin Maechler. 2021. Matrix: Sparse and Dense Matrix Classes and Methods. http://Matrix.R-forge.R-project.org/.
Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2021. Lme4: Linear Mixed-Effects Models Using Eigen and S4. https://github.com/lme4/lme4/.
Fox, John. 2003. “Effect Displays in R for Generalised Linear Models.” Journal of Statistical Software 8 (15): 1–27. http://www.jstatsoft.org/v08/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. http://www.jstatsoft.org/v32/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. http://tinyurl.com/carbook.
Fox, John, Sanford Weisberg, and Brad Price. 2020. 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. 2019. 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. 2019. 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.
———. 2021. Multcomp: Simultaneous Inference in General Parametric Models. https://CRAN.R-project.org/package=multcomp.
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. 2019. 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. 2021. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Ripley, Brian. 2021. MASS: Support Functions and Datasets for Venables and Ripley’s Mass. http://www.stats.ox.ac.uk/pub/MASS4/.
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.
———. 2019. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.
———. 2021a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.
———. 2021b. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.
———. 2021c. 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. 2021. 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. 2021. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.
Wickham, Hadley, and Jim Hester. 2020. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.