Fix spring analysis

source/_verkenning.Rmd

@@ -125,6 +125,7 @@ venn({{group}}, "month")
 ### Soortenaccumulatie
 
 #### Saturatie pan traps in tijd
+
 <!-- spell-check: ignore:start -->
 >Moreover, some rarefaction and extrapolation curves may cross one or more times, so that the rank order of diversity measured among samples could change depending on the sampling effort that is used for standardized comparisons (Chao and Jost 2012). 
 <!-- spell-check: ignore:end -->
@@ -409,17 +410,27 @@ rarefaction_result_3 %>%
 cat("### Functionele groepen\n")
 
 nat_ap <- {{group}} %>%
-  left_join(naturalhistorytraits, join_by(species_nm == species))
+  left_join(naturalhistorytraits, join_by(species_nm == species)) %>%
+  mutate(
+    traits_combined = paste(
+      foraging_prefence,
+      sociality,
+      nesting_type,
+      tongue_length
+    )
+  )
 ```
 
 ```{r, fig.cap = "Unieke functionele groepen apoidea voor transecten (TS) vs. pan traps (PT).", eval="{{group}}"=="apoidea"}
-venn(nat_ap, "method_cd", "optional_functional_group")
+
+venn(nat_ap, "method_cd", "traits_combined")
 ```
 
 ```{r, fig.cap = "Unieke functionele groepen apoidea per submethode.", eval="{{group}}"=="apoidea"}
-venn(nat_ap, "method_combi", "optional_functional_group")
+venn(nat_ap, "method_combi", "traits_combined")
 ```
 
 ```{r, fig.cap = "Unieke functionele groepen apoidea per maand.", eval="{{group}}"=="apoidea"}
-venn(nat_ap, "month", "optional_functional_group")
+venn(nat_ap, "month", "traits_combined")
 ```
+

source/data_analysis_spring_2023.Rmd

@@ -16,8 +16,6 @@
 ```{r setup, include=FALSE}
 library(knitr)
 opts_chunk$set(message = FALSE, warning = FALSE)
-opts_knit$set(
-  root.dir = here::here("source"))
 options(scipen = 10)
 
 library(RODBC)
@@ -68,7 +66,8 @@
   ggVennDiagram::ggVennDiagram(x, label_alpha = 0, edge_size = 0.75) +
     scale_fill_distiller(palette = "Greens", direction = 1) +
     scale_x_continuous(expand = c(0.2, 0.2)) +
-    theme(legend.position = "none")
+    theme(legend.position = "none") +
+    labs(title = paste0("Aantal unieke van ", count_id))
 }
 ```
 
@@ -106,7 +105,7 @@
 Van de eerste twee tabellen wordt een versimpelde versie gemaakt.
 
 ```{r load-data}
-conn <- odbcConnectAccess2007(here("data/SPRING.accdb"))
+conn <- odbcConnectAccess2007("../data/SPRING.accdb")
 
 identifications <- sqlFetch(conn, "identifications")
 samples <- sqlFetch(conn, "samples")
