diff --git a/LICENSE.md b/LICENSE.md
index 7e5a52a..bed44e0 100644
--- a/LICENSE.md
+++ b/LICENSE.md
@@ -1,6 +1,4 @@
-## creative commons
-
-# Attribution 4.0 International
+# Creative Commons Attribution 4.0 International
 
 Creative Commons Corporation (“Creative Commons”) is not a law firm and does not provide legal services or legal advice. Distribution of Creative Commons public licenses does not create a lawyer-client or other relationship. Creative Commons makes its licenses and related information available on an “as-is” basis. Creative Commons gives no warranties regarding its licenses, any material licensed under their terms and conditions, or any related information. Creative Commons disclaims all liability for damages resulting from their use to the fullest extent possible.
 
diff --git a/checklist.yml b/checklist.yml
index d0dd988..010abd9 100644
--- a/checklist.yml
+++ b/checklist.yml
@@ -18,3 +18,5 @@ spelling:
   other:
     en-GB:
       - source/r/geocomputations.R
+      - source/rmarkdown/analyses_diversity/test_breakaway_diversity.Rmd
+      - source/rmarkdown/compositional_analysis/test_rademu.Rmd
diff --git a/inst/en_gb.dic b/inst/en_gb.dic
new file mode 100644
index 0000000..a0e2923
--- /dev/null
+++ b/inst/en_gb.dic
@@ -0,0 +1,6 @@
+β
+Akker
+Mangiola
+modelmisspecification
+poisson
+sccomp
diff --git a/inst/nl_be.dic b/inst/nl_be.dic
index 6940d95..7e760cd 100644
--- a/inst/nl_be.dic
+++ b/inst/nl_be.dic
@@ -1,36 +1,92 @@
+'by
 +/-
 -SNE
 .Rdata
+.txt
+18S
+18s
+@example
+AICc
 Achaeta
+Aitchinson
 Annelida
+Beta
 Cmon
+Dicranocentrus
+Entomobryidae
+Folsomia
+GRTS
+HPC
 Henlea
+Isotomidae
+Isotomurus
 Lumbricus
 MBAG
 Marionina
+Megalothorax
+Names
 OTU
 OTUs
 Ordinaties
+Protaphorura
 Rarefaction
+Scientific
+WestEurope
+archaea
+arthropoda
+ascomycota
+basidiomycota
+be
 beta-binomiaal
 bimodaliteit
 bodembiodiversiteit
 bodemkoolstofmeetnet
+bodemstaal
+bodemstaalnamelocaties
+can
 centered-log-ratio
+chytridiomycota
+collembola
+collineaire
 counts
 covariaat
 covariaten
 credibiliteitsinterval
 csv
+design'
 design-variabelen
+diversity
 eDNA
+from
+fysico-chemische
+fysicochemische
 gepaardheid
 intercept
 macro-invertebraten
+meso-
 metabarcoding
+metazoa
+mollusca
 multivariaat
+names
 natuurgrasland
+natuurgraslanden
+nematoda
+nematodenstalen
+nolint
+nulwaarnemingen
+ordinatieruimte
+ordinaties
 phyloseq
 phylum
+platyhelminthes
+presence
+randomisatieprocedure
 rarefaction
+rarefy
+read
 reads
+result
+sub-sampled
+tardigrada
+which
diff --git a/source/r/check_presence.R b/source/r/check_presence.R
index 5fa58d5..5cef3c6 100644
--- a/source/r/check_presence.R
+++ b/source/r/check_presence.R
@@ -1,5 +1,5 @@
 #' Deze functie gaat na of een soort volgens GBIF voorkomt in een gekozen set
-#' van landen
+#' van landen # nolint start
 #'
 #' Default zijn West-Europese landen ingesteld: België, Frankrijk, Duitsland,
 #' Luxemburg, Nederland, Zwitserland, Oostenrijk
@@ -10,7 +10,7 @@
 #'
 #' @example
 #' result_WestEurope <- check_presence_from_file()
-#' print(result_WestEurope)
+#' print(result_WestEurope) # nolint end
 check_presence <- function(
     input = NULL,
     countries = c("BE", "FR", "DE", "LU", "NL", "CH", "AT")) {
@@ -35,20 +35,23 @@ check_presence <- function(
   for (country in countries) {
     # Function to check occurrence data and handle errors
     check_occurrence <- function(name) {
-      data <- tryCatch({
-        rgbif::occ_data(
-          scientificName = name,
-          country = country,  # Specify one country at a time
-          limit = 1
-        )
-      }, error = function(e) {
-        return(list(error = TRUE))
-      })
+      data <- tryCatch(
+        {
+          rgbif::occ_data(
+            scientificName = name,
+            country = country, # Specify one country at a time
+            limit = 1
+          )
+        },
+        error = function(e) {
+          return(list(error = TRUE))
+        }
+      )
 
       if (inherits(data, "error") || is.null(data$data)) {
-        return(FALSE)  # Occurrence data is not present or error occurred
+        return(FALSE) # Occurrence data is not present or error occurred
       } else {
-        return(TRUE)   # Occurrence data is present
+        return(TRUE) # Occurrence data is present
       }
     }
 
diff --git a/source/r/geocomputations.R b/source/r/geocomputations.R
index 157a889..03df253 100644
--- a/source/r/geocomputations.R
+++ b/source/r/geocomputations.R
@@ -10,11 +10,13 @@
 #'
 #' @examples
 point_to_gridcell <- function(
-  xy,
-  cell_width_m = 500,
-  point_position = c("center", "lowerleft", "upperleft", "lowerright",
-                     "upperright"),
-  crs = 31370) {
+    xy,
+    cell_width_m = 500,
+    point_position = c(
+      "center", "lowerleft", "upperleft", "lowerright",
+      "upperright"
+    ),
+    crs = 31370) {
   point_position <- match.arg(point_position)
 
   if (point_position != "center") stop(point_position, " not yet implemented")
@@ -24,13 +26,15 @@ point_to_gridcell <- function(
   xy <- sf::st_geometry(xy)
 
   # buffer with 1 point per quandrant
-  xy_buffer <- sf::st_buffer(x = xy,
-                             dist = cell_width_m / 2,
-                             nQuadSegs = 1)
+  xy_buffer <- sf::st_buffer(
+    x = xy,
+    dist = cell_width_m / 2,
+    nQuadSegs = 1
+  )
 
   # rotate 45 degrees around centroid
   rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
-  pl <- (xy_buffer - xy) * rot(pi/4) + xy
+  pl <- (xy_buffer - xy) * rot(pi / 4) + xy
   pl <- sf::st_sf(data.frame(xy_df, pl), crs = crs)
   return(pl)
 }
@@ -51,23 +55,24 @@ point_to_gridcell <- function(
 #'
 #' @examples
 landusemetrics_grid_cell <- function(
-  grid_cell,
-  layer,
-  grid_group_by_col = "POINT_ID",
-  layer_group_by_col = "",
-  progress = FALSE
-) {
+    grid_cell,
+    layer,
+    grid_group_by_col = "POINT_ID",
+    layer_group_by_col = "",
+    progress = FALSE) {
   require(duckdb)
-  if (inherits(layer, "SpatRaster") | inherits(layer, "RasterLayer")) {
+  require(dplyr)
+  if (inherits(layer, "SpatRaster") || inherits(layer, "RasterLayer")) {
     crs_grid <- gsub("^((.*?),\\n\\s*?){2}", "", sf::st_crs(grid_cell)$wkt)
     crs_layer <- gsub("^((.*?),\\n\\s*?){2}", "", terra::crs(layer))
     assertthat::assert_that(crs_grid == crs_layer)
 
     landcoverfraction <- function(df) {
       df %>%
-        mutate(frac_total = coverage_fraction / sum(coverage_fraction)) %>%
-        group_by(!!!syms(grid_group_by_col), value) %>%
-        summarize(freq = sum(frac_total), .groups = "drop_last")
+        mutate(frac_total = .data$coverage_fraction /
+          sum(.data$coverage_fraction)) %>%
+        group_by(!!!syms(grid_group_by_col), .data$value) %>%
+        summarize(freq = sum(.data$frac_total), .groups = "drop_last")
     }
 
     res <- exactextractr::exact_extract(
@@ -76,20 +81,20 @@ landusemetrics_grid_cell <- function(
       fun = landcoverfraction,
       summarize_df = TRUE,
       include_cols = grid_group_by_col,
-      progress = progress)
+      progress = progress
+    )
 
     return(res)
-
   }
 
   if (inherits(layer, "sf")) {
     assertthat::assert_that(sf::st_crs(grid_cell)$wkt == sf::st_crs(layer)$wkt)
 
-    int <- st_intersection(layer, grid_cell)
+    int <- sf::st_intersection(layer, grid_cell)
 
     cell_areas <- grid_cell %>%
       select(!!!syms(grid_group_by_col)) %>%
-      mutate(cell_area = sf::st_area(geometry)) %>%
+      mutate(cell_area = sf::st_area(.data$geometry)) %>%
       sf::st_drop_geometry()
 
     temparrow <- tempfile(fileext = ".parquet")
@@ -102,14 +107,15 @@ landusemetrics_grid_cell <- function(
 
     int <- arrow::open_dataset(temparrow) %>%
       arrow::to_duckdb() %>%
-      group_by(!!!syms(grid_group_by_col),
-               !!!syms(layer_group_by_col),
-               cell_area) %>%
-      summarise(area_m2 = sum(area)) %>%
-      mutate(area_prop = area_m2 / cell_area) %>%
+      group_by(
+        !!!syms(grid_group_by_col),
+        !!!syms(layer_group_by_col),
+        .data$cell_area
+      ) %>%
+      summarise(area_m2 = sum(.data$area)) %>%
+      mutate(area_prop = .data$area_m2 / .data$cell_area) %>%
       collect()
 
     return(int)
   }
 }
-
diff --git a/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd b/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd
index ba18c1c..8fb1526 100644
--- a/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd
+++ b/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd
@@ -16,11 +16,11 @@ data_subset <- combined %>%
 form_richness <- formula(
   observed ~
     log(total_count + 1)
-  + landgebruik
-  + diepte
-  + landgebruik:diepte
-  + (1 | cmon_plot_id)
-  )
+    + landgebruik
+    + diepte
+    + landgebruik:diepte
+    + (1 | cmon_plot_id)
+)
 
 ziform <- formula(~ 1 + (1 | cmon_plot_id))
 
@@ -28,35 +28,41 @@ cat("Fitting Poisson model\n")
 m_pois <- glmmTMB(
   formula = form_richness,
   data = data_subset,
-  family = poisson())
+  family = poisson()
+)
 cat("Fitting hurdle Poisson model\n")
 m_hpois <- glmmTMB(
   formula = form_richness,
   ziformula = ziform,
   data = data_subset,
-  family = truncated_poisson())
+  family = truncated_poisson()
+)
 cat("Fitting negative binomial model\n")
 m_nbinom <- glmmTMB(
   formula = form_richness,
   data = data_subset,
-  family = nbinom2())
+  family = nbinom2()
+)
 cat("Fitting hurdle negative binomial model\n")
 m_hnbinom <- glmmTMB(
   formula = form_richness,
   ziformula = ziform,
   data = data_subset,
-  family = truncated_nbinom2())
+  family = truncated_nbinom2()
+)
 cat("Fitting generalized Poisson model\n")
 m_genpois <- glmmTMB(
   formula = form_richness,
   data = data_subset,
-  family = genpois())
+  family = genpois()
+)
 cat("Fitting hurdle generalized Poisson model\n")
 m_hgenpois <- glmmTMB(
   formula = form_richness,
   ziformula = ziform,
   data = data_subset,
-  family = truncated_genpois())
+  family = truncated_genpois()
+)
 
 mlist <- list(
   pois = m_pois,
@@ -68,8 +74,8 @@ mlist <- list(
 )
 ```
 
-Het model met de laagste AIC waarde (hoogste AIC gewicht: AIC_wt) heeft de beste fit volgens AIC criterium.
-Dit betekent niet noodzakelijk dat dit model geen enkel probleem meer heeft met betrekking tot zero-inflation (deflation) of over- (onder-) dispersie.
+Het model met de laagste AIC waarde (hoogste AIC gewicht: `AIC_wt`) heeft de beste fit volgens AIC criterium.
+Dit betekent niet noodzakelijk dat dit model geen enkel probleem meer heeft met betrekking tot `zero-inflation` (`deflation`) of over- (onder-) dispersie.
 
 Modellen waarvan de AIC niet berekend kan worden, worden verwijderd.
 Deze modellen zijn wellicht niet geconvergeerd.
@@ -120,7 +126,8 @@ marginaleffects::plot_predictions(
   condition = c("total_count"),
   re.form = NA,
   vcov = TRUE,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Voorspeld aantal taxa")
 ```
 
@@ -130,25 +137,26 @@ marginaleffects::plot_predictions(
   condition = c("landgebruik", "diepte"),
   re.form = NA,
   vcov = TRUE,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Voorspeld aantal taxa")
 ```
 
 
-#### Met physicochemische data
+#### Met fysicochemische data
 
-Stepwise model selection:
+Stapsgewijze modelselectie:
 
 ```{r}
 variables_no_collinearity <- c(
-  "total_count", 
-  "landgebruik", 
-  "diepte", 
-  "cmon_plot_id", 
-  "c_density", 
-  "cn_stockbased", 
-  "swc_vol", 
-  "ph_kcl", 
+  "total_count",
+  "landgebruik",
+  "diepte",
+  "cmon_plot_id",
+  "c_density",
+  "cn_stockbased",
+  "swc_vol",
+  "ph_kcl",
   "bd"
 )
 
@@ -157,11 +165,12 @@ cat(
   paste(
     "De originele dataset heeft",
     nr1,
-  "rijen.\n"
+    "rijen.\n"
   )
 )
 data_subset <- data_subset[
-  complete.cases(data_subset[ , variables_no_collinearity]), ]
+  complete.cases(data_subset[, variables_no_collinearity]),
+]
 nr2 <- nrow(data_subset)
 
 cat(
@@ -169,7 +178,7 @@ cat(
     "Na verwijderen van rijen waarvan er ontbrekende gegevens zijn voor
     één van de covariaten, heeft de dataset",
     nr2,
-  "rijen.\n"
+    "rijen.\n"
   )
 )
 ```
@@ -182,13 +191,13 @@ Dit wordt het `full model`.
 ```{r extend-model-additional-predictors-{{group}}-{{primerset}}-{{unit}}}
 form_richness_full <- update(
   form_richness,
-  . ~ . 
-  + swc_vol
-  + c_density
-  + cn_stockbased
-  + swc_vol
-  + poly(ph_kcl, 2)
-  + bd
+  . ~ .
+    + swc_vol
+    + c_density
+    + cn_stockbased
+    + swc_vol
+    + poly(ph_kcl, 2)
+    + bd
 )
 
 m_richness_full <- update(
@@ -209,7 +218,6 @@ coefs_df <- rbind(
   data.frame(Model = "Base", base_model_coefs),
   data.frame(Model = "Full", full_model_coefs)
 )
-
 ```
 
 ```{r plot-coefficient-estimates-{{group}}-{{primerset}}-{{unit}}}
@@ -221,14 +229,18 @@ ggplot(
   geom_errorbar(
     aes(
       ymin = estimate - std.error,
-      ymax = estimate + std.error),
-    width = 0.2, position = position_dodge(0.9)) +
+      ymax = estimate + std.error
+    ),
+    width = 0.2, position = position_dodge(0.9)
+  ) +
   coord_flip() +
   theme_minimal() +
-  labs(title = "Coefficient Estimates and Standard Errors",
-       x = "Predictor",
-       y = "Estimate",
-       fill = "Model")
+  labs(
+    title = "Coefficient Estimates and Standard Errors",
+    x = "Predictor",
+    y = "Estimate",
+    fill = "Model"
+  )
 ```
 
 
@@ -241,9 +253,10 @@ model_set <- dredge(
   m.lim = c(NA, 7),
   fixed = ~
     cond(log(total_count + 1))
-  + cond(landgebruik)
-  + cond(diepte)
-  + (1 | cmon_plot_id))
+    + cond(landgebruik)
+      + cond(diepte)
+      + (1 | cmon_plot_id)
+)
 ```
 
 Per variabele berekenen we de som van modelgewichten over alle modellen en tellen we in hoeveel modellen de variabele voorkwam:
@@ -258,7 +271,7 @@ Visualisatie van modellen die minder dan 4 AICc verschillen van het model met la
 ```{r}
 dd <- subset(model_set, delta < 4)
 
-par(mar = c(3,5,6,4))
+par(mar = c(3, 5, 6, 4))
 plot(dd, labAsExpr = TRUE)
 ```
 
@@ -293,13 +306,17 @@ interactingvars <- stringr::str_split(interactingvars, ":") %>%
   unique()
 efvars <- efvars[!efvars %in% interactingvars]
 efvars <- stringr::str_replace_all(
-  efvars, "^(\\w+\\()?(.+?)(\\))?$", "\\2")
+  efvars, "^(\\w+\\()?(.+?)(\\))?$", "\\2"
+)
 efvars <- ifelse(
-  efvars == "diepte:landgebruik", "landgebruik:diepte", efvars)
+  efvars == "diepte:landgebruik", "landgebruik:diepte", efvars
+)
 efvars <- ifelse(
-  efvars == "total_count + 1", "total_count", efvars)
+  efvars == "total_count + 1", "total_count", efvars
+)
 efvars <- ifelse(
-  efvars == "ph_kcl, 2", "ph_kcl", efvars)
+  efvars == "ph_kcl, 2", "ph_kcl", efvars
+)
 
 p <- vector(mode = "list", length = length(efvars))
 p <- set_names(p, efvars)
