From b63052732af7e6dbe7bc1b4a240a4fa282eea9a4 Mon Sep 17 00:00:00 2001
From: hansvancalster <hans.vancalster@inbo.be>
Date: Wed, 10 Jul 2024 15:39:13 +0200
Subject: [PATCH 1/3] add rmarkdowns

---
 .../_child_model_observed_richness.Rmd        |  77 +++++++
 .../_child_model_shannon.Rmd                  |  78 +++++++
 .../_child_model_simpson.Rmd                  |  66 ++++++
 .../analyses_diversity/analyses_diversity.Rmd | 210 ++++++++++++++++++
 4 files changed, 431 insertions(+)
 create mode 100644 source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd
 create mode 100644 source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd
 create mode 100644 source/rmarkdown/analyses_diversity/_child_model_simpson.Rmd
 create mode 100644 source/rmarkdown/analyses_diversity/analyses_diversity.Rmd

diff --git a/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd b/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd
new file mode 100644
index 0000000..4723176
--- /dev/null
+++ b/source/rmarkdown/analyses_diversity/_child_model_observed_richness.Rmd
@@ -0,0 +1,77 @@
+### `r stringr::str_to_sentence("{{group}} - {{primerset}} - {{unit}}")`
+
+```{r}
+data_subset <- combined %>%
+  filter(
+    group == "{{group}}",
+    primerset == "{{primerset}}",
+    unit == "{{unit}}"
+  )
+```
+
+
+```{r m-richness-{{group}}-{{primerset}}-{{unit}}}
+form_richness <- formula(
+  observed ~
+    log(total_count)
+  + landgebruik
+  + diepte
+  + landgebruik:diepte
+  + (1 | cmon_plot_id)
+  )
+
+cat("Fitting Poisson model")
+
+m_richness <- glmmTMB(
+  formula = form_richness,
+  data = data_subset,
+  family = poisson())
+
+test <- check_overdispersion(m_richness)
+test <- test$p_value < 0.05
+
+if (test) {
+  cat("Overdispersion detected: fitting Negative binomial model instead.")
+  m_richness <- glmmTMB(
+    formula = form_richness,
+    data = data_subset,
+    family = nbinom2())
+}
+```
+
+
+```{r m-richness-checks-{{group}}-{{primerset}}-{{unit}}}
+performance::check_model(m_richness)
+```
+
+
+```{r m-richness-summary-{{group}}-{{primerset}}-{{unit}}}
+summary(m_richness) |> print(digits = 2)
+```
+
+```{r m-richness-anova-{{group}}-{{primerset}}-{{unit}}}
+car::Anova(m_richness) |> print(digits = 2)
+```
+
+```{r m-richness-preds-total-count-{{group}}-{{primerset}}-{{unit}}}
+marginaleffects::plot_predictions(
+  m_richness,
+  condition = c("total_count"),
+  re.form = NA,
+  vcov = TRUE,
+  type = "response") +
+  labs(y = "Voorspeld aantal taxa") +
+  scale_x_log10()
+```
+
+```{r m-richness-preds-landgebruikxdiepte-{{group}}-{{primerset}}-{{unit}}}
+marginaleffects::plot_predictions(
+  m_richness,
+  condition = c("landgebruik", "diepte"),
+  re.form = NA,
+  vcov = TRUE,
+  type = "response") +
+  labs(y = "Voorspeld aantal taxa")
+```
+
+
diff --git a/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd b/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd
new file mode 100644
index 0000000..174997d
--- /dev/null
+++ b/source/rmarkdown/analyses_diversity/_child_model_shannon.Rmd
@@ -0,0 +1,78 @@
+### `r stringr::str_to_sentence("{{group}} - {{primerset}} - {{unit}}")`
+
+```{r}
+data_subset <- combined %>%
+  filter(
+    group == "{{group}}",
+    primerset == "{{primerset}}",
+    unit == "{{unit}}"
+  )
+```
+
+
+```{r m-shannon-{{group}}-{{primerset}}-{{unit}}}
+form_shannon <- formula(
+  shannon ~
+    log(total_count)
+  + landgebruik
+  + diepte
+  + landgebruik:diepte
+  + (1 | cmon_plot_id)
+  )
+
+zeroes <- min(data_subset$shannon) == 0
+
+if (!zeroes) {
+  cat("Fitting Gamma model")
+
+  m_shannon <- glmmTMB(
+    formula = form_shannon,
+    data = data_subset,
+    family = Gamma(link = "log"))
+} else {
+  cat("Zeroes detected (only one taxum observed): Fitting zero-inflated Gamma model")
+
+  m_shannon <- glmmTMB(
+    formula = form_shannon,
+    ziformula = ~ 1,
+    data = data_subset,
+    family = ziGamma(link = "log"))
+}
+```
+
+
+```{r m-shannon-checks-{{group}}-{{primerset}}-{{unit}}}
+performance::check_model(m_shannon)
+```
+
+
+```{r m-shannon-summary-{{group}}-{{primerset}}-{{unit}}}
+summary(m_shannon) |> print(digits = 2)
+```
+
+```{r m-shannon-anova-{{group}}-{{primerset}}-{{unit}}}
+car::Anova(m_shannon) |> print(digits = 2)
+```
+
+```{r m-shannon-preds-total-count-{{group}}-{{primerset}}-{{unit}}}
+marginaleffects::plot_predictions(
+  m_shannon,
+  condition = c("total_count"),
+  re.form = NA,
+  vcov = TRUE,
+  type = "response") +
+  labs(y = "Shannon index") +
+  scale_x_log10()
+```
+
+```{r m-shannon-preds-landgebruikxdiepte-{{group}}-{{primerset}}-{{unit}}}
+marginaleffects::plot_predictions(
+  m_shannon,
+  condition = c("landgebruik", "diepte"),
+  re.form = NA,
+  vcov = TRUE,
+  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
new file mode 100644
index 0000000..29c57ef
--- /dev/null
+++ b/source/rmarkdown/analyses_diversity/_child_model_simpson.Rmd
@@ -0,0 +1,66 @@
+### `r stringr::str_to_sentence("{{group}} - {{primerset}} - {{unit}}")`
+
+```{r}
+data_subset <- combined %>%
+  filter(
+    group == "{{group}}",
+    primerset == "{{primerset}}",
+    unit == "{{unit}}"
+  )
+```
+
+
+```{r m-simpson-{{group}}-{{primerset}}-{{unit}}}
+form_simpson <- formula(
+  simpson ~
+    log(total_count)
+  + landgebruik
+  + diepte
+  + landgebruik:diepte
+  + (1 | cmon_plot_id)
+  )
+
+cat("Fitting ordinal beta model")
+
+m_simpson <- glmmTMB(
+  formula = form_simpson,
+  data = data_subset,
+  family = ordbeta())
+```
+
+
+```{r m-simpson-checks-{{group}}-{{primerset}}-{{unit}}, error = TRUE}
+performance::check_model(m_simpson)
+```
+
+
+```{r m-simpson-summary-{{group}}-{{primerset}}-{{unit}}}
+summary(m_simpson) |> print(digits = 2)
+```
+
+```{r m-simpson-anova-{{group}}-{{primerset}}-{{unit}}, error = TRUE}
+car::Anova(m_simpson) |> print(digits = 2)
+```
+
+```{r m-simpson-preds-total-count-{{group}}-{{primerset}}-{{unit}}, error = TRUE}
+marginaleffects::plot_predictions(
+  m_simpson,
+  condition = c("total_count"),
+  re.form = NA,
+  vcov = TRUE,
+  type = "response") +
+  labs(y = "Simpson index") +
+  scale_x_log10()
+```
+
+```{r m-simpson-preds-landgebruikxdiepte-{{group}}-{{primerset}}-{{unit}}, error = TRUE}
+marginaleffects::plot_predictions(
+  m_simpson,
+  condition = c("landgebruik", "diepte"),
+  re.form = NA,
+  vcov = TRUE,
+  type = "response") +
+  labs(y = "Simpson index")
+```
+
+
diff --git a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
new file mode 100644
index 0000000..46d2fb6
--- /dev/null
+++ b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
@@ -0,0 +1,210 @@
+---
+title: "Analyses diversiteit"
+author: "Sam Lambrechts, Sylke De Backer, Io Deflem, Emma Cartuyvels, Hans Van Calster"
+date: "`r Sys.Date()`"
+output:
+  bookdown::html_document2:
+    toc: true
+    toc_float: true
+    code_folding: hide
+editor_options:
+  markdown:
+    wrap: sentence
+---
+
+```{r setup, include=FALSE}
+library(knitr)
+library(here)
+opts_chunk$set(echo = TRUE, out.width = "100%")
+opts_knit$set(root.dir = here::here())
+
+library(ggplot2)
+library(ggforce)
+library(dplyr)
+library(tidyr)
+library(purrr)
+library(glmmTMB)
+library(marginaleffects)
+library(performance)
+
+mbag_bodem_folder <- "G:/Gedeelde drives/PRJ_MBAG/4c_bodembiodiversiteit" # nolint
+```
+
+# Inlezen data
+
+```{r inlezen}
+diversiteit <- readr::read_csv(
+  file.path(
+    mbag_bodem_folder,
+    "data", "statistiek", "dataframe_overkoepelend",
+    "mbag_combined_dataframe.csv")
+)
+
+metadata <- readxl::read_excel(
+  file.path(
+    mbag_bodem_folder,
+    "data", 
+    "Stratificatie_MBAG_plots",
+    "MBAG_stratfile_v1_cleaned.xlsx")
+) %>%
+  janitor::clean_names() %>%
+  rename(landgebruik_ruimtebeslag = landgebruik) %>%
+  mutate(
+    landgebruik = factor(
+      landgebruik_mbag,
+      levels = c(
+        "Akker", "Tijdelijk grasland", "Blijvend grasland",
+        "Residentieel grasland", "Natuurgrasland")
+    ),
+    diepte = gsub("-", "_", diepte) |> factor()
+  )
+
+combined <- diversiteit %>%
+  inner_join(metadata, by = "sample")
+```
+
+```{r}
+glimpse(combined)
+```
+
+# Verkenning
+
+Het geobserveerde aantal taxa neemt toe met het aantal reads in een staal:
+
+```{r}
+combined %>%
+  filter(total_count > 0) %>%
+  ggplot() +
+  geom_point(aes(x = total_count, y = observed)) +
+  scale_x_log10() +
+  labs(y = "Observed number of taxa") +
+  facet_wrap(group ~ primerset ~ unit, scales = "free")
+```
+
+Relaties tussen de drie diversiteitsmaten:
+
+```{r}
+combined %>%
+  ggplot() +
+  geom_autopoint(
+    aes(
+      colour = group, shape = primerset
+    ),
+    alpha = 0.5
+  ) +
+  geom_autodensity(
+    aes(
+      fill = group
+    ),
+    alpha = 0.5
+  ) +
+  facet_matrix(
+    rows = vars(observed, shannon, simpson), layer.diag = 2
+  )
+```
+
+# Analyses
+
+## Geobserveerde soortenrijkdom
+
+We fitten een model waarbij we de log van het totaal aantal reads in een staal als covariaat toevoegen samen met landgebruik en diepte van het staal (de laatste twee in interactie).
+We voegen ook een random intercept toe die aangeeft op welke locatie een staal werd genomen.
+Hiermee geven we aan dat stalen op verschillende diepte, maar dezelfde locatie, gecorreleerd zijn.
+
+
+```{r child-model-observed-richness, echo=FALSE, results="asis"}
+inputs <- combined %>%
+  distinct(group, primerset, unit) %>%
+  arrange(
+    group, primerset, unit
+  )
+
+pmap(
+  inputs,
+  function(
+    group = group,
+    primerset = primerset,
+    unit = unit) {
+    knit_expand(
+      here("source", "rmarkdown", "analyses_diversity",
+           "_child_model_observed_richness.Rmd"),
+      group = group,
+      primerset = primerset,
+      unit = unit
+    )
+  }
+) %>%
+  paste(collapse = "\n") -> rmd
+
+#clipr::write_clip(rmd)
+
+knit_child(text = rmd, quiet = TRUE) %>%
+  cat()
+```
+
+## Shannon index
+
+```{r child-model-shannon, echo=FALSE, results="asis"}
+inputs <- combined %>%
+  distinct(group, primerset, unit) %>%
+  arrange(
+    group, primerset, unit
+  )
+
+pmap(
+  inputs,
+  function(
+    group = group,
+    primerset = primerset,
+    unit = unit) {
+    knit_expand(
+      here("source", "rmarkdown", "analyses_diversity",
+           "_child_model_shannon.Rmd"),
+      group = group,
+      primerset = primerset,
+      unit = unit
+    )
+  }
+) %>%
+  paste(collapse = "\n") -> rmd
+
+#clipr::write_clip(rmd)
+
+knit_child(text = rmd, quiet = TRUE) %>%
+  cat()
+```
+
+## Simpson index
+
+```{r child-model-simpson, echo=FALSE, results="asis"}
+inputs <- combined %>%
+  distinct(group, primerset, unit) %>%
+  arrange(
+    group, primerset, unit
+  )
+
+pmap(
+  inputs,
+  function(
+    group = group,
+    primerset = primerset,
+    unit = unit) {
+    knit_expand(
+      here("source", "rmarkdown", "analyses_diversity",
+           "_child_model_simpson.Rmd"),
+      group = group,
+      primerset = primerset,
+      unit = unit
+    )
+  }
+) %>%
+  paste(collapse = "\n") -> rmd
+
+#clipr::write_clip(rmd)
+
+knit_child(text = rmd, quiet = TRUE) %>%
+  cat()
+```
+
+
+

From 685be38d99ee106f2fb5409cca2dcfc5da787be5 Mon Sep 17 00:00:00 2001
From: hansvancalster <hans.vancalster@inbo.be>
Date: Fri, 23 Aug 2024 13:34:37 +0200
Subject: [PATCH 2/3] new metadata

---
 .../rmarkdown/analyses_diversity/analyses_diversity.Rmd  | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
index 46d2fb6..9b23af6 100644
--- a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
+++ b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
@@ -15,7 +15,7 @@ editor_options:
 ```{r setup, include=FALSE}
 library(knitr)
 library(here)