@@ -480,7 +479,7 @@
   function(
       title = title, group = group, familyvec = familyvec) {
     knit_expand(
-      here::here("source", "_verkenning.Rmd"),
+      "_verkenning.Rmd",
       title = title,
       group = group,
       familyvec = familyvec
@@ -494,8 +493,21 @@
 # enkel nodig indien je interactief werkt en de code van deze chunks nodig hebt
 # clipr::write_clip(rmd) # nolint
 
-knit(text = rmd, quiet = TRUE) %>%
-  cat()
+
+# Function to execute R code chunks from the generated rmd
+execute_rmd_chunks <- function(rmd_text) {
+  # Extract R code from the rmd content
+  r_code <- knitr::purl(text = rmd, quiet = TRUE)
+  eval(parse(text = r_code), envir = .GlobalEnv)
+}
+
+# Execute the generated rmd chunks
+if (interactive()) {
+  execute_rmd_chunks(rmd)
+} else {
+  knit(text = rmd, quiet = TRUE) %>%
+    cat()
+}
 ```
 
 
@@ -749,7 +761,7 @@
 ## Effect aantal bloemen
 
 ```{r}
-flowers <- read_xlsx(here::here("data", "Number of flowers.xlsx")) |>
+flowers <- read_xlsx("../data/Number of flowers.xlsx") |>
    janitor::clean_names()
 ```
 
@@ -767,45 +779,107 @@
             by = join_by(sampling_site_cd == sample_code))
 ```
 
+
+```{r}
+flowers <- flowers |>
+  mutate(location_code = str_sub(sampling_site_cd, start = 1L, end = 8L),
+         maand = str_sub(sampling_site_cd, start = 20L, end = 20L))
+```
+
+Aantal soorten bijen ifv aantal bloemen:
+
 ```{r}
 flowers |>
   ggplot(aes(x = number_of_floral_units, y = n_species_ap)) +
   geom_point() +
-  geom_smooth(method = "lm")
+  geom_smooth(
+    method = "glm", method.args = list(family = "poisson"), se = FALSE) +
+  scale_x_continuous(trans = "pseudo_log") +
+  facet_grid(location_code ~ maand, scales = "free")
 ```
 
+Aantal individuen van bijen ifv aantal bloemen:
+
 ```{r}
-model1 <- lm(n_ind_ap ~ number_of_floral_units, data = flowers |>
-               filter(number_of_floral_units < 4000))
-summary(model1)
+flowers |>
+  ggplot(aes(x = number_of_floral_units, y = n_ind_ap)) +
+  geom_point() +
+  geom_smooth(
+    method = "glm", method.args = list(family = "poisson"), se = FALSE) +
+  scale_x_continuous(trans = "pseudo_log") +
+  facet_grid(location_code ~ maand, scales = "free")
+```
 
-model2 <- lm(n_species_ap ~ number_of_floral_units, data = flowers |>
-               filter(number_of_floral_units < 4000))
-summary(model2)
+Aantal soorten zweefvliegen ifv aantal bloemen:
 
-model3 <- lm(n_ind_syr ~ number_of_floral_units, data = flowers |>
-               filter(number_of_floral_units < 4000))
-summary(model3)
+```{r}
+flowers |>
+  ggplot(aes(x = number_of_floral_units, y = n_species_syr)) +
+  geom_point() +
+  geom_smooth(
+    method = "glm", method.args = list(family = "poisson"), se = FALSE) +
+  scale_x_continuous(trans = "pseudo_log") +
+  facet_grid(location_code ~ maand, scales = "free")
+```
 
-model4 <- lm(n_species_syr ~ number_of_floral_units, data = flowers |>
-               filter(number_of_floral_units < 4000))
-summary(model4)
+Aantal individuen van zweefvliegen ifv aantal bloemen:
+
+```{r}
+flowers |>
+  ggplot(aes(x = number_of_floral_units, y = n_ind_syr)) +
+  geom_point() +
+  geom_smooth(
+    method = "glm", method.args = list(family = "poisson"), se = FALSE) +
+  scale_x_continuous(trans = "pseudo_log") +
+  facet_grid(location_code ~ maand, scales = "free")
 ```
 
-Geen significant effect van het aantal bloemeenheden op het aantal soorten of het aantal gevangen individuen; het effect wordt nog minder significant wanneer één sterke uitschieter wordt verwijderd.
+```{r}
+flowers |>
+  ggplot(
+        aes(
+      x = maand,
+      y = number_of_floral_units
+    )
+  ) +
+  geom_violin() +
+  ggforce::geom_sina() +
+  scale_y_continuous(
+    trans = "pseudo_log",
+    breaks = c(0, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000))
+```
 
-### Uitgebreid model
 ```{r}
-flowers <- flowers |>
-  mutate(location_code = str_sub(sampling_site_cd, start = 1L, end = 8L),
-         maand = str_sub(sampling_site_cd, start = 20L, end = 20L))
+flowers |>
+  ggplot(
+    aes(
+      x = n_ind_ap,
+      y = n_species_ap
+    )
+  ) +
+  stat_sum() +
+  geom_abline()
+
+flowers |>
+  ggplot(
+    aes(
+      x = n_ind_syr,
+      y = n_species_syr
+    )
+  ) +
+  stat_sum() +
+  geom_abline()
 ```
 
 
 ```{r}
 model1 <- glmmTMB(
-  n_ind_ap ~ log(number_of_floral_units + 1) + location_code + maand,
-  ziformula = ~ 1,
+  n_species_ap ~
+    log(number_of_floral_units + 1) 
+  + maand
+  + location_code
+  ,
+  ziformula = ~ maand,
   family = "poisson",
   na.action = na.exclude,
   data = flowers
@@ -817,13 +891,27 @@
 ```
 
 ```{r}
+performance::check_overdispersion(model1)
 performance::check_model(model1)
 ```
 
+```{r}
+marginaleffects::plot_predictions(
+  model1, "number_of_floral_units",
+  vcov = TRUE,
+  re.form = NA
+) +
+  scale_x_continuous(trans = "pseudo_log")
+```
+
+
 ```{r}
 model2 <- glmmTMB(
-  n_species_syr ~ number_of_floral_units + location_code + maand,
-  ziformula = ~ 1,
+  n_species_syr ~
+    log(number_of_floral_units + 1)
+  + location_code
+  + maand,
+  ziformula = ~ maand,
   family = "poisson",
   na.action = na.exclude,
   data = flowers
@@ -835,9 +923,81 @@
 ```
 
 ```{r}
+performance::check_overdispersion(model2)
 performance::check_model(model2)
 ```
 
+```{r}
+marginaleffects::plot_predictions(
+  model2, "number_of_floral_units",
+  vcov = TRUE,
+  re.form = NA
+) +
+  scale_x_continuous(trans = "pseudo_log")
+```
+
+```{r}
+model3 <- glmmTMB(
+  n_ind_ap ~
+    log(number_of_floral_units + 1)
+  + location_code
+  + maand,
+  ziformula = ~ maand,
+  family = "nbinom2",
+  na.action = na.exclude,
+  data = flowers
+)
+```
+
+```{r}
+summary(model3)
+```
+
+```{r}
+performance::check_overdispersion(model3)
+performance::check_model(model3)
+```
+
+```{r}
+marginaleffects::plot_predictions(
+  model3, "number_of_floral_units",
+  vcov = TRUE,
+  re.form = NA
+) +
+  scale_x_continuous(trans = "pseudo_log")
+```
+
+```{r}
+model4 <- glmmTMB(
+  n_ind_syr ~
+  + log(number_of_floral_units + 1)
+  + location_code
+  + maand,
+  ziformula = ~ maand,
+  family = "nbinom2",
+  na.action = na.exclude,
+  data = flowers
+)
+```
+
+```{r}
+summary(model4)
+```
+
+```{r}
+performance::check_overdispersion(model4)
+performance::check_model(model4)
+```
+
+```{r}
+marginaleffects::plot_predictions(
+  model4, "number_of_floral_units",
+  vcov = TRUE,
+  re.form = NA
+) +
+  scale_x_continuous(trans = "pseudo_log")
+```
+
 # Kosten
 
 ```{r}