@@ -309,7 +326,8 @@ for (i in efvars) {
     nd <- insight::get_datagrid(
       lowest_aic,
       by = conds,
-      length = 50)
+      length = 50
+    )
     pl <- marginaleffects::plot_predictions(
       model = lowest_aic,
       by = conds,
@@ -319,11 +337,12 @@ for (i in efvars) {
       points = 0.1,
       re.form = NA
     )
-    } else {
+  } else {
     nd <- insight::get_datagrid(
       lowest_aic,
       by = i,
-      length = 50)
+      length = 50
+    )
     pl <- marginaleffects::plot_predictions(
       lowest_aic,
       by = i,
@@ -333,7 +352,7 @@ for (i in efvars) {
       points = 0.1,
       re.form = NA
     )
-    }
+  }
   p[[i]] <- pl
 }
 
diff --git a/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd b/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd
index 25384f7..cf3f790 100644
--- a/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd
+++ b/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd
@@ -14,11 +14,11 @@ data_subset <- combined %>%
 form_shannon <- formula(
   shannon ~
     log(total_count)
-  + landgebruik
-  + diepte
-  + landgebruik:diepte
-  + (1 | cmon_plot_id)
-  )
+    + landgebruik
+    + diepte
+    + landgebruik:diepte
+    + (1 | cmon_plot_id)
+)
 
 zeroes <- min(data_subset$shannon, na.rm = TRUE) == 0
 
@@ -28,15 +28,20 @@ if (!zeroes) {
   m_shannon <- glmmTMB(
     formula = form_shannon,
     data = data_subset,
-    family = Gamma(link = "log"))
+    family = Gamma(link = "log")
+  )
 } else {
-  cat("Zeroes detected (only one taxum observed): Fitting zero-inflated Gamma model")
+  cat(
+    "Zeroes detected (only one taxum observed):
+    Fitting zero-inflated Gamma model"
+  )
 
   m_shannon <- glmmTMB(
     formula = form_shannon,
-    ziformula = ~ 1,
+    ziformula = ~1,
     data = data_subset,
-    family = ziGamma(link = "log"))
+    family = ziGamma(link = "log")
+  )
 }
 ```
 
@@ -60,7 +65,8 @@ marginaleffects::plot_predictions(
   condition = c("total_count"),
   re.form = NA,
   vcov = TRUE,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Shannon index") +
   scale_x_log10()
 ```
@@ -71,7 +77,8 @@ marginaleffects::plot_predictions(
   condition = c("landgebruik", "diepte"),
   re.form = NA,
   vcov = TRUE,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Shannon index")
 ```
 
diff --git a/source/rmarkdown/analyses_diversity/_child_model_simpson.Rmd b/source/rmarkdown/analyses_diversity/_child_model_simpson.Rmd
index 29c57ef..0457750 100644
--- a/source/rmarkdown/analyses_diversity/_child_model_simpson.Rmd
+++ b/source/rmarkdown/analyses_diversity/_child_model_simpson.Rmd
@@ -14,18 +14,19 @@ data_subset <- combined %>%
 form_simpson <- formula(
   simpson ~
     log(total_count)
-  + landgebruik
-  + diepte
-  + landgebruik:diepte
-  + (1 | cmon_plot_id)
-  )
+    + landgebruik
+    + diepte
+    + landgebruik:diepte
+    + (1 | cmon_plot_id)
+)
 
 cat("Fitting ordinal beta model")
 
 m_simpson <- glmmTMB(
   formula = form_simpson,
   data = data_subset,
-  family = ordbeta())
+  family = ordbeta()
+)
 ```
 
 
@@ -48,7 +49,8 @@ marginaleffects::plot_predictions(
   condition = c("total_count"),
   re.form = NA,
   vcov = TRUE,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Simpson index") +
   scale_x_log10()
 ```
@@ -59,7 +61,8 @@ marginaleffects::plot_predictions(
   condition = c("landgebruik", "diepte"),
   re.form = NA,
   vcov = TRUE,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Simpson index")
 ```
 
diff --git a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
index b06a783..49e47de 100644
--- a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
+++ b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
@@ -34,14 +34,10 @@ mbag_bodem_folder <- "G:/Gedeelde drives/PRJ_MBAG/4c_bodembiodiversiteit" # noli
 
 ```{r helper-zi-disp}
 inflation_dispersion <- function(model) {
-  test2 <- check_zeroinflation(model, tolerance = 0.1)
-  zeroes_within_tolerance <- between(
-    test2$ratio,
-    1 - test2$tolerance,
-    1 + test2$tolerance)
+  test2 <- performance::check_zeroinflation(model, tolerance = 0.1)
   zero_inflation <- test2$ratio < 1 - test2$tolerance
   zero_deflation <- test2$ratio > 1 + test2$tolerance
-  test <- check_overdispersion(model)
+  test <- performance::check_overdispersion(model)
   is_overdispersed <- test$p_value < 0.05 & test$dispersion_ratio > 1
   is_underdispersed <- test$p_value < 0.05 & test$dispersion_ratio < 1
   out <- list(
@@ -62,15 +58,17 @@ diversiteit <- readr::read_csv(
   file.path(
     mbag_bodem_folder,
     "data", "statistiek", "dataframe_overkoepelend",
-    "mbag_combined_dataframe_v2.csv")
+    "mbag_combined_dataframe_v2.csv"
+  )
 )
 
 metadata <- readr::read_csv(
   file.path(
     mbag_bodem_folder,
-    "data", 
+    "data",
     "Stratificatie_MBAG_plots",
-    "MBAG_stratfile_v2_cleaned_13.csv")
+    "MBAG_stratfile_v2_cleaned_13.csv"
+  )
 ) %>%
   janitor::clean_names() %>%
   rename(
@@ -86,11 +84,11 @@ metadata <- readr::read_csv(
       landgebruik_mbag,
       levels = c(
         "Akker", "Tijdelijk grasland", "Blijvend grasland",
-        "Residentieel grasland", "Natuurgrasland", "Heide", "Moeras")
+        "Residentieel grasland", "Natuurgrasland", "Heide", "Moeras"
+      )
     ),
     diepte = gsub("_|/", "-", diepte) |> factor()
   )
-
 ```
 
 In `mbag_combined_dataframe_v2.csv` zitten de samples zonder taxa niet in (na te vragen)!
@@ -99,7 +97,7 @@ In `mbag_combined_dataframe_v2.csv` zitten de samples zonder taxa niet in (na te
 sum(diversiteit$observed == 0)
 ```
 
-Deze nulwaarnemingen moeten we terug toevoegen (observed wordt dan 0, maar Shannon en Simpson zijn dan niet gedefinieerd).
+Deze nulwaarnemingen moeten we terug toevoegen (`observed` wordt dan 0, maar Shannon en Simpson zijn dan niet gedefinieerd).
 
 ```{r}
 groups <- diversiteit %>%
@@ -108,14 +106,15 @@ groups <- diversiteit %>%
 samples <- metadata %>%
   distinct(sample)
 
-#nrow(groups) * nrow(samples)
 samples_groups <- expand_grid(
-  samples, groups)
+  samples, groups
+)
 
 combined <- samples_groups %>%
   full_join(
     diversiteit,
-    by = c("sample", "primerset", "group", "unit")) %>%
+    by = c("sample", "primerset", "group", "unit")
+  ) %>%
   mutate(
     observed = ifelse(is.na(observed), 0, observed),
     total_count = ifelse(is.na(total_count), 0, total_count)
@@ -207,7 +206,8 @@ metadata %>%
     swc_vol, c_density, cn_stockbased, swc_vol, ph_kcl, bd
   ) %>%
   pivot_longer(
-    cols = c(swc_vol, c_density, cn_stockbased, swc_vol, ph_kcl, bd)) %>%
+    cols = c(swc_vol, c_density, cn_stockbased, swc_vol, ph_kcl, bd)
+  ) %>%
   ggplot(aes(x = landgebruik, y = value)) +
   geom_violin() +
   geom_sina(alpha = 0.3) +
@@ -234,12 +234,14 @@ inputs <- combined %>%
 pmap(
   inputs,
   function(
-    group = group,
-    primerset = primerset,
-    unit = unit) {
+      group = group,
+      primerset = primerset,
+      unit = unit) {
     knit_expand(
-      here("source", "rmarkdown", "analyses_diversity",
-           "_child_model_observed_richness.Rmd"),
+      here(
+        "source", "rmarkdown", "analyses_diversity",
+        "_child_model_observed_richness.Rmd"
+      ),
       group = group,
       primerset = primerset,
       unit = unit
@@ -248,7 +250,7 @@ pmap(
 ) %>%
   paste(collapse = "\n") -> rmd
 
-#clipr::write_clip(rmd)
+# clipr::write_clip(rmd) # nolint
 
 execute_code <- function(rmd) {
   # Extract R code from the rmd content
@@ -276,12 +278,14 @@ inputs <- combined %>%
 pmap(
   inputs,
   function(
-    group = group,
-    primerset = primerset,
-    unit = unit) {
+      group = group,
+      primerset = primerset,
+      unit = unit) {
     knit_expand(
-      here("source", "rmarkdown", "analyses_diversity",
-           "_child_model_shannon.Rmd"),
+      here(
+        "source", "rmarkdown", "analyses_diversity",
+        "_child_model_shannon.Rmd"
+      ),
       group = group,
       primerset = primerset,
       unit = unit
@@ -290,7 +294,7 @@ pmap(
 ) %>%
   paste(collapse = "\n") -> rmd
 
-#clipr::write_clip(rmd)
+# clipr::write_clip(rmd) # nolint
 
 knit_child(text = rmd, quiet = TRUE) %>%
   cat()
@@ -308,12 +312,14 @@ inputs <- combined %>%
 pmap(
   inputs,
   function(
-    group = group,
-    primerset = primerset,
-    unit = unit) {
+      group = group,
+      primerset = primerset,
+      unit = unit) {
     knit_expand(
-      here("source", "rmarkdown", "analyses_diversity",
-           "_child_model_simpson.Rmd"),
+      here(
+        "source", "rmarkdown", "analyses_diversity",
+        "_child_model_simpson.Rmd"
+      ),
       group = group,
       primerset = primerset,
       unit = unit
@@ -322,7 +328,7 @@ pmap(
 ) %>%
   paste(collapse = "\n") -> rmd
 
-#clipr::write_clip(rmd)
+# clipr::write_clip(rmd) # nolint
 
 knit_child(text = rmd, quiet = TRUE) %>%
   cat()
diff --git a/source/rmarkdown/analyses_diversity/test_breakaway_diversity.Rmd b/source/rmarkdown/analyses_diversity/test_breakaway_diversity.Rmd
index 4c00723..e3e5be5 100644
--- a/source/rmarkdown/analyses_diversity/test_breakaway_diversity.Rmd
+++ b/source/rmarkdown/analyses_diversity/test_breakaway_diversity.Rmd
@@ -29,15 +29,16 @@ mbag_bodem_folder <- "G:/Gedeelde drives/PRJ_MBAG/4c_bodembiodiversiteit" # noli
 library(breakaway)
 ```
 
-# Inlezen data
+# Read data
 
 ```{r inlezen}
 metadata <- readr::read_csv(
   file.path(
     mbag_bodem_folder,
-    "data", 
+    "data",
     "Stratificatie_MBAG_plots",
-    "MBAG_stratfile_v2_cleaned_13.csv")
+    "MBAG_stratfile_v2_cleaned_13.csv"
+  )
 ) %>%
   janitor::clean_names() %>%
   rename(
@@ -53,7 +54,8 @@ metadata <- readr::read_csv(
       landgebruik_mbag,
       levels = c(
         "Akker", "Tijdelijk grasland", "Blijvend grasland",
-        "Residentieel grasland", "Natuurgrasland", "Heide", "Moeras")
+        "Residentieel grasland", "Natuurgrasland", "Heide", "Moeras"
+      )
     ),
     diepte = gsub("_|/", "-", diepte) |> factor()
   )
@@ -63,19 +65,18 @@ load(
     mbag_bodem_folder,
     "data", "statistiek", "Annelida", "phyloseq",
     "physeq_Olig01_Annelida_species.Rdata"
-    )
   )
+)
 
-physeq_Olig01_Annelida_species <- physeq_Olig01_Annelida_species |> 
+physeq_olig01_annelida_species <- physeq_Olig01_Annelida_species |> # nolint
   phyloseq::subset_samples(
     !Landgebruik_MBAG %in% c("Moeras", "Heide")
   )
-
 ```
 
 
 ```{r}
-br <- breakaway(physeq_Olig01_Annelida_species)
+br <- breakaway(physeq_olig01_annelida_species)
 br_tbl <- summary(br)
 ```
 
@@ -85,24 +86,26 @@ br_tbl
 
 
 ```{r}
-combined <- br_tbl |> 
-  inner_join(metadata, by = join_by(sample_names == sample)) |> 
+combined <- br_tbl |>
+  inner_join(metadata, by = join_by(sample_names == sample)) |>
   filter(
     !landgebruik_mbag %in% c("Moeras", "Heide")
-  ) |> 
+  ) |>
   mutate(
     landgebruik_mbag = factor(
       landgebruik_mbag,
-      levels = c("Akker", "Tijdelijk grasland", "Blijvend grasland",
-                 "Residentieel grasland", "Natuurgrasland")
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
     )
   )
 ```
 
 
 ```{r}
-combined |> 
-  filter(error < 1e-7 ) %>%
+combined |>
+  filter(error < 1e-7) %>%
   count(landgebruik_mbag)
 ```
 
@@ -118,12 +121,12 @@ ma <- betta(
   formula = estimate ~ landgebruik_mbag,
   ses = error,
   data = combined
-  )
+)
 ma$table
 ```
 
 
-Instead, using meta-analytic capabilities of brms package:
+Instead, using meta-analytic capabilities of `brms` package:
 
 ```{r}
 library(brms)
@@ -135,12 +138,7 @@ ma_brms <- brm(
   family = "skew_normal",
   backend = "cmdstanr",
   cores = 4
-#  ,
-#  algorithm = "pathfinder",
-#  history_size = 100,
-#  psis_resample = TRUE
 )
-
 ```
 
 ```{r}
@@ -151,7 +149,7 @@ summary(ma_brms)
 conditional_effects(
   ma_brms,
   "landgebruik_mbag:diepte"
-  )
+)
 ```
 
-The result is very much comparable to what we obtained via glmmTMB with an offset term - so probably not much added value.
+The result is very much comparable to what we obtained via `glmmTMB` with an offset term - so probably not much added value.
diff --git a/source/rmarkdown/annelida_compositional_differences.Rmd b/source/rmarkdown/annelida_compositional_differences.Rmd
index 02c4a46..c130c36 100644
--- a/source/rmarkdown/annelida_compositional_differences.Rmd
+++ b/source/rmarkdown/annelida_compositional_differences.Rmd
@@ -38,9 +38,10 @@ mbag_folder <- "G:/Gedeelde drives/PRJ_MBAG" # nolint
 
 
 if (!"GenomeInfoDbData" %in% rownames(installed.packages())) {
-install.packages(
-  "GenomeInfoDbData",
-  repos = "https://bioconductor.org/packages/3.17/data/annotation")
+  install.packages(
+    "GenomeInfoDbData",
+    repos = "https://bioconductor.org/packages/3.17/data/annotation"
+  )
 }
 if (!"sccomp" %in% rownames(installed.packages())) {
   remotes::install_github("stemangiola/sccomp")
@@ -60,14 +61,14 @@ Deze tabelletjes kunnen best gewoon csv bestanden zijn die we inlezen.
 ```
 ```{r load-rdata-file}
 path_naar_bestand <- file.path(
-    mbag_folder,
-    "4c_bodembiodiversiteit",
-    "data",
-    "statistiek",
-    "Annelida",
-    "phyloseq",
-    "physeq_Olig01_Annelida_species.Rdata"
-  )
+  mbag_folder,
+  "4c_bodembiodiversiteit",
+  "data",
+  "statistiek",
+  "Annelida",
+  "phyloseq",
+  "physeq_Olig01_Annelida_species.Rdata"
+)
 
 load(path_naar_bestand)
 ```
@@ -88,10 +89,12 @@ tidy_physeq <- tidytacos::from_phyloseq(physeq_Olig01_Annelida_species) %>%
   tidytacos::mutate_samples(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 ```
 
@@ -113,7 +116,7 @@ tidy_physeq %>%
       paste(collapse = "; ") %>%
       stringr::str_sub(start = 1, end = 13),
     period = interval(start = min(Datum_staalname), end = max(Datum_staalname)),
-    ) %>%
+  ) %>%
   kable()
 ```
 
@@ -164,7 +167,8 @@ tidy_physeq %>%
   geom_point(aes(x = Datum_staalname, y = total_count, colour = Staalnemer)) +
   geom_smooth(aes(
     x = Datum_staalname, y = total_count,
-    colour = Staalnemer, fill = Staalnemer)) +
+    colour = Staalnemer, fill = Staalnemer
+  )) +
   scale_y_log10(breaks = c(1e+04, 1e+05, 3e+05, 1e+06))
 ```
 
@@ -177,7 +181,7 @@ We definiëren een functie die volgende stappen uitvoert:
     -   nadeel: indien zeldzame taxa geclusterd zitten in een groep, missen we dit patroon als we ze verwijderen
 -   (optioneel) "rarefy" de samples tot een gelijk aantal reads (instelbaar)
     -   voordeel: eventueel vertekening door verschillen in totaal aantal reads tussen stalen wordt weggewerkt; kan zeker van belang zijn als er systematische verschillen in totaal aantal reads zitten tussen de groepen die we willen vergelijken (zie vorige sectie)
-    -   nadeel: dit is een randomisatieprocedure die in feite 10 - 100-tal keer zou herhaald moeten worden waarna de resultaten van elk van de random datasets gecombineerd worden (vergelijkbaar met "multiple imputation"); bij een ordinatie is er echter geen manier om tot een consensus ordinatie te komen
+    -   nadeel: dit is een randomisatieprocedure die in feite 10 - 100-tal keer zou herhaald moeten worden waarna de resultaten van elk van de random datasets gecombineerd worden (vergelijkbaar met `"multiple imputation"`); bij een ordinatie is er echter geen manier om tot een consensus ordinatie te komen
 -   We transformeren met de robuuste "centered-log-ratio" transformatie omdat we werken met compositionele data waar enkel relatieve verschillen in abundantie betekenisvol zijn
 -   Eerst voeren we een principale componenten analyse uit op deze data
     -   de afstanden tussen samples in de ordinatieruimte zijn dan Euclidische afstanden tussen "centered-log-ratio" getransformeerde variabelen - dit is in de literatuur ook bekend als de Aitchinson afstand.
@@ -191,8 +195,7 @@ calc_ordinations <- function(
     rarefy_even_reads = NULL,
     frequency_threshold = 0,
     cumulative_prop_explained = 0.9,
-    use_rank = NULL
-) {
+    use_rank = NULL) {
   set.seed(20315)
 
   data <- tacos_data %>%
@@ -204,9 +207,10 @@ calc_ordinations <- function(
       }
     } %>%
     tidytacos::filter_counts(
-      .by = "taxon_id", n() > frequency_threshold)
+      .by = "taxon_id", n() > frequency_threshold
+    )
+
 
-  
   if (!is.null(rarefy_even_reads)) {
     message(paste("samples rarefied to", rarefy_even_reads, "reads"))
     data <- data %>%
@@ -220,17 +224,21 @@ calc_ordinations <- function(
       id_cols = "sample_id",
       names_from = "taxon_id",
       values_from = "count",
-      values_fill = 0)
+      values_fill = 0
+    )
 
   pca1 <- prcomp(
     vegan::decostand(
       counts %>% select(-"sample_id"),
-      method = "rclr"))
+      method = "rclr"
+    )
+  )
   d <- pca1$sdev^2
   truncate_axis <- sum(cumsum(d) / sum(d) < cumulative_prop_explained)
   message(sprintf(
     "Calculating tsne based on first %s axes from PCA solution",
-    truncate_axis))
+    truncate_axis
+  ))
   tsne1 <- tsne::tsne(pca1$x[, 1:truncate_axis])
 
   plotdata <-