-opts_chunk$set(echo = TRUE, out.width = "100%")
+opts_chunk$set(echo = TRUE, error = TRUE, out.width = "100%")
 opts_knit$set(root.dir = here::here())
 
 library(ggplot2)
@@ -40,21 +40,20 @@ diversiteit <- readr::read_csv(
     "mbag_combined_dataframe.csv")
 )
 
-metadata <- readxl::read_excel(
+metadata <- readr::read_csv(
   file.path(
     mbag_bodem_folder,
     "data", 
     "Stratificatie_MBAG_plots",
-    "MBAG_stratfile_v1_cleaned.xlsx")
+    "MBAG_stratfile_v2_cleaned_12.csv")
 ) %>%
   janitor::clean_names() %>%
-  rename(landgebruik_ruimtebeslag = landgebruik) %>%
   mutate(
     landgebruik = factor(
       landgebruik_mbag,
       levels = c(
         "Akker", "Tijdelijk grasland", "Blijvend grasland",
-        "Residentieel grasland", "Natuurgrasland")
+        "Residentieel grasland", "Natuurgrasland", "Heide")
     ),
     diepte = gsub("-", "_", diepte) |> factor()
   )

From bf5a9326cd902b1b559bd7624cf1bd739ee2a9fb Mon Sep 17 00:00:00 2001
From: hansvancalster <hans.vancalster@inbo.be>
Date: Fri, 23 Aug 2024 13:41:59 +0200
Subject: [PATCH 3/3] fix gsub

---
 source/rmarkdown/analyses_diversity/analyses_diversity.Rmd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
index 9b23af6..9fe9c2f 100644
--- a/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
+++ b/source/rmarkdown/analyses_diversity/analyses_diversity.Rmd
@@ -55,7 +55,7 @@ metadata <- readr::read_csv(
         "Akker", "Tijdelijk grasland", "Blijvend grasland",
         "Residentieel grasland", "Natuurgrasland", "Heide")
     ),
-    diepte = gsub("-", "_", diepte) |> factor()
+    diepte = gsub("_|/", "-", diepte) |> factor()
   )
 
 combined <- diversiteit %>%