@@ -238,7 +246,8 @@ calc_ordinations <- function(
       as_tibble(vegan::scores(pca1, display = "sites", choices = 1:2)),
       as_tibble(tsne1),
       counts %>%
-        dplyr::select(.data$sample_id)) %>%
+        dplyr::select(.data$sample_id)
+    ) %>%
     left_join(data$samples, by = "sample_id")
 
   return(list(pca = pca1, tsne = tsne1, plotdata = plotdata))
@@ -248,30 +257,39 @@ plot_pca <- function(plotdata, colour = NULL) {
   ggpca <- plotdata %>%
     ggplot() +
     geom_point(
-      aes(x = .data$PC1, y = .data$PC2, colour = {{colour}})
+      aes(x = .data$PC1, y = .data$PC2, colour = {{ colour }})
     ) +
     geom_vline(xintercept = 0, alpha = 0.3) +
     geom_hline(yintercept = 0, alpha = 0.3) +
-    stat_ellipse(aes(
-      x = .data$PC1, y = .data$PC2, colour = {{colour}}),
-      level = 0.95) +
-    labs(x = "Eerste principale component",
-         y = "Tweede principale component")
-    return(ggpca)
+    stat_ellipse(
+      aes(
+        x = .data$PC1, y = .data$PC2, colour = {{ colour }}
+      ),
+      level = 0.95
+    ) +
+    labs(
+      x = "Eerste principale component",
+      y = "Tweede principale component"
+    )
+  return(ggpca)
 }
 
 plot_tsne <- function(plotdata, colour = NULL) {
   ggtsne <- plotdata %>%
     ggplot() +
     geom_point(
-      aes(x = .data$V1, y = .data$V2, colour = {{colour}})) +
+      aes(x = .data$V1, y = .data$V2, colour = {{ colour }})
+    ) +
     geom_vline(xintercept = 0, alpha = 0.3) +
     geom_hline(yintercept = 0, alpha = 0.3) +
     stat_ellipse(aes(
-      x = .data$V1, y = .data$V2, colour = {{colour}}), level = 0.95) +
-    labs(x = "Eerste t-SNE as",
-         y = "Tweede t-SNE as")
-    return(ggtsne)
+      x = .data$V1, y = .data$V2, colour = {{ colour }}
+    ), level = 0.95) +
+    labs(
+      x = "Eerste t-SNE as",
+      y = "Tweede t-SNE as"
+    )
+  return(ggtsne)
 }
 ```
 
@@ -281,7 +299,8 @@ plot_tsne <- function(plotdata, colour = NULL) {
 ordination <- calc_ordinations(
   tacos_data = tidy_physeq,
   frequency_threshold = 0,
-  rarefy_even_reads = 100000)
+  rarefy_even_reads = 100000
+)
 ```
 
 ```{r}
@@ -312,7 +331,8 @@ ordination <- calc_ordinations(
   tacos_data = tidy_physeq,
   use_rank = "species",
   frequency_threshold = 0,
-  rarefy_even_reads = 100000)
+  rarefy_even_reads = 100000
+)
 ```
 
 ```{r}
@@ -342,7 +362,8 @@ ordination <- calc_ordinations(
   tacos_data = tidy_physeq,
   use_rank = "genus",
   frequency_threshold = 0,
-  rarefy_even_reads = 30780)
+  rarefy_even_reads = 30780
+)
 ```
 
 ```{r}
@@ -384,8 +405,10 @@ sccomp_data_annelida <- tidy_physeq %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, genus, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, genus, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -401,9 +424,9 @@ if (!file.exists(here("source", "brms_models", "s_annelida.rds"))) {
     .data = sccomp_data_annelida,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -411,8 +434,10 @@ if (!file.exists(here("source", "brms_models", "s_annelida.rds"))) {
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_annelida,
-          here("source", "brms_models", "s_annelida.rds"))
+  saveRDS(
+    s_annelida,
+    here("source", "brms_models", "s_annelida.rds")
+  )
 }
 
 s_annelida <- readRDS(here("source", "brms_models", "s_annelida.rds"))
@@ -483,25 +508,30 @@ s_annelida_predictions <- sccomp::sccomp_predict(
   s_annelida,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_annelida %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_annelida %>%
-              distinct(taxon_id, genus),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_annelida %>%
+      distinct(taxon_id, genus),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_annelida_predictions <- s_annelida_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_annelida %>%
-              distinct(taxon_id, genus),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_annelida %>%
+      distinct(taxon_id, genus),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     genus = factor(genus, levels = sortering$genus)
@@ -521,11 +551,12 @@ s_annelida_predictions %>%
       y = proportion_mean,
       ymin = proportion_lower,
       ymax = proportion_upper,
-      colour = Diepte),
+      colour = Diepte
+    ),
     position = position_dodge(width = 0.4),
     alpha = 0.7
   ) +
-  facet_wrap(~ genus, scales = "free_y") +
+  facet_wrap(~genus, scales = "free_y") +
   scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -536,9 +567,12 @@ We kunnen dit (zonder credibiliteitsinterval) ook als volgt voorstellen:
 s_annelida_predictions %>%
   ggplot() +
   geom_col(
-    aes(x = Landgebruik_MBAG,
-        y = proportion_mean,
-        fill = genus)) +
+    aes(
+      x = Landgebruik_MBAG,
+      y = proportion_mean,
+      fill = genus
+    )
+  ) +
   facet_grid(~Diepte) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -562,8 +596,10 @@ sccomp_data_annelida <- tidy_physeq %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, species, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, species, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -572,7 +608,7 @@ sccomp_data_annelida <- tidy_physeq %>%
 Voor deze analyse, aggregeren we de taxa tot `r used_rank` niveau en beperken we de analyse tot de `r max_taxa` taxa met hoogste prevalentie.
 Dit is goed voor `r round(sum(sccomp_data_annelida$count) / sum(tidy_physeq$counts$count) * 100)` % van alle reads.
 
-Hier dus eens getest met de 293 meest prevalente OTUs
+Hier dus eens getest met de 293 meest voorkomende OTUs
 
 ```{r sccomp-annelida-all-species}
 fs::dir_create(here("source", "brms_models"))
@@ -581,9 +617,9 @@ if (!file.exists(here("source", "brms_models", "s_annelida_293_species.rds"))) {
     .data = sccomp_data_annelida,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -591,11 +627,15 @@ if (!file.exists(here("source", "brms_models", "s_annelida_293_species.rds"))) {
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_annelida_293_species,
-          here("source", "brms_models", "s_annelida_293_species.rds"))
+  saveRDS(
+    s_annelida_293_species,
+    here("source", "brms_models", "s_annelida_293_species.rds")
+  )
 }
 
-s_annelida_293_species <- readRDS(here("source", "brms_models", "s_annelida_293_species.rds"))
+s_annelida_293_species <- readRDS(
+  here("source", "brms_models", "s_annelida_293_species.rds")
+)
 ```
 
 ### Model evaluatie
@@ -622,8 +662,8 @@ Als de blauwe boxen goed overeenkomen met de geobserveerde boxen en data, dan is
 Deze boxplot figuren leren ons dat het model de data goed kan beschrijven, maar dat er mogelijkheden zijn om het model te verbeteren:
 
 -   er zijn mogelijk verschillen in variantie tussen de groepen: we kunnen dus eventueel de variantie ook modelleren in functie van deze factoren
--   Sprake van bimodaliteit? In 't292' ( _xx_ _xx_) akker? Maar wel verdwenen voor _xx_ ( _xx_) dus in natuurgraslanden?: de boxplot zit in de zone met proporties tussen 0 en 0.25, terwijl er een groep van punten is die clustert tussen 0.25 en 1, en vooral tusen 0.5 en 1. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen)?
--   Idem voor _xx_ _xx_ 't298' in tijdelijk grasland? de boxplot zit in de zone met proporties tussen x en x, terwijl er een groep van punten is die clustert tussen y en y
+-   Sprake van bimodaliteit? In `'t292'` (`_xx_ _xx_`) akker? Maar wel verdwenen voor `_xx_` (`_xx_`) dus in natuurgraslanden?: de boxplot zit in de zone met proporties tussen 0 en 0.25, terwijl er een groep van punten is die clustert tussen 0.25 en 1, en vooral tussen 0.5 en 1. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen)?
+-   Idem voor `_xx_ _xx_` `'t298'` in tijdelijk grasland? de boxplot zit in de zone met proporties tussen x en x, terwijl er een groep van punten is die clustert tussen y en y
 
 ```{r sccomp-boxplot-landgebruik-species}
 plots$boxplot[[1]]
@@ -664,25 +704,30 @@ s_annelida_predictions <- sccomp::sccomp_predict(
   s_annelida_293_species,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_annelida_293_species %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_annelida %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_annelida %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_annelida_predictions <- s_annelida_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_annelida %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_annelida %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     species = factor(species, levels = sortering$species)
@@ -690,8 +735,8 @@ s_annelida_predictions <- s_annelida_predictions %>%
 ```
 
 De volgende figuur toont de modelvoorspellingen.
-Duidelijk geen toename van _xx_ _xx_ met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. De voorheen geobserveerde toename van _xx_ met afnemende landgebruiksintensiteit moet dus aan andere _xx_ soorten te wijten zijn? Duidelijke afname van _xx_ _xx_ en _xx_ _xx_ (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit
-_xx_ _xx_ komt duidelijk meer voor in 10-30 cm dan 0-10 cm.
+Duidelijk geen toename van `_xx_ _xx_` met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. De voorheen geobserveerde toename van `_xx_` met afnemende landgebruiksintensiteit moet dus aan andere `_xx_` soorten te wijten zijn? Duidelijke afname van `_xx_ _xx_` en `_xx_ _xx_` (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit
+`_xx_ _xx_` komt duidelijk meer voor in 10-30 cm dan 0-10 cm.
 
 ```{r plot-sccomp-predictions-species}
 # Filter out species with "_otu" in their names
@@ -700,15 +745,15 @@ filtered_species_data <- s_annelida_predictions %>%
 
 # Find unique species from the filtered data and split them into chunks of 9
 unique_species <- unique(filtered_species_data$species)
-species_chunks <- split(unique_species, ceiling(seq_along(unique_species)/9))
+species_chunks <- split(unique_species, ceiling(seq_along(unique_species) / 9))
 
 # Loop through each chunk and create a plot
 plots <- list() # To store the plots
 for (i in seq_along(species_chunks)) {
   # Subset the data for the current chunk of species
-  subset_data <- filtered_species_data %>% 
+  subset_data <- filtered_species_data %>%
     filter(species %in% species_chunks[[i]])
-  
+
   # Create the plot for the current subset
   p <- subset_data %>%
     ggplot() +
@@ -718,22 +763,20 @@ for (i in seq_along(species_chunks)) {
         y = proportion_mean,
         ymin = proportion_lower,
         ymax = proportion_upper,
-        colour = Diepte),
+        colour = Diepte
+      ),
       position = position_dodge(width = 0.4),
       alpha = 0.7
     ) +
-    facet_wrap(~ species, scales = "free_y", ncol = 3) +
+    facet_wrap(~species, scales = "free_y", ncol = 3) +
     scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
     theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
-  
+
   # Save the plot in the list
   plots[[i]] <- p
 }
 
 # Now you can print or save each plot separately
-# For example, to view the first set of plots, you can do:
-# print(plots[[1]])
-# And so on for each plot in 'plots'
 ```
 
 ```{r}
@@ -778,9 +821,12 @@ We kunnen dit (zonder credibiliteitsinterval) ook als volgt voorstellen:
 s_annelida_predictions %>%
   ggplot() +
   geom_col(
-    aes(x = Landgebruik_MBAG,
-        y = proportion_mean,
-        fill = species)) +
+    aes(
+      x = Landgebruik_MBAG,
+      y = proportion_mean,
+      fill = species
+    )
+  ) +
   facet_grid(~Diepte) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
diff --git a/source/rmarkdown/annelida_data_analyse.Rmd b/source/rmarkdown/annelida_data_analyse.Rmd
index a32a996..ff8e410 100644
--- a/source/rmarkdown/annelida_data_analyse.Rmd
+++ b/source/rmarkdown/annelida_data_analyse.Rmd
@@ -65,14 +65,14 @@ Deze tabelletjes kunnen best gewoon csv bestanden zijn die we inlezen.
 
 ```{r load-rdata-file}
 path_naar_bestand <- file.path(
-    mbag_folder,
-    "4c_bodembiodiversiteit",
-    "data",
-    "statistiek",
-    "Annelida",
-    "phyloseq",
-    "physeq_Olig01_Annelida_species.Rdata"
-  )
+  mbag_folder,
+  "4c_bodembiodiversiteit",
+  "data",
+  "statistiek",
+  "Annelida",
+  "phyloseq",
+  "physeq_Olig01_Annelida_species.Rdata"
+)
 
 load(path_naar_bestand)
 ```
@@ -97,10 +97,10 @@ Filter taxongroep:
 ```{r filter-taxongroep}
 physeq <- subset_taxa(
   physeq,
-  phylum == "Annelida")
+  phylum == "Annelida"
+)
 
 physeq_rarefied <- rarefy_even_depth(physeq, sample.size = 30780)
-
 ```
 
 
@@ -120,7 +120,6 @@ psotu2veg <- function(physeq) {
   }
   return(as(otu, "matrix"))
 }
-
 ```
 
 
@@ -147,7 +146,6 @@ tidy_physeq_rarefied <- tidy_physeq_rarefied %>%
   ) %>%
   tidytacos::add_alpha() %>%
   tidytacos::add_total_count()
-
 ```
 
 
@@ -206,27 +204,24 @@ tidy_physeq$samples %>%
 ## OTU tabel
 
 ```{r verken-otu-data}
-
-
-glimpse(otu_table(physeq_rarefied) %>% as.data.frame %>% as_tibble())
+glimpse(otu_table(physeq_rarefied) %>% as.data.frame() %>% as_tibble())
 ```
 
 ## Taxonomie tabel
 
 
 ```{r verken-taxonomie-data}
-glimpse(tax_table(physeq_rarefied) %>% as.data.frame %>% as_tibble())
+glimpse(tax_table(physeq_rarefied) %>% as.data.frame() %>% as_tibble())
 ```
 
 ```{r GBIF-check-presence}
-
 species <- as.vector(tax_table(physeq)[, "species"])
 species <- species[!grepl("_otu", species)]
 species <- gsub("_", " ", species)
 source(here::here("source/r/check_presence.R"))
 gbif_check <- check_presence(species)
 gbif_check %>%
-filter(!present) %>%
+  filter(!present) %>%
   kable(caption = "Not present according to GBIF in Western-Europe")
 ```
 
@@ -251,7 +246,10 @@ tidy_physeq_rarefied %>%
 
 ```{r}
 # Create a new combined column
-combined_column <- paste(sample_data(physeq_rarefied)$Landgebruik_MBAG, sample_data(physeq_rarefied)$Diepte, sep = "_")
+combined_column <- paste(sample_data(physeq_rarefied)$Landgebruik_MBAG,
+  sample_data(physeq_rarefied)$Diepte,
+  sep = "_"
+)
 
 # Add the new combined column to the sample metadata with the desired name
 sample_data(physeq_rarefied)$Landgebruik_MBAG_diepte <- combined_column
@@ -260,12 +258,17 @@ sample_data(physeq_rarefied)$Landgebruik_MBAG_diepte <- combined_column
 physeq_rarefied <- merge_phyloseq(physeq_rarefied, sample_data(physeq_rarefied))
 
 
-if (system(command = "which ktImportText",  intern = FALSE, ignore.stderr = TRUE,  ignore.stdout = TRUE) != 1) {
-  psadd::plot_krona(physeq_rarefied, "MBAG_Olig01_Annelida_rar_species_alle_stalen_Landgebruik_MBAG_per_diepte", "Landgebruik_MBAG_diepte", trim = T)
+if (system(
+  command = "which ktImportText", intern = FALSE, ignore.stderr = TRUE,
+  ignore.stdout = TRUE
+) != 1) {
+  psadd::plot_krona(
+    physeq_rarefied,
+    "MBAG_Olig01_Annelida_rar_species_alle_stalen_Landgebruik_MBAG_per_diepte",
+    "Landgebruik_MBAG_diepte",
+    trim = TRUE
+  )
 }
-
-
-
 ```
 
 
@@ -290,7 +293,8 @@ Het geobserveerde aantal taxa neemt toe met het aantal reads in een staal:
 tidy_physeq$samples %>%
   filter(total_count > 0) %>%
   ggplot(
-    aes(x = total_count, y = observed)) +
+    aes(x = total_count, y = observed)
+  ) +
   geom_point() +
   scale_x_log10() +
   labs(y = "Observed number of taxa")
@@ -310,16 +314,17 @@ samples_totcount_10k <- tidy_physeq$samples %>%
 form_annelida_richness <- formula(
   observed ~
     log(total_count)
-  + Landgebruik_MBAG
-  + Diepte
-  + Landgebruik_MBAG:Diepte
-  + (1 | Cmon_PlotID)
-  )
+    + Landgebruik_MBAG
+    + Diepte
+    + Landgebruik_MBAG:Diepte
+    + (1 | Cmon_PlotID)
+)
 
 m_annelida_richness <- glmmTMB::glmmTMB(
   formula = form_annelida_richness,
   data = samples_totcount_10k,
-  family = glmmTMB::nbinom2())
+  family = glmmTMB::nbinom2()
+)
 ```
 
 Model validatie toont dat er problemen zijn met de residuele variabiliteit.
@@ -345,7 +350,8 @@ Zoals te verwachten, is het belangrijk om te corrigeren voor het totaal aantal r
 marginaleffects::plot_predictions(
   m_annelida_richness,
   condition = c("total_count"),
-  vcov = vcov(m_annelida_richness)) +
+  vcov = vcov(m_annelida_richness)
+) +
   labs(y = "Voorspeld aantal taxa") +
   scale_x_log10()
 ```
@@ -359,7 +365,8 @@ marginaleffects::plot_predictions(
   condition = c("Landgebruik_MBAG", "Diepte"),
   vcov = vcov(m_annelida_richness),
   re.form = NA,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Voorspelde aantal taxa")
 ```
 
@@ -377,16 +384,20 @@ estimated_species_richness <- vegan::rarefy(
   veganobject,
   30780,
   se = TRUE,
-  MARGIN = 1)
+  MARGIN = 1
+)
 
 
 true_rarefaction_results <- cbind(
   sam_new %>% select(Cmon_PlotID, Diepte, Landgebruik_MBAG),
   est = estimated_species_richness["S", ],
-  se = estimated_species_richness["se", ]) %>%
+  se = estimated_species_richness["se", ]
+) %>%
   as_tibble() %>%
-  mutate(est_lower = est - 2 * se,
-         est_upper = est + 2 * se)
+  mutate(
+    est_lower = est - 2 * se,
+    est_upper = est + 2 * se
+  )
 ```
 
 
@@ -402,17 +413,17 @@ library(brms)
 
 form_rarefaction <- bf(
   est | se(se, sigma = TRUE) ~
-  + Landgebruik_MBAG
-  + Diepte
-  + Landgebruik_MBAG:Diepte
-  + (1 | Cmon_PlotID)
+    +Landgebruik_MBAG
+    + Diepte
+    + Landgebruik_MBAG:Diepte
+    + (1 | Cmon_PlotID)
 )
 
 m_rarefaction <- brm(
   formula = form_rarefaction,
   data = true_rarefaction_results,
-  family = gaussian())
-
+  family = gaussian()
+)
 ```
 
 
@@ -426,8 +437,9 @@ Dit model geeft gelijkaardige resultaten als het model waarbij we rechtstreeks w
 ```{r m-rarefied-predictions}
 conditional_effects(
   m_rarefaction,
-    effects = c("Landgebruik_MBAG:Diepte"),
-  re.form = NA)
+  effects = c("Landgebruik_MBAG:Diepte"),
+  re.form = NA
+)
 ```
 
 ## OTU counts 
@@ -486,7 +498,6 @@ tidy_physeq %>%
 ## Ordinatie
 
 ```{r physec-rarefied}
-
 ```
 
 ```{r ordination-vegan}
@@ -505,8 +516,6 @@ samples <- sample_data(physeq_rarefied) %>%
 rda <- rda(ordmat)
 
 biplot(rda)
-
-
 ```
 
 
@@ -522,7 +531,9 @@ pcoa_result <- ordinate(physeq, method = "PCoA", distance = "bray")
 
 # Create the ordination plot with color based on land_use and shape
 # based on diepte
-plot_ordination(physeq, pcoa_result, type = "samples",
-                color = "Landgebruik_MBAG", shape = "Diepte") +
+plot_ordination(physeq, pcoa_result,
+  type = "samples",
+  color = "Landgebruik_MBAG", shape = "Diepte"
+) +
   theme_minimal()
 ```
diff --git a/source/rmarkdown/collembola_compositional_differences.Rmd b/source/rmarkdown/collembola_compositional_differences.Rmd
index 4eeb1fc..055d647 100644
--- a/source/rmarkdown/collembola_compositional_differences.Rmd
+++ b/source/rmarkdown/collembola_compositional_differences.Rmd
@@ -1,5 +1,5 @@
 ---
-title: "Collembola: verschillen in taxonomische samenstelling tussen groepen van stalen"
+title: "`Collembola`: verschillen in taxonomische samenstelling tussen groepen van stalen"
 author: "Sam Lambrechts, Io Deflem, Emma Cartuyvels, Hans Van Calster"
 date: "`r Sys.Date()`"
 output:
@@ -39,9 +39,10 @@ mbag_folder <- "G:/Gedeelde drives/PRJ_MBAG" # nolint
 
 
 if (!"GenomeInfoDbData" %in% rownames(installed.packages())) {
-install.packages(
-  "GenomeInfoDbData",
-  repos = "https://bioconductor.org/packages/3.17/data/annotation")
+  install.packages(
+    "GenomeInfoDbData",
+    repos = "https://bioconductor.org/packages/3.17/data/annotation"
+  )
 }
 if (!"sccomp" %in% rownames(installed.packages())) {
   remotes::install_github("stemangiola/sccomp")
@@ -61,15 +62,15 @@ Deze tabelletjes kunnen best gewoon csv bestanden zijn die we inlezen.
 ```
 ```{r load-rdata-file}
 path_naar_bestand <- file.path(
-    mbag_folder,
-    "4c_bodembiodiversiteit",
-    "data",
-    "statistiek",
-    "Coll01",
-    "met_otu_legetaxonomie",
-    "phyloseq",
-    "physeq_Coll01_Collembola_species.Rdata"
-  )
+  mbag_folder,
+  "4c_bodembiodiversiteit",
+  "data",
+  "statistiek",
+  "Coll01",
+  "met_otu_legetaxonomie",
+  "phyloseq",
+  "physeq_Coll01_Collembola_species.Rdata"
+)
 
 load(path_naar_bestand)
 ```
@@ -81,9 +82,13 @@ Het bestand dat wordt ingelezen is `r basename(path_naar_bestand)`.
 ```{r create-tidytacos-objects}
 physeq <- physeq_Coll01_Collembola.species
 
-physeq_known_species <- subset_taxa(physeq, !grepl("_otu", tax_table(physeq)[, "species"]))
+physeq_known_species <- subset_taxa(
+  physeq, !grepl("_otu", tax_table(physeq)[, "species"])
+)
 
-physeq_known_genera <- subset_taxa(physeq, !grepl("_otu", tax_table(physeq)[, "genus"]))
+physeq_known_genera <- subset_taxa(
+  physeq, !grepl("_otu", tax_table(physeq)[, "genus"])
+)
 
 
 
@@ -98,10 +103,12 @@ tidy_physeq <- tidytacos::from_phyloseq(physeq) %>%
   tidytacos::mutate_samples(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 
 ranks <- physeq_known_species %>%
@@ -115,10 +122,12 @@ tidy_physeq_known_species <- tidytacos::from_phyloseq(physeq_known_species) %>%
   tidytacos::mutate_samples(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 
 ranks <- physeq_known_genera %>%
@@ -132,10 +141,12 @@ tidy_physeq_known_genera <- tidytacos::from_phyloseq(physeq_known_genera) %>%
   tidytacos::mutate_samples(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 ```
 
@@ -157,7 +168,7 @@ tidy_physeq %>%
       paste(collapse = "; ") %>%
       stringr::str_sub(start = 1, end = 13),
     period = interval(start = min(Datum_staalname), end = max(Datum_staalname)),
-    ) %>%
+  ) %>%
   kable()
 ```
 
@@ -208,7 +219,8 @@ tidy_physeq %>%
   geom_point(aes(x = Datum_staalname, y = total_count, colour = Staalnemer)) +
   geom_smooth(aes(
     x = Datum_staalname, y = total_count,
-    colour = Staalnemer, fill = Staalnemer)) +
+    colour = Staalnemer, fill = Staalnemer
+  )) +
   scale_y_log10(breaks = c(1e+03, 1e+04, 1e+05, 3e+05, 1e+06))
 ```
 
@@ -219,9 +231,9 @@ We definiëren een functie die volgende stappen uitvoert:
 -   (optioneel) We verwijderen taxa die minder dan `frequency_threshold` (instelbaar) voorkomen in de stalen
     -   voordeel: dit zijn taxa die weinig informatie dragen
     -   nadeel: indien zeldzame taxa geclusterd zitten in een groep, missen we dit patroon als we ze verwijderen
--   (optioneel) "rarefy" de samples tot een gelijk aantal reads (instelbaar)
+-   (optioneel) `"rarefy"` de samples tot een gelijk aantal reads (instelbaar)
     -   voordeel: eventueel vertekening door verschillen in totaal aantal reads tussen stalen wordt weggewerkt; kan zeker van belang zijn als er systematische verschillen in totaal aantal reads zitten tussen de groepen die we willen vergelijken (zie vorige sectie)
-    -   nadeel: dit is een randomisatieprocedure die in feite 10 - 100-tal keer zou herhaald moeten worden waarna de resultaten van elk van de random datasets gecombineerd worden (vergelijkbaar met "multiple imputation"); bij een ordinatie is er echter geen manier om tot een consensus ordinatie te komen
+    -   nadeel: dit is een randomisatieprocedure die in feite 10 - 100-tal keer zou herhaald moeten worden waarna de resultaten van elk van de random datasets gecombineerd worden (vergelijkbaar met `"multiple imputation"`); bij een ordinatie is er echter geen manier om tot een consensus ordinatie te komen
 -   We transformeren met de robuuste "centered-log-ratio" transformatie omdat we werken met compositionele data waar enkel relatieve verschillen in abundantie betekenisvol zijn
 -   Eerst voeren we een principale componenten analyse uit op deze data
     -   de afstanden tussen samples in de ordinatieruimte zijn dan Euclidische afstanden tussen "centered-log-ratio" getransformeerde variabelen - dit is in de literatuur ook bekend als de Aitchinson afstand.
@@ -235,8 +247,7 @@ calc_ordinations <- function(
     rarefy_even_reads = NULL,
     frequency_threshold = 0,
     cumulative_prop_explained = 0.9,
-    use_rank = NULL
-) {
+    use_rank = NULL) {
   set.seed(20315)
 
   data <- tacos_data %>%
@@ -248,9 +259,10 @@ calc_ordinations <- function(
       }
     } %>%
     tidytacos::filter_counts(
-      .by = "taxon_id", n() > frequency_threshold)
+      .by = "taxon_id", n() > frequency_threshold
+    )
+
 
-  
   if (!is.null(rarefy_even_reads)) {
     message(paste("samples rarefied to", rarefy_even_reads, "reads"))
     data <- data %>%
@@ -264,17 +276,21 @@ calc_ordinations <- function(
       id_cols = "sample_id",
       names_from = "taxon_id",
       values_from = "count",
-      values_fill = 0)
+      values_fill = 0
+    )
 
   pca1 <- prcomp(
     vegan::decostand(
       counts %>% select(-"sample_id"),
-      method = "rclr"))
+      method = "rclr"
+    )
+  )
   d <- pca1$sdev^2
   truncate_axis <- sum(cumsum(d) / sum(d) < cumulative_prop_explained)
   message(sprintf(
     "Calculating tsne based on first %s axes from PCA solution",
-    truncate_axis))
+    truncate_axis
+  ))
   tsne1 <- tsne::tsne(pca1$x[, 1:truncate_axis])
 
   plotdata <-
@@ -282,7 +298,8 @@ calc_ordinations <- function(
       as_tibble(vegan::scores(pca1, display = "sites", choices = 1:2)),
       as_tibble(tsne1),
       counts %>%
-        dplyr::select(.data$sample_id)) %>%
+        dplyr::select(.data$sample_id)
+    ) %>%
     left_join(data$samples, by = "sample_id")
 
   return(list(pca = pca1, tsne = tsne1, plotdata = plotdata))
@@ -292,30 +309,39 @@ plot_pca <- function(plotdata, colour = NULL) {
   ggpca <- plotdata %>%
     ggplot() +
     geom_point(
-      aes(x = .data$PC1, y = .data$PC2, colour = {{colour}})
+      aes(x = .data$PC1, y = .data$PC2, colour = {{ colour }})
     ) +
     geom_vline(xintercept = 0, alpha = 0.3) +
     geom_hline(yintercept = 0, alpha = 0.3) +
-    stat_ellipse(aes(
-      x = .data$PC1, y = .data$PC2, colour = {{colour}}),
-      level = 0.95) +
-    labs(x = "Eerste principale component",
-         y = "Tweede principale component")
-    return(ggpca)
+    stat_ellipse(
+      aes(
+        x = .data$PC1, y = .data$PC2, colour = {{ colour }}
+      ),
+      level = 0.95
+    ) +
+    labs(
+      x = "Eerste principale component",
+      y = "Tweede principale component"
+    )
+  return(ggpca)
 }
 
 plot_tsne <- function(plotdata, colour = NULL) {
   ggtsne <- plotdata %>%
     ggplot() +
     geom_point(
-      aes(x = .data$V1, y = .data$V2, colour = {{colour}})) +
+      aes(x = .data$V1, y = .data$V2, colour = {{ colour }})
+    ) +
     geom_vline(xintercept = 0, alpha = 0.3) +
     geom_hline(yintercept = 0, alpha = 0.3) +
     stat_ellipse(aes(
-      x = .data$V1, y = .data$V2, colour = {{colour}}), level = 0.95) +
-    labs(x = "Eerste t-SNE as",
-         y = "Tweede t-SNE as")
-    return(ggtsne)
+      x = .data$V1, y = .data$V2, colour = {{ colour }}
+    ), level = 0.95) +
+    labs(
+      x = "Eerste t-SNE as",
+      y = "Tweede t-SNE as"
+    )
+  return(ggtsne)
 }
 ```
 
@@ -325,7 +351,8 @@ plot_tsne <- function(plotdata, colour = NULL) {
 ordination <- calc_ordinations(
   tacos_data = tidy_physeq,
   frequency_threshold = 0,
-  rarefy_even_reads = 1e+05)
+  rarefy_even_reads = 1e+05
+)
 ```
 
 ```{r}
@@ -356,7 +383,8 @@ ordination <- calc_ordinations(
   tacos_data = tidy_physeq,
   use_rank = "genus",
   frequency_threshold = 0,
-  rarefy_even_reads = 1e+05)
+  rarefy_even_reads = 1e+05
+)
 ```
 
 ```{r}
@@ -398,8 +426,10 @@ sccomp_data_collembola <- tidy_physeq %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, genus, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, genus, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -415,9 +445,9 @@ if (!file.exists(here("source", "brms_models", "s_collembola_genus.rds"))) {
     .data = sccomp_data_collembola,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -425,11 +455,15 @@ if (!file.exists(here("source", "brms_models", "s_collembola_genus.rds"))) {
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_collembola_genus,
-          here("source", "brms_models", "s_collembola_genus.rds"))
+  saveRDS(
+    s_collembola_genus,
+    here("source", "brms_models", "s_collembola_genus.rds")
+  )
 }
 
-s_collembola_genus <- readRDS(here("source", "brms_models", "s_collembola_genus.rds"))
+s_collembola_genus <- readRDS(
+  here("source", "brms_models", "s_collembola_genus.rds")
+)
 ```
 
 ### Model evaluatie
@@ -498,25 +532,30 @@ s_collembola_predictions <- sccomp::sccomp_predict(
   s_collembola_genus,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_collembola_genus %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, genus),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, genus),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_collembola_predictions <- s_collembola_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, genus),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, genus),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     genus = factor(genus, levels = sortering$genus)
@@ -535,11 +574,12 @@ s_collembola_predictions %>%
       y = proportion_mean,
       ymin = proportion_lower,
       ymax = proportion_upper,
-      colour = Diepte),
+      colour = Diepte
+    ),
     position = position_dodge(width = 0.4),
     alpha = 0.7
   ) +
-  facet_wrap(~ genus, scales = "free_y") +
+  facet_wrap(~genus, scales = "free_y") +
   scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -552,9 +592,12 @@ We zien een duidelijke toename van _Megalothorax_, en een duidelijke afname van
 s_collembola_predictions %>%
   ggplot() +
   geom_col(
-    aes(x = Landgebruik_MBAG,
-        y = proportion_mean,
-        fill = genus)) +
+    aes(
+      x = Landgebruik_MBAG,
+      y = proportion_mean,
+      fill = genus
+    )
+  ) +
   facet_grid(~Diepte) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -580,8 +623,10 @@ sccomp_data_collembola <- tidy_physeq_known_genera %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, genus, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, genus, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -592,14 +637,16 @@ Dit is goed voor `r round(sum(sccomp_data_collembola$count) / sum(tidy_physeq_kn
 
 ```{r sccomp-Collembola-known-genus}
 fs::dir_create(here("source", "brms_models"))
-if (!file.exists(here("source", "brms_models", "s_collembola_known_genera.rds"))) {
+if (
+  !file.exists(here("source", "brms_models", "s_collembola_known_genera.rds"))
+) {
   s_collembola_known_genera <- sccomp_estimate(
     .data = sccomp_data_collembola,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -607,11 +654,15 @@ if (!file.exists(here("source", "brms_models", "s_collembola_known_genera.rds"))
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_collembola_known_genera,
-          here("source", "brms_models", "s_collembola_known_genera.rds"))
+  saveRDS(
+    s_collembola_known_genera,
+    here("source", "brms_models", "s_collembola_known_genera.rds")
+  )
 }
 
-s_collembola_known_genera <- readRDS(here("source", "brms_models", "s_collembola_known_genera.rds"))
+s_collembola_known_genera <- readRDS(
+  here("source", "brms_models", "s_collembola_known_genera.rds")
+)
 ```
 
 ### Model evaluatie
@@ -679,25 +730,30 @@ s_collembola_predictions <- sccomp::sccomp_predict(
   s_collembola_known_genera,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_collembola_known_genera %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, genus),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, genus),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_collembola_predictions <- s_collembola_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, genus),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, genus),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     genus = factor(genus, levels = sortering$genus)
@@ -716,11 +772,12 @@ s_collembola_predictions %>%
       y = proportion_mean,
       ymin = proportion_lower,
       ymax = proportion_upper,
-      colour = Diepte),
+      colour = Diepte
+    ),
     position = position_dodge(width = 0.4),
     alpha = 0.7
   ) +
-  facet_wrap(~ genus, scales = "free_y") +
+  facet_wrap(~genus, scales = "free_y") +
   scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -733,9 +790,12 @@ We zien een duidelijke toename van _Megalothorax_, _Protaphorura_, en een duidel
 s_collembola_predictions %>%
   ggplot() +
   geom_col(
-    aes(x = Landgebruik_MBAG,
-        y = proportion_mean,
-        fill = genus)) +
+    aes(
+      x = Landgebruik_MBAG,
+      y = proportion_mean,
+      fill = genus
+    )
+  ) +
   facet_grid(~Diepte) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -761,8 +821,10 @@ sccomp_data_collembola <- tidy_physeq %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, species, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, species, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -778,9 +840,9 @@ if (!file.exists(here("source", "brms_models", "s_collembola_species.rds"))) {
     .data = sccomp_data_collembola,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -788,11 +850,15 @@ if (!file.exists(here("source", "brms_models", "s_collembola_species.rds"))) {
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_collembola_species,
-          here("source", "brms_models", "s_collembola_species.rds"))
+  saveRDS(
+    s_collembola_species,
+    here("source", "brms_models", "s_collembola_species.rds")
+  )
 }
 
-s_collembola_species <- readRDS(here("source", "brms_models", "s_collembola_species.rds"))
+s_collembola_species <- readRDS(
+  here("source", "brms_models", "s_collembola_species.rds")
+)
 ```
 
 ### Model evaluatie
@@ -819,8 +885,8 @@ Als de blauwe boxen goed overeenkomen met de geobserveerde boxen en data, dan is
 Deze boxplot figuren leren ons dat het model de data goed kan beschrijven, maar dat er mogelijkheden zijn om het model te verbeteren:
 
 -   er zijn mogelijk verschillen in variantie tussen de groepen: we kunnen dus eventueel de variantie ook modelleren in functie van deze factoren
--   Sprake van bimodaliteit? In 't292' ( _Megalothorax_ _incertus_) akker? Maar wel verdwenen voor _Megalothorax_ ( _incertus_) dus in natuurgraslanden?: de boxplot zit in de zone met proporties tussen 0 en 0.25, terwijl er een groep van punten is die clustert tussen 0.25 en 1, en vooral tusen 0.5 en 1. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen)?
--   Idem voor _Folsomia_ _candida_ 't298' in tijdelijk grasland? de boxplot zit in de zone met proporties tussen x en x, terwijl er een groep van punten is die clustert tussen y en y
+-   Sprake van bimodaliteit? In `'t292'` (`_Megalothorax_ _incertus_`) akker? Maar wel verdwenen voor `_Megalothorax_` (`_incertus_`) dus in natuurgraslanden?: de boxplot zit in de zone met proporties tussen 0 en 0.25, terwijl er een groep van punten is die clustert tussen 0.25 en 1, en vooral tussen 0.5 en 1. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen)?
+-   Idem voor `_Folsomia_ _candida_` `'t298'` in tijdelijk grasland? de boxplot zit in de zone met proporties tussen x en x, terwijl er een groep van punten is die clustert tussen y en y
 
 ```{r sccomp-boxplot-landgebruik-species}
 plots$boxplot[[1]]
@@ -861,25 +927,30 @@ s_collembola_predictions <- sccomp::sccomp_predict(
   s_collembola_species,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_collembola_species %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_collembola_predictions <- s_collembola_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     species = factor(species, levels = sortering$species)
@@ -887,8 +958,8 @@ s_collembola_predictions <- s_collembola_predictions %>%
 ```
 
 De volgende figuur toont de modelvoorspellingen.
-Duidelijk geen toename van _Megalothorax_ _incertus_ met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. De voorheen geobserveerde toename van _Megalothorax_ met afnemende landgebruiksintensiteit moet dus aan andere _Megalothorax_ soorten te wijten zijn? Duidelijke afname van _Folsomia_ _candida_ en _Isotomurus_ _maculatus_ (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit
-_Megalothorax_ _incertus_ komt duidelijk meer voor in 10-30 cm dan 0-10 cm.
+Duidelijk geen toename van `_Megalothorax_ _incertus_` met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. De voorheen geobserveerde toename van `_Megalothorax_` met afnemende landgebruiksintensiteit moet dus aan andere `_Megalothorax_` soorten te wijten zijn? Duidelijke afname van `_Folsomia_ _candida_` en `_Isotomurus_ _maculatus_` (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit
+`_Megalothorax_ _incertus_` komt duidelijk meer voor in 10-30 cm dan 0-10 cm.
 
 ```{r plot-sccomp-predictions-species}
 s_collembola_predictions %>%
@@ -899,11 +970,12 @@ s_collembola_predictions %>%
       y = proportion_mean,
       ymin = proportion_lower,
       ymax = proportion_upper,
-      colour = Diepte),
+      colour = Diepte
+    ),
     position = position_dodge(width = 0.4),
     alpha = 0.7
   ) +
-  facet_wrap(~ species, scales = "free_y") +
+  facet_wrap(~species, scales = "free_y") +
   scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -914,9 +986,12 @@ We kunnen dit (zonder credibiliteitsinterval) ook als volgt voorstellen:
 s_collembola_predictions %>%
   ggplot() +
   geom_col(
-    aes(x = Landgebruik_MBAG,
-        y = proportion_mean,
-        fill = species)) +
+    aes(
+      x = Landgebruik_MBAG,
+      y = proportion_mean,
+      fill = species
+    )
+  ) +
   facet_grid(~Diepte) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -943,8 +1018,10 @@ sccomp_data_collembola <- tidy_physeq_known_species %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, species, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, species, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -955,14 +1032,16 @@ Dit is goed voor `r round(sum(sccomp_data_collembola$count) / sum(tidy_physeq_kn
 
 ```{r sccomp-Collembola-known-species}
 fs::dir_create(here("source", "brms_models"))
-if (!file.exists(here("source", "brms_models", "s_collembola_known_species.rds"))) {
+if (
+  !file.exists(here("source", "brms_models", "s_collembola_known_species.rds"))
+) {
   s_collembola_known_species <- sccomp_estimate(
     .data = sccomp_data_collembola,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -970,11 +1049,15 @@ if (!file.exists(here("source", "brms_models", "s_collembola_known_species.rds")
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_collembola_known_species,
-          here("source", "brms_models", "s_collembola_known_species.rds"))
+  saveRDS(
+    s_collembola_known_species,
+    here("source", "brms_models", "s_collembola_known_species.rds")
+  )
 }
 
-s_collembola_known_species <- readRDS(here("source", "brms_models", "s_collembola_known_species.rds"))
+s_collembola_known_species <- readRDS(
+  here("source", "brms_models", "s_collembola_known_species.rds")
+)
 ```
 
 ### Model evaluatie
@@ -1001,7 +1084,7 @@ Als de blauwe boxen goed overeenkomen met de geobserveerde boxen en data, dan is
 Deze boxplot figuren leren ons dat het model de data goed kan beschrijven, maar dat er mogelijkheden zijn om het model te verbeteren:
 
 -   er zijn mogelijk verschillen in variantie tussen de groepen: we kunnen dus eventueel de variantie ook modelleren in functie van deze factoren
--   Op het eerste zicht geen sprake van bimodaliteit??? : de boxplot zit in de zone met proporties tussen x en x.x, terwijl er een groep van punten is die clustert tussen y.y en y. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen).
+-   Op het eerste zicht geen sprake van bimodaliteit??? : de boxplot zit in de zone met proporties tussen `x` en `x.x`, terwijl er een groep van punten is die clustert tussen `y.y` en `y`. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen).
 
 ```{r sccomp-boxplot-landgebruik-known-species}
 plots$boxplot[[1]]
@@ -1042,25 +1125,30 @@ s_collembola_predictions <- sccomp::sccomp_predict(
   s_collembola_known_species,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_collembola_known_species %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_collembola_predictions <- s_collembola_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     species = factor(species, levels = sortering$species)
@@ -1068,8 +1156,8 @@ s_collembola_predictions <- s_collembola_predictions %>%
 ```
 
 De volgende figuur toont de modelvoorspellingen.
-Duidelijk geen toename van _Megalothorax_ _incertus_ met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. Wel een duidelijke toename van _Megalothorax_ _minimus_, _Isotoma_ _viridis_, _Protaphorura_ _armata_, en _Megalothorax_ _willemi_ (deze laatste enkel in de 10-30) met afnemende landgebruiksintensiteit. Duidelijke afname van _Folsomia_ _candida_ en _Isotomurus_ _maculatus_ (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit.
-_Megalothorax_ _incertus_, _Megalothorax_ _willemi_, en _Folsomia_ _candida_ lijken meer voor te komen in 10-30 cm dan 0-10 cm, terwijl _Willowsia_ _nigromaculata_, _Isotomurus_ _maculatus_, en _Sminthurus_ _viridis_ meer voor lijken te komen in 0-10 cm, dan 10-30 cm.
+Duidelijk geen toename van `_Megalothorax_ _incertus_` met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. Wel een duidelijke toename van `_Megalothorax_ _minimus_`, `_Isotoma_ _viridis_`, `_Protaphorura_ _armata_`, en `_Megalothorax_ _willemi_` (deze laatste enkel in de 10-30) met afnemende landgebruiksintensiteit. Duidelijke afname van `_Folsomia_ _candida_` en `_Isotomurus_ _maculatus_` (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit.
+`_Megalothorax_ _incertus_`, `_Megalothorax_ _willemi_`, en `_Folsomia_ _candida_` lijken meer voor te komen in 10-30 cm dan 0-10 cm, terwijl `_Willowsia_ _nigromaculata_`, `_Isotomurus_ _maculatus_`, en `_Sminthurus_ _viridis_` meer voor lijken te komen in 0-10 cm, dan 10-30 cm.
 
 ```{r plot-sccomp-predictions-known-species}
 s_collembola_predictions %>%
@@ -1080,11 +1168,12 @@ s_collembola_predictions %>%
       y = proportion_mean,
       ymin = proportion_lower,
       ymax = proportion_upper,
-      colour = Diepte),
+      colour = Diepte
+    ),
     position = position_dodge(width = 0.4),
     alpha = 0.7
   ) +
-  facet_wrap(~ species, scales = "free_y") +
+  facet_wrap(~species, scales = "free_y") +
   scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
@@ -1095,15 +1184,18 @@ We kunnen dit (zonder credibiliteitsinterval) ook als volgt voorstellen:
 s_collembola_predictions %>%
   ggplot() +
   geom_col(
-    aes(x = Landgebruik_MBAG,
-        y = proportion_mean,
-        fill = species)) +
+    aes(
+      x = Landgebruik_MBAG,
+      y = proportion_mean,
+      fill = species
+    )
+  ) +
   facet_grid(~Diepte) +
   theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
 ```
 
 
-## `sccomp`: differentiële samenstelling en variabiliteitsanalyse op soortniveau met alle OTUs met een minimum occurence van 5
+## `sccomp`: differentiële samenstelling en variabiliteitsanalyse op soortniveau met alle OTUs met een minimum voorkomen van 5
 
 ```{r sccomp-inputs-all-species}
 max_taxa <- 207
@@ -1120,8 +1212,10 @@ sccomp_data_collembola <- tidy_physeq %>%
   ) %>%
   tidytacos::everything() %>%
   mutate(count = as.integer(count)) %>%
-  select(taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
-         Landgebruik_MBAG, species, occurrence) %>%
+  select(
+    taxon_id, sample_id, count, sample, Cmon_PlotID, Diepte,
+    Landgebruik_MBAG, species, occurrence
+  ) %>%
   filter(
     complete.cases(.)
   )
@@ -1130,18 +1224,20 @@ sccomp_data_collembola <- tidy_physeq %>%
 Voor deze analyse, aggregeren we de taxa tot `r used_rank` niveau en beperken we de analyse tot de `r max_taxa` taxa met hoogste prevalentie.
 Dit is goed voor `r round(sum(sccomp_data_collembola$count) / sum(tidy_physeq$counts$count) * 100)` % van alle reads.
 
-Hier dus eens getest met alle OTUs met een occurrence van minstens 5 (in dit geval de 207 meest prevalente OTUs)
+Hier dus eens getest met alle OTUs met een voorkomen van minstens 5 (in dit geval de 207 meest voorkomende OTUs)
 
 ```{r sccomp-collembola-all-species, eval=FALSE}
 fs::dir_create(here("source", "brms_models"))
-if (!file.exists(here("source", "brms_models", "s_collembola_207_species.rds"))) {
+if (
+  !file.exists(here("source", "brms_models", "s_collembola_207_species.rds"))
+) {
   s_collembola_207_species <- sccomp_estimate(
     .data = sccomp_data_collembola,
     formula_composition = ~
       Landgebruik_MBAG
-    + Diepte
-    + Landgebruik_MBAG:Diepte
-    + (1 | Cmon_PlotID),
+      + Diepte
+        + Landgebruik_MBAG:Diepte
+        + (1 | Cmon_PlotID),
     .sample = sample_id,
     .cell_group = taxon_id,
     .count = count,
@@ -1150,15 +1246,17 @@ if (!file.exists(here("source", "brms_models", "s_collembola_207_species.rds")))
     max_sampling_iterations = 2000,
     bimodal_mean_variability_association = FALSE
   )
-  saveRDS(s_collembola_207_species,
-          here("source", "brms_models", "s_collembola_207_species.rds"))
+  saveRDS(
+    s_collembola_207_species,
+    here("source", "brms_models", "s_collembola_207_species.rds")
+  )
 }
-
-
 ```
 
 ```{r sccomp-collembola-all-species-load, eval=file.exists("here("source", "brms_models", "s_collembola_207_species.rds")")}
-s_collembola_207_species <- readRDS(here("source", "brms_models", "s_collembola_207_species.rds"))
+s_collembola_207_species <- readRDS(
+  here("source", "brms_models", "s_collembola_207_species.rds")
+)
 ```
 
 
@@ -1177,25 +1275,30 @@ s_collembola_predictions <- sccomp::sccomp_predict(
   s_collembola_207_species,
   formula_composition = ~
     Landgebruik_MBAG +
-    Diepte +
-    Landgebruik_MBAG:Diepte,
+      Diepte +
+      Landgebruik_MBAG:Diepte,
   new_data = my_new_data,
   number_of_draws = 2000,
-  mcmc_seed = 1254)
+  mcmc_seed = 1254
+)
 
 sortering <- s_collembola_207_species %>%
   filter(parameter == "Landgebruik_MBAGNatuurgrasland") %>%
   select(taxon_id, c_lower) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   arrange(c_lower)
 
 s_collembola_predictions <- s_collembola_predictions %>%
   left_join(my_new_data, by = join_by(sample_id)) %>%
-  left_join(sccomp_data_collembola %>%
-              distinct(taxon_id, species),
-            by = join_by(taxon_id)) %>%
+  left_join(
+    sccomp_data_collembola %>%
+      distinct(taxon_id, species),
+    by = join_by(taxon_id)
+  ) %>%
   mutate(
     taxon_id = factor(taxon_id, levels = sortering$taxon_id),
     species = factor(species, levels = sortering$species)
@@ -1203,8 +1306,8 @@ s_collembola_predictions <- s_collembola_predictions %>%
 ```
 
 De volgende figuur toont de modelvoorspellingen.
-Duidelijk geen toename van _xx_ _xx_ met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. De voorheen geobserveerde toename van _xx_ met afnemende landgebruiksintensiteit moet dus aan andere _xx_ soorten te wijten zijn? Duidelijke afname van _xx_ _xx_ en _xx_ _xx_ (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit
-_xx_ _xx_ komt duidelijk meer voor in 10-30 cm dan 0-10 cm.
+Duidelijk geen toename van `_xx_ _xx_` met afnemende landgebruiksintensiteit. Integendeel, duidelijk lager in natuurgrasland. De voorheen geobserveerde toename van `_xx_` met afnemende landgebruiksintensiteit moet dus aan andere `_xx_` soorten te wijten zijn? Duidelijke afname van `_xx_ _xx_` en `_xx_ _xx_` (deze laatste enkel op basis van de 0-10 cm) met afnemende landgebruiksintensiteit
+`_xx_ _xx_` komt duidelijk meer voor in 10-30 cm dan 0-10 cm.
 
 ```{r plot-sccomp-predictions-all-species, eval=file.exists("here("source", "brms_models", "s_collembola_207_species.rds")")}
 # Filter out species with "_otu" in their names
@@ -1213,15 +1316,15 @@ filtered_species_data <- s_collembola_predictions %>%
 
 # Find unique species from the filtered data and split them into chunks of 9
 unique_species <- unique(filtered_species_data$species)
-species_chunks <- split(unique_species, ceiling(seq_along(unique_species)/9))
+species_chunks <- split(unique_species, ceiling(seq_along(unique_species) / 9))
 
 # Loop through each chunk and create a plot
 plots <- list() # To store the plots
 for (i in seq_along(species_chunks)) {
   # Subset the data for the current chunk of species
-  subset_data <- filtered_species_data %>% 
+  subset_data <- filtered_species_data %>%
     filter(species %in% species_chunks[[i]])
-  
+
   # Create the plot for the current subset
   p <- subset_data %>%
     ggplot() +
@@ -1231,22 +1334,20 @@ for (i in seq_along(species_chunks)) {
         y = proportion_mean,
         ymin = proportion_lower,
         ymax = proportion_upper,
-        colour = Diepte),
+        colour = Diepte
+      ),
       position = position_dodge(width = 0.4),
       alpha = 0.7
     ) +
-    facet_wrap(~ species, scales = "free_y", ncol = 3) +
+    facet_wrap(~species, scales = "free_y", ncol = 3) +
     scale_y_log10(breaks = c(0.001, 0.01, 0.05, 0.1, 0.25, 0.5, 1)) +
     theme(axis.text.x = element_text(angle = -90, vjust = 0.5, hjust = 0))
-  
+
   # Save the plot in the list
   plots[[i]] <- p
 }
 
 # Now you can print or save each plot separately
-# For example, to view the first set of plots, you can do:
-# print(plots[[1]])
-# And so on for each plot in 'plots'
 ```
 
 ```{r print-all-plots, eval=file.exists("here("source", "brms_models", "s_collembola_207_species.rds")")}
@@ -1279,8 +1380,8 @@ Als de blauwe boxen goed overeenkomen met de geobserveerde boxen en data, dan is
 Deze boxplot figuren leren ons dat het model de data goed kan beschrijven, maar dat er mogelijkheden zijn om het model te verbeteren:
 
 -   er zijn mogelijk verschillen in variantie tussen de groepen: we kunnen dus eventueel de variantie ook modelleren in functie van deze factoren
--   Sprake van bimodaliteit? In 't292' ( _xx_ _xx_) akker? Maar wel verdwenen voor _xx_ ( _xx_) dus in natuurgraslanden?: de boxplot zit in de zone met proporties tussen 0 en 0.25, terwijl er een groep van punten is die clustert tussen 0.25 en 1, en vooral tusen 0.5 en 1. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen)?
--   Idem voor _xx_ _xx_ 't298' in tijdelijk grasland? de boxplot zit in de zone met proporties tussen x en x, terwijl er een groep van punten is die clustert tussen y en y
+-   Sprake van bimodaliteit? In `'t292'` (`_xx_ _xx_`) akker? Maar wel verdwenen voor `_xx_` (`_xx_`) dus in natuurgraslanden?: de boxplot zit in de zone met proporties tussen 0 en 0.25, terwijl er een groep van punten is die clustert tussen 0.25 en 1, en vooral tussen 0.5 en 1. Na te gaan of deze outliers kunnen verklaard worden met een ontbrekende covariaat (indien niet, kunnen we ze eventueel uit de analyse verwijderen)?
+-   Idem voor `_xx_ _xx_` `'t298'` in tijdelijk grasland? de boxplot zit in de zone met proporties tussen x en x, terwijl er een groep van punten is die clustert tussen y en y
 
 ```{r sccomp-boxplot-landgebruik-all-species, eval=file.exists("here("source", "brms_models", "s_collembola_207_species.rds")")}
 plots$boxplot[[1]]
diff --git a/source/rmarkdown/Collembola_data_analyse.Rmd b/source/rmarkdown/collembola_data_analyse.Rmd
similarity index 79%
rename from source/rmarkdown/Collembola_data_analyse.Rmd
rename to source/rmarkdown/collembola_data_analyse.Rmd
index 48cb582..c56df5b 100644
--- a/source/rmarkdown/Collembola_data_analyse.Rmd
+++ b/source/rmarkdown/collembola_data_analyse.Rmd
@@ -1,5 +1,5 @@
 ---
-title: "Collembola data analyse"
+title: "`Collembola` data analyse"
 author: "Sam Lambrechts, Io Deflem, Emma Cartuyvels, Hans Van Calster"
 date: "`r Sys.Date()`"
 output:
@@ -25,15 +25,17 @@ library(phyloseq)
 library(vegan)
 library(ggplot2)
 library(ggrepel)
-library(ggpubr) #installatie was nodig op hpc
+library(ggpubr) # installatie was nodig op hpc
 library(gridExtra)
-library(ggbiplot) #installatie was nodig op hpc
-library(factoextra) #installatie was nodig op hpc
+library(ggbiplot) # installatie was nodig op hpc
+library(factoextra) # installatie was nodig op hpc
 library(dplyr)
 library(rstatix)
 library(dunn.test) # installatie was nodig op hpc
 library(compositions) # installatie was nodig op hpc
-library(sf) # installatie nodig op hpc, alsook van de dependcy "units", voor units was ook installatie texinfo nodig, installatie units en dus ook sf niet gelukt op hpc
+library(sf) # installatie nodig op hpc, alsook van de dependcy "units",
+# voor units was ook installatie texinfo nodig,
+# installatie units en dus ook sf niet gelukt op hpc
 
 #### Handig pakket met veel phyloseq add-ons
 if (!"psadd" %in% rownames(installed.packages())) {
@@ -50,18 +52,21 @@ if (!"ggVennDiagram" %in% rownames(installed.packages())) {
 
 if (!"microViz" %in% rownames(installed.packages())) {
   install.packages(
-  "microViz",
-  repos = c(davidbarnett = "https://david-barnett.r-universe.dev",
-            getOption("repos"))
-)
+    "microViz",
+    repos = c(
+      davidbarnett = "https://david-barnett.r-universe.dev",
+      getOption("repos")
+    )
+  )
 }
 
-# voor microviz was het ook nodig om de microbiome en ComplexHeatmap packages te installeren via biocmanager op hpc
+# voor microviz was het ook nodig om de microbiome en ComplexHeatmap packages te
+# installeren via biocmanager op hpc
 
-library(tidytacos) #installatie was nodig op hpc
-library(glmmTMB) #installatie was nodig op hpc
-library(marginaleffects) #installatie was nodig op hpc
-library(performance) #installatie was nodig op hpc
+library(tidytacos) # installatie was nodig op hpc
+library(glmmTMB) # installatie was nodig op hpc
+library(marginaleffects) # installatie was nodig op hpc
+library(performance) # installatie was nodig op hpc
 
 mbag_folder <- "G:/Gedeelde drives/PRJ_MBAG" # nolint
 ```
@@ -76,15 +81,15 @@ Deze tabelletjes kunnen best gewoon csv bestanden zijn die we inlezen.
 
 ```{r load-rdata-file}
 path_naar_bestand <- file.path(
-    mbag_folder,
-    "4c_bodembiodiversiteit",
-    "data",
-    "statistiek",
-    "Coll01",
-    "zonder_otu_legetaxonomie",
-    "phyloseq",
-    "physeq_Coll01_Collembola.Rdata"
-  )
+  mbag_folder,
+  "4c_bodembiodiversiteit",
+  "data",
+  "statistiek",
+  "Coll01",
+  "zonder_otu_legetaxonomie",
+  "phyloseq",
+  "physeq_Coll01_Collembola.Rdata"
+)
 
 load(path_naar_bestand)
 ```
@@ -98,8 +103,12 @@ Het bestand dat wordt ingelezen is `r basename(path_naar_bestand)`.
 ```{r hernoem-object}
 physeq <- physeq_Coll01_Collembola
 
-physeq_known_species <- subset_taxa(physeq, !grepl("_otu", tax_table(physeq)[, "species"]))
-physeq_known_genera <- subset_taxa(physeq, !grepl("_otu", tax_table(physeq)[, "genus"]))
+physeq_known_species <- subset_taxa(
+  physeq, !grepl("_otu", tax_table(physeq)[, "species"])
+)
+physeq_known_genera <- subset_taxa(
+  physeq, !grepl("_otu", tax_table(physeq)[, "genus"])
+)
 ```
 
 ```{r}
@@ -107,20 +116,24 @@ physeq <- physeq %>%
   microViz::ps_mutate(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 
 physeq_known_species <- physeq_known_species %>%
   microViz::ps_mutate(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 
 
@@ -128,10 +141,12 @@ physeq_known_genera <- physeq_known_genera %>%
   microViz::ps_mutate(
     Datum_staalname = lubridate::mdy_hm(Datum_staalname),
     Landgebruik_MBAG = factor(
-    Landgebruik_MBAG,
-    levels = c(
-      "Akker", "Tijdelijk grasland", "Blijvend grasland",
-      "Residentieel grasland", "Natuurgrasland"))
+      Landgebruik_MBAG,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland"
+      )
+    )
   )
 ```
 
@@ -145,10 +160,10 @@ Filter taxongroep:
 ```{r filter-taxongroep}
 physeq <- subset_taxa(
   physeq,
-  class == "Collembola")
+  class == "Collembola"
+)
 
 physeq_rarefied <- rarefy_even_depth(physeq, sample.size = 23033)
-
 ```
 
 
@@ -168,7 +183,6 @@ psotu2veg <- function(physeq) {
   }
   return(as(otu, "matrix"))
 }
-
 ```
 
 
@@ -227,7 +241,7 @@ tidy_physeq$samples %>%
   kable()
 ```
 
-En ook de verdeling van de samples over de belangrijkste design-variabelen voor de subsampled (23033 reads) data, om te zien tot welke landgebruikstypes de 28 stalen die zijn afgevallen (door te weinig reads) behoren:
+En ook de verdeling van de samples over de belangrijkste design-variabelen voor de sub-sampled (23033 reads) data, om te zien tot welke landgebruikstypes de 28 stalen die zijn afgevallen (door te weinig reads) behoren:
 
 ```{r designvars_rar}
 tidy_physeq_rarefied$samples %>%
@@ -244,7 +258,7 @@ tidy_physeq$samples %>%
   kable()
 ```
 
-Aantal bodemlocaties subsampled dataset (23033 reads):
+Aantal bodemlocaties sub-sampled dataset (23033 reads):
 
 ```{r}
 tidy_physeq_rarefied$samples %>%
@@ -269,20 +283,17 @@ tidy_physeq$samples %>%
 ## OTU tabel
 
 ```{r verken-otu-data}
-
-
-glimpse(otu_table(physeq_rarefied) %>% as.data.frame %>% as_tibble())
+glimpse(otu_table(physeq_rarefied) %>% as.data.frame() %>% as_tibble())
 ```
 
 ## Taxonomie tabel
 
 
 ```{r verken-taxonomie-data}
-glimpse(tax_table(physeq_rarefied) %>% as.data.frame %>% as_tibble())
+glimpse(tax_table(physeq_rarefied) %>% as.data.frame() %>% as_tibble())
 ```
 
 ```{r GBIF-check-presence}
-
 species <- as.vector(tax_table(physeq)[, "species"])
 species <- species[!grepl("_otu", species)]
 species <- gsub("_", " ", species)
@@ -290,7 +301,7 @@ species <- species[nchar(species) > 0]
 source(here::here("source/r/check_presence.R"))
 gbif_check <- check_presence(species)
 gbif_check %>%
-filter(!present) %>%
+  filter(!present) %>%
   kable(caption = "Not present according to GBIF in Western-Europe")
 ```
 
@@ -315,7 +326,10 @@ tidy_physeq_rarefied %>%
 
 ```{r krona-plot}
 # Create a new combined column
-combined_column <- paste(sample_data(physeq_rarefied)$Landgebruik_MBAG, sample_data(physeq_rarefied)$Diepte, sep = "_")
+combined_column <- paste(sample_data(physeq_rarefied)$Landgebruik_MBAG,
+  sample_data(physeq_rarefied)$Diepte,
+  sep = "_"
+)
 
 # Add the new combined column to the sample metadata with the desired name
 sample_data(physeq_rarefied)$Landgebruik_MBAG_diepte <- combined_column
@@ -324,8 +338,16 @@ sample_data(physeq_rarefied)$Landgebruik_MBAG_diepte <- combined_column
 physeq_rarefied <- merge_phyloseq(physeq_rarefied, sample_data(physeq_rarefied))
 
 
-if (system(command = "which ktImportText",  intern = FALSE, ignore.stderr = TRUE,  ignore.stdout = TRUE) != 1) {
-  psadd::plot_krona(physeq_rarefied, "MBAG_Olig01_collembola_rar_species_alle_stalen_Landgebruik_MBAG_per_diepte", "Landgebruik_MBAG_diepte", trim = T)
+if (system(
+  command = "which ktImportText", intern = FALSE, ignore.stderr = TRUE,
+  ignore.stdout = TRUE
+) != 1) {
+  psadd::plot_krona(
+    physeq_rarefied,
+    "MBAG_Olig01_collembola_rar_species_alle_stalen_Landgebruik_MBAG_per_diepte"
+    , "Landgebruik_MBAG_diepte",
+    trim = TRUE
+  )
 }
 ```
 
@@ -351,7 +373,8 @@ Het geobserveerde aantal taxa neemt toe met het aantal reads in een staal:
 tidy_physeq$samples %>%
   filter(total_count > 0) %>%
   ggplot(
-    aes(x = total_count, y = observed)) +
+    aes(x = total_count, y = observed)
+  ) +
   geom_point() +
   scale_x_log10() +
   labs(y = "Observed number of taxa")
@@ -371,16 +394,17 @@ samples_totcount_filtered <- tidy_physeq$samples %>%
 form_collembola_richness <- formula(
   observed ~
     log(total_count)
-  + Landgebruik_MBAG
-  + Diepte
-  + Landgebruik_MBAG:Diepte
-  + (1 | Cmon_PlotID)
-  )
+    + Landgebruik_MBAG
+    + Diepte
+    + Landgebruik_MBAG:Diepte
+    + (1 | Cmon_PlotID)
+)
 
 m_collembola_richness <- glmmTMB::glmmTMB(
   formula = form_collembola_richness,
   data = samples_totcount_filtered,
-  family = glmmTMB::nbinom2())
+  family = glmmTMB::nbinom2()
+)
 ```
 
 Zijn er problemen met de residuele variabiliteit op basis van de model validatie?
@@ -406,7 +430,8 @@ Zoals te verwachten, is het belangrijk om te corrigeren voor het totaal aantal r
 marginaleffects::plot_predictions(
   m_collembola_richness,
   condition = c("total_count"),
-  vcov = vcov(m_collembola_richness)) +
+  vcov = vcov(m_collembola_richness)
+) +
   labs(y = "Voorspeld aantal taxa") +
   scale_x_log10()
 ```
@@ -420,21 +445,20 @@ marginaleffects::plot_predictions(
   condition = c("Landgebruik_MBAG", "Diepte"),
   vcov = vcov(m_collembola_richness),
   re.form = NA,
-  type = "response") +
+  type = "response"
+) +
   labs(y = "Voorspelde aantal taxa")
 ```
 
 
-### Shannon diversity index analyse
+### Shannon diversiteit index analyse
 
 Dit is een alternatief voor voorgaande analyse.
 
-Shannon diversity index modelleren als een Negatief Binomiaal model zoals hierboven?
+Shannon diversiteitsindex modelleren als een Negatief Binomiaal model zoals hierboven?
 
 ```{r m-collembola-shannon-diversity-predictions}
 shannon_diversity_data <- estimate_richness(physeq, measures = "Shannon")
-
-
 ```
 
 
@@ -451,16 +475,20 @@ estimated_species_richness <- vegan::rarefy(
   veganobject,
   23033,
   se = TRUE,
-  MARGIN = 1)
+  MARGIN = 1
+)
 
 
 true_rarefaction_results <- cbind(
   sam_new %>% select(Cmon_PlotID, Diepte, Landgebruik_MBAG),
   est = estimated_species_richness["S", ],
-  se = estimated_species_richness["se", ]) %>%
+  se = estimated_species_richness["se", ]
+) %>%
   as_tibble() %>%
-  mutate(est_lower = est - 2 * se,
-         est_upper = est + 2 * se)
+  mutate(
+    est_lower = est - 2 * se,
+    est_upper = est + 2 * se
+  )
 ```
 
 
@@ -476,17 +504,17 @@ library(brms)
 
 form_rarefaction <- bf(
   est | se(se, sigma = TRUE) ~
-  + Landgebruik_MBAG
-  + Diepte
-  + Landgebruik_MBAG:Diepte
-  + (1 | Cmon_PlotID)
+    +Landgebruik_MBAG
+    + Diepte
+    + Landgebruik_MBAG:Diepte
+    + (1 | Cmon_PlotID)
 )
 
 m_rarefaction <- brm(
   formula = form_rarefaction,
   data = true_rarefaction_results,
-  family = gaussian())
-
+  family = gaussian()
+)
 ```
 
 
@@ -500,8 +528,9 @@ Geeft dit model gelijkaardige resultaten als het model waarbij we rechtstreeks w
 ```{r m-rarefied-predictions}
 conditional_effects(
   m_rarefaction,
-    effects = c("Landgebruik_MBAG:Diepte"),
-  re.form = NA)
+  effects = c("Landgebruik_MBAG:Diepte"),
+  re.form = NA
+)
 ```
 
 ## OTU counts 
@@ -578,8 +607,6 @@ samples <- sample_data(physeq_rarefied) %>%
 rda <- rda(ordmat)
 
 biplot(rda)
-
-
 ```
 
 
@@ -595,7 +622,9 @@ pcoa_result <- ordinate(physeq, method = "PCoA", distance = "bray")
 
 # Create the ordination plot with color based on land_use and shape
 # based on diepte
-plot_ordination(physeq, pcoa_result, type = "samples",
-                color = "Landgebruik_MBAG", shape = "Diepte") +
+plot_ordination(physeq, pcoa_result,
+  type = "samples",
+  color = "Landgebruik_MBAG", shape = "Diepte"
+) +
   theme_minimal()
 ```
diff --git a/source/rmarkdown/compositional_analysis/test_rademu.Rmd b/source/rmarkdown/compositional_analysis/test_rademu.Rmd
index 4298697..be3eaab 100644
--- a/source/rmarkdown/compositional_analysis/test_rademu.Rmd
+++ b/source/rmarkdown/compositional_analysis/test_rademu.Rmd
@@ -28,15 +28,16 @@ library(radEmu)
 library(phyloseq)
 ```
 
-# Inlezen data
+# Read data
 
 ```{r inlezen}
 metadata <- readr::read_csv(
   file.path(
     mbag_bodem_folder,
-    "data", 
+    "data",
     "Stratificatie_MBAG_plots",
-    "MBAG_stratfile_v2_cleaned_13.csv")
+    "MBAG_stratfile_v2_cleaned_13.csv"
+  )
 ) %>%
   janitor::clean_names() %>%
   rename(
@@ -52,7 +53,8 @@ metadata <- readr::read_csv(
       landgebruik_mbag,
       levels = c(
         "Akker", "Tijdelijk grasland", "Blijvend grasland",
-        "Residentieel grasland", "Natuurgrasland", "Heide", "Moeras")
+        "Residentieel grasland", "Natuurgrasland", "Heide", "Moeras"
+      )
     ),
     diepte = gsub("_|/", "-", diepte) |> factor()
   )
@@ -61,9 +63,9 @@ load(
     mbag_bodem_folder,
     "data", "statistiek", "Annelida", "phyloseq",
     "physeq_Olig01_Annelida_species.Rdata"
-    )
   )
-physeq_Olig01_Annelida_species <- physeq_Olig01_Annelida_species |> 
+)
+physeq_olig01_annelida_species <- physeq_Olig01_Annelida_species |> # nolint
   phyloseq::subset_samples(
     !Landgebruik_MBAG %in% c("Moeras", "Heide")
   )
@@ -72,26 +74,26 @@ physeq_Olig01_Annelida_species <- physeq_Olig01_Annelida_species |>
 
 
 ```{r}
-sample_data(physeq_Olig01_Annelida_species)$Landgebruik_MBAG <- factor(
-  sample_data(physeq_Olig01_Annelida_species)$Landgebruik_MBAG,
+sample_data(physeq_olig01_annelida_species)$Landgebruik_MBAG <- factor(
+  sample_data(physeq_olig01_annelida_species)$Landgebruik_MBAG,
   levels = c(
     "Akker", "Tijdelijk grasland", "Blijvend grasland",
-    "Residentieel grasland", "Natuurgrasland")
+    "Residentieel grasland", "Natuurgrasland"
+  )
 )
-
 ```
 
 >Next, we want to confirm that all samples have at least one non-zero count across the categories we’ve chosen and that all categories have at least one non-zero count across the samples we’ve chosen.
 
 ```{r}
-sum(rowSums(otu_table(physeq_Olig01_Annelida_species)) == 0)
-sum(colSums(otu_table(physeq_Olig01_Annelida_species)) == 0)
+sum(rowSums(otu_table(physeq_olig01_annelida_species)) == 0)
+sum(colSums(otu_table(physeq_olig01_annelida_species)) == 0)
 ```
 
 ```{r}
 # reduce data size a bit, because emuFit on full dataset takes a long time
 # here aggregating to genus
-physeq_Olig01_Annelida_genus <- physeq_Olig01_Annelida_species |>
+physeq_olig01_annelida_genus <- physeq_olig01_annelida_species |>
   phyloseq::tax_glom(taxrank = "genus")
 ```
 
@@ -99,10 +101,11 @@ physeq_Olig01_Annelida_genus <- physeq_Olig01_Annelida_species |>
 
 ```{r}
 m_fit <- emuFit(
-  formula = ~ Landgebruik_MBAG + Diepte + Landgebruik_MBAG:Diepte, 
-  Y = physeq_Olig01_Annelida_genus,
-  cluster = sample_data(physeq_Olig01_Annelida_genus)$Cmon_PlotID,
-  run_score_tests = FALSE) 
+  formula = ~ Landgebruik_MBAG + Diepte + Landgebruik_MBAG:Diepte,
+  Y = physeq_olig01_annelida_genus,
+  cluster = sample_data(physeq_olig01_annelida_genus)$Cmon_PlotID,
+  run_score_tests = FALSE
+)
 
 m_fit$estimation_converged
 ```
@@ -112,14 +115,15 @@ Visualise:
 ```{r}
 coeftab <- as_tibble(m_fit$coef)
 
-genustab <- phyloseq::tax_table(physeq_Olig01_Annelida_genus)
+genustab <- phyloseq::tax_table(physeq_olig01_annelida_genus)
 
 genustbl <- as_tibble(
   genustab@.Data
-) |> 
+) |>
   mutate(
-    category = dimnames(genustab)[[1]]) |> 
-  rowwise() |> 
+    category = dimnames(genustab)[[1]]
+  ) |>
+  rowwise() |>
   mutate(
     across(
       all_of(c("phylum", "class", "order", "family", "genus")),
@@ -128,9 +132,9 @@ genustbl <- as_tibble(
     ),
     index = sum(c_across(phylum_:genus_)),
     resolved_rank = c("phylum", "class", "order", "family", "genus")[index]
-  ) |> 
-  ungroup() |> 
-  dplyr::select(phylum:category, resolved_rank) |> 
+  ) |>
+  ungroup() |>
+  dplyr::select(phylum:category, resolved_rank) |>
   mutate(
     resolved_name = case_when(
       resolved_rank == "genus" ~ genus,
@@ -149,24 +153,24 @@ genustbl <- as_tibble(
   )
 
 coeftab |>
-  inner_join(genustbl) |> 
+  inner_join(genustbl) |>
   mutate(
     genus = factor(genus, levels = unique(genus[order(resolved_higher)]))
-  ) |> 
+  ) |>
   ggplot() +
   geom_pointrange(
     aes(
       x = genus,
       y = estimate, ymin = lower, ymax = upper,
-      colour = resolved_higher)
-  ) + 
+      colour = resolved_higher
+    )
+  ) +
   facet_wrap(~covariate)
-
 ```
 
 Identify taxa for which you want robust score tests:
 
-For instance, which taxa at depth 0 - 10 are more likely to occur or not occur in Natuurgrasland compared to the reference category Akker:
+For instance, which taxa at depth 0 - 10 are more likely to occur or not occur in `Natuurgrasland` compared to the reference category `Akker`:
 
 The following code chunk takes a very long time to run and is therefore not evaluated.
 But it can be run in parallel: https://statdivlab.github.io/radEmu/articles/parallel_radEmu.html
@@ -178,7 +182,7 @@ This would need to be multiplied by the number of taxa of interest to get the ti
 
 
 ```{r eval=FALSE}
-coefselection <- coeftab |> 
+coefselection <- coeftab |>
   filter(
     covariate == "Landgebruik_MBAGNatuurgrasland",
     sign(lower) == sign(upper)
@@ -186,16 +190,18 @@ coefselection <- coeftab |>
 
 
 covariate_to_test <- which(
-  rownames(m_fit$B) == "Landgebruik_MBAGNatuurgrasland")
+  rownames(m_fit$B) == "Landgebruik_MBAGNatuurgrasland"
+)
 taxa_to_test <- which(
   coefselection$category %in% colnames(m_fit$B)
-  )
+)
 m_refit <- emuFit(
-  formula = ~ Landgebruik_MBAG + Diepte + Landgebruik_MBAG:Diepte, 
-  Y = physeq_Olig01_Annelida_genus,
+  formula = ~ Landgebruik_MBAG + Diepte + Landgebruik_MBAG:Diepte,
+  Y = physeq_olig01_annelida_genus,
   test_kj = data.frame(
-    k = covariate_to_test, 
-    j = taxa_to_test[1]),
+    k = covariate_to_test,
+    j = taxa_to_test[1]
+  ),
   fitted_model = m_fit,
   refit = FALSE,
   run_score_tests = TRUE
@@ -216,13 +222,14 @@ as_tibble(m_refit$coef) |>
   ) %>%
   mutate(
     genus = factor(genus, levels = unique(genus[order(resolved_higher)]))
-  ) |> 
+  ) |>
   ggplot() +
   geom_pointrange(
     aes(
       x = genus,
       y = estimate, ymin = lower, ymax = upper,
-      colour = resolved_higher)
+      colour = resolved_higher
+    )
   ) +
   labs(
     title = "Natuurgrasland vs Akker (diepte 0-10)"
@@ -231,13 +238,13 @@ as_tibble(m_refit$coef) |>
 
 
 
-More details of radEmu approach:
+More details of `radEmu` approach:
 
-- ‘we introduce a Firth penalty on β derived from a formal equivalence between our model and the multinomial logistic model’ (Clausen en Willis, 2024, p. 6)
+- 'we introduce a Firth penalty on $\beta$ derived from a formal equivalence between our model and the multinomial logistic model' (Clausen en Willis, 2024, p. 6)
 
-- ‘we first introduce a Poisson log likelihood for our model.’ (Clausen en Willis, 2024, p. 6)
+- 'we first introduce a Poisson log likelihood for our model.' (Clausen en Willis, 2024, p. 6)
 
-- ‘Note that the profile likelihood (8) is equal to a multinomial log likelihood with a logistic link (up to a constant). This accords with both known results about the marginal Poisson distribution of multinomial random variables (Birch, 1963), and our observations on identifiability of β (only p × (J − 1) parameters can be identified in a multinomial logistic regression of a J-dimensional outcome on p regressors)’ (Clausen en Willis, 2024, p. 6)
+- 'Note that the profile likelihood (8) is equal to a multinomial log likelihood with a logistic link (up to a constant). This accords with both known results about the marginal Poisson distribution of multinomial random variables (Birch, 1963), and our observations on identifiability of $\beta$ (only p × (J − 1) parameters can be identified in a multinomial logistic regression of a J-dimensional outcome on p regressors)' (Clausen en Willis, 2024, p. 6)
 
 - they show that the model is robust against model-misspecification: simulated data generated under a zero-inflated negative binomial (instead of poisson) still had favourable type I error
     - but a drawback is that you cannot rely on the confidence intervals to judge significance: the score test needs to be computed. In case of modelmisspecification, I think the confidence intervals will likely be too narrow (given that data are usually overdispersed compared to the Poisson)
@@ -247,6 +254,6 @@ More details of radEmu approach:
 		- independent beta-binomial models under constraint that sum of probabilities across taxa equals 1
 		- logit-based fold differences
 		- compared to the multinomial logistic model, sccomp can cope directly with excess variance through the beta-binomial (instead of relying on 'robust' procedures to make inferences). This is important for this type of data.
-		- see Table 1 of [Mangiola et al 2023](https://doi.org/10.1073/pnas.2203828120) where the radEmu approach can be fit in (most similar to ANCOM-BC2 according to Clausen en Willis 2024 paper).
+		- see Table 1 of [Mangiola et al 2023](https://doi.org/10.1073/pnas.2203828120) where the `radEmu` approach can be fit in (most similar to `ANCOM-BC2` according to Clausen en Willis 2024 paper).
 
 
diff --git a/source/rmarkdown/metadata_bodemstalen.Rmd b/source/rmarkdown/metadata_bodemstalen.Rmd
index 35bd5ec..e6c14a5 100644
--- a/source/rmarkdown/metadata_bodemstalen.Rmd
+++ b/source/rmarkdown/metadata_bodemstalen.Rmd
@@ -27,7 +27,7 @@ conflicted::conflicts_prefer(dplyr::select, dplyr::filter)
 # Inlezen data
 
 ```{r inlezen}
-mbag_shared_drive <- "G:/Gedeelde drives/PRJ_MBAG"
+mbag_shared_drive <- "G:/Gedeelde drives/PRJ_MBAG" # nolint
 bodem_meta <- read_excel(
   file.path(
     mbag_shared_drive,
@@ -38,8 +38,10 @@ bodem_meta <- read_excel(
   sheet = "Alle_stalen"
 ) %>%
   janitor::clean_names() %>%
-  mutate(mbag_e_dna_staal = factor(mbag_e_dna_staal),
-         mbag_nematodenstaal = factor(mbag_nematodenstaal)) 
+  mutate(
+    mbag_e_dna_staal = factor(mbag_e_dna_staal),
+    mbag_nematodenstaal = factor(mbag_nematodenstaal)
+  )
 
 nematodenstalen <- read_excel(
   file.path(
@@ -51,13 +53,15 @@ nematodenstalen <- read_excel(
   sheet = "Subset_nematodenstalen_ILVO"
 ) %>%
   janitor::clean_names()
-
 ```
 
 
 ```{r inlezen-gis}
-
-filepath <- "Z:/Projects/PRJ_MBAG/4c_bodembiodiversiteit/steekproef/MBAG_eDNA_sampling"
+filepath <-
+  file.path(
+    "Z:",
+    "Projects/PRJ_MBAG/4c_bodembiodiversiteit/steekproef/MBAG_eDNA_sampling"
+  )
 
 bodem_meta_sf <- read_sf(filepath) %>%
   janitor::clean_names() %>%
@@ -83,9 +87,11 @@ Dit komt neer op alle locaties met grasland of akker als landgebruik waarvan er
 bodem_meta %>%
   filter(mbag_e_dna_staal == 1) %>%
   count(staalopslag, diepte, mbag_luc) %>%
-  pivot_wider(values_from = n,
-              names_from = diepte,
-              names_prefix = "diepte_")
+  pivot_wider(
+    values_from = n,
+    names_from = diepte,
+    names_prefix = "diepte_"
+  )
 ```
 
 De missing data bij `mbag_luc`, kunnen we ondervangen via `cmon_lu_text`:
@@ -93,13 +99,16 @@ De missing data bij `mbag_luc`, kunnen we ondervangen via `cmon_lu_text`:
 ```{r}
 bodem_meta %>%
   filter(mbag_e_dna_staal == 1) %>%
-  mutate(mbag_luc = ifelse(is.na(mbag_luc), cmon_l_utext, mbag_luc),
-         mbag_luc = ifelse(mbag_luc == "Akkerland", "Akker", mbag_luc)) %>%
+  mutate(
+    mbag_luc = ifelse(is.na(mbag_luc), cmon_l_utext, mbag_luc),
+    mbag_luc = ifelse(mbag_luc == "Akkerland", "Akker", mbag_luc)
+  ) %>%
   count(staalopslag, diepte, mbag_luc) %>%
-  pivot_wider(values_from = n,
-              names_from = diepte,
-              names_prefix = "diepte_")
-
+  pivot_wider(
+    values_from = n,
+    names_from = diepte,
+    names_prefix = "diepte_"
+  )
 ```
 
 Hoe belangrijk zijn de 0-10 cm versus 10-30 cm stalen voor eDNA?
@@ -115,13 +124,16 @@ glimpse(nematodenstalen)
 
 ```{r}
 bodem_locs_meta_sf <- bodem_meta_sf %>%
-  mutate(mbag_luc = ifelse(is.na(mbag_luc), cmon_l_utext, mbag_luc),
-         mbag_luc = ifelse(mbag_luc == "Akkerland", "Akker", mbag_luc)) %>%
+  mutate(
+    mbag_luc = ifelse(is.na(mbag_luc), cmon_l_utext, mbag_luc),
+    mbag_luc = ifelse(mbag_luc == "Akkerland", "Akker", mbag_luc)
+  ) %>%
   select(plot_id, diepte, starts_with("mbag")) %>%
   group_by(plot_id, mbag_luc, mbag_e_dna) %>%
-  summarise(dieptes = paste(diepte, collapse = " en "),
-            .groups = "drop")
-
+  summarise(
+    dieptes = paste(diepte, collapse = " en "),
+    .groups = "drop"
+  )
 ```
 
 
@@ -134,7 +146,8 @@ Een deel van de plots zijn momenteel nog in analyse, maar binnenkort hebben we w
 
 ```{r landbouwstreken}
 landbouwstreken <- read_sf(
-  "S:/Vlaanderen/Landbouw/Landbouwstreken_België/Lbstrbel.shp") %>%
+  "S:/Vlaanderen/Landbouw/Landbouwstreken_België/Lbstrbel.shp" # nolint
+) %>%
   janitor::clean_names() %>%
   st_transform(crs = 31370)
 
@@ -145,13 +158,17 @@ bodem_locs_meta_sf <- bodem_locs_meta_sf %>%
 ```{r}
 bodem_locs_meta_sf %>%
   ggplot() +
-  geom_sf(data = landbouwstreken %>%
-            st_filter(bodem_locs_meta_sf),
-          aes(fill = naam),
-          alpha = 0.2) +
-  geom_sf_text(data = landbouwstreken %>%
-            st_filter(bodem_locs_meta_sf),
-          aes(label = naam)) +
+  geom_sf(
+    data = landbouwstreken %>%
+      st_filter(bodem_locs_meta_sf),
+    aes(fill = naam),
+    alpha = 0.2
+  ) +
+  geom_sf_text(
+    data = landbouwstreken %>%
+      st_filter(bodem_locs_meta_sf),
+    aes(label = naam)
+  ) +
   geom_sf(aes(colour = mbag_e_dna)) +
   guides(fill = "none")
 ```
@@ -167,11 +184,13 @@ bodem_locs_meta_sf %>%
   count(mbag_e_dna, landbouwstreek, name = "aantal_locaties") %>%
   mutate(mbag_e_dna = factor(
     mbag_e_dna,
-    labels = c("niet geselecteerd", "wel geselecteeerd"))) %>%
+    labels = c("niet geselecteerd", "wel geselecteeerd")
+  )) %>%
   pivot_wider(
     names_from = mbag_e_dna,
     values_from = aantal_locaties,
-    values_fill = 0) %>%
+    values_fill = 0
+  ) %>%
   kable()
 ```
 
@@ -180,13 +199,16 @@ Wanneer we deze data (enkel de wel geselecteerde locaties en zonder Duinen en We
 ```{r}
 bodem_locs_meta_sf %>%
   st_drop_geometry() %>%
-  filter(mbag_e_dna == 1,
-         !landbouwstreek %in% c("Duinen", "Weidestreek (Luik)")) %>%
+  filter(
+    mbag_e_dna == 1,
+    !landbouwstreek %in% c("Duinen", "Weidestreek (Luik)")
+  ) %>%
   count(mbag_e_dna, mbag_luc, landbouwstreek, name = "aantal_locaties") %>%
   pivot_wider(
     names_from = mbag_luc,
     values_from = aantal_locaties,
-    values_fill = 0) %>%
+    values_fill = 0
+  ) %>%
   select(-mbag_e_dna) %>%
   kable()
 ```
@@ -205,23 +227,23 @@ lbg_binding <- arrow::open_dataset(
 )
 
 bodem_locs_lbg <- landusemetrics_grid_cell(
-    grid_cell = bodem_locs_meta_sf %>%
-      st_buffer(dist = 10),
-    layer = lbg_binding %>%
-      select(LBLHFDTLT, geometry) %>%
-      sfarrow::read_sf_dataset() %>%
-      st_transform(31370),
-    grid_group_by_col = "plot_id",
-    layer_group_by_col = "LBLHFDTLT")
+  grid_cell = bodem_locs_meta_sf %>%
+    st_buffer(dist = 10),
+  layer = lbg_binding %>%
+    select(LBLHFDTLT, geometry) %>%
+    sfarrow::read_sf_dataset() %>%
+    st_transform(31370),
+  grid_group_by_col = "plot_id",
+  layer_group_by_col = "LBLHFDTLT"
+)
 
 mapping <- lbg_binding %>%
-    select(GEWASGROEP, LBLHFDTLT) %>%
-    collect() %>%
-    distinct()
+  select(GEWASGROEP, LBLHFDTLT) %>%
+  collect() %>%
+  distinct()
 
 bodem_locs_lbg <- bodem_locs_lbg %>%
-    left_join(mapping)
-
+  left_join(mapping)
 ```
 
 De kolom samenstelling geeft de samenstelling van gewasgroepen rond de eDNA bodemstaalnamelocaties.
@@ -229,12 +251,12 @@ De kolom samenstelling geeft de samenstelling van gewasgroepen rond de eDNA bode
 ```{r landbouwhoofdteelt}
 bodem_locs_meta_sf %>%
   left_join(bodem_locs_lbg %>%
-              group_by(plot_id, GEWASGROEP) %>%
-              summarise(area_prop = sum(area_prop)) %>%
-              summarise(samenstelling = paste(
-                paste0(GEWASGROEP, " (", round(area_prop, 2), ")"),
-                collapse = " - "
-              ))) %>%
+    group_by(plot_id, GEWASGROEP) %>%
+    summarise(area_prop = sum(area_prop)) %>%
+    summarise(samenstelling = paste(
+      paste0(GEWASGROEP, " (", round(area_prop, 2), ")"),
+      collapse = " - "
+    ))) %>%
   filter(mbag_e_dna == 1) %>%
   count(mbag_luc, samenstelling) %>%
   arrange(desc(n)) %>%