From f1baee1e66fead4ab2ef9b6b5bad067af9893f4a Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Wed, 13 Nov 2024 16:35:32 +0100
Subject: [PATCH 01/13] initial commit

---
 .../markdown/test_densiteitsmodellering.Rmd   | 62 +++++++++++++++++++
 1 file changed, 62 insertions(+)
 create mode 100644 source/markdown/test_densiteitsmodellering.Rmd

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
new file mode 100644
index 00000000..5fc4b07e
--- /dev/null
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -0,0 +1,62 @@
+---
+title: "Test density estimation and modelling"
+author: "Ward Langeraert"
+date: "`r Sys.Date()`"
+output:
+  bookdown::html_document2:
+    code_folding: hide
+    toc: true
+    toc_float: true
+    toc_collapsed: true
+knit: (function(inputFile, encoding) {
+  rmarkdown::render(inputFile, encoding = encoding, output_dir = "../../output/rapporten/markdown/2024") })
+editor_options: 
+  chunk_output_type: console
+---
+
+```{r setup, include=FALSE}
+# set up
+library(knitr)
+opts_chunk$set(
+  echo = TRUE,
+  warning = FALSE,
+  message = FALSE,
+  out.width = "100%"
+)
+opts_knit$set(root.dir = here::here())
+```
+
+```{r}
+# packages
+library(tidyverse)
+library(sf)
+library(targets)
+
+# Conflicts
+conflicted::conflicts_prefer(dplyr::filter)
+conflicted::conflicts_prefer(dplyr::select)
+
+# Source
+source(here::here("source", "R", "tar_read_sf.R"))
+
+# Paths
+mbag_dir <- here::here()
+targets_store <- here::here("source", "targets", "data_preparation", "_targets")
+```
+
+# Achtergrond
+
+In kader van het meetnet agrarische soorten (MAS), worden vogels en Haas gemonitord volgens een vast protocol. Hierbij geven tellers waarnemingen in vanaf een telpunt tot 300 meter ver. Dit doen ze 4x per jaar. Doordat we de afstand van de teller tot de waarneming kennen, kunnen we de detectiekans per soort schatten en vervolgens de werkelijke densiteit berekenen (incl. onzekerheid).
+
+Het meetnet dekt verschillende strata: landbouwregio’s, soortbeschermingsplan (SBP; binnen of buiten) en openheid landschap (OL: open landschap, HOL: half-open landschap).
+
+We willen de volgende onderzoeksvragen beantwoorden.
+
+- Wat is de detectiekans van de soorten binnen het MAS?
+- Wat is de densiteit van soorten binnen het MAS?
+  - Per stratum
+  - Per landbouwregio
+  - Vlaanderen (steekproefkader)
+- Hoe is de densiteit van doelsoorten verspreid in Vlaanderen?
+- Kunnen we kerngebieden afbakenen voor doelsoorten binnen Vlaanderen?
+

From b4d6898d46cb3980055eabd3ce7ccc792f163516 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Wed, 13 Nov 2024 17:11:00 +0100
Subject: [PATCH 02/13] add breeding period file

---
 data/SOVON/README.txt | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/data/SOVON/README.txt b/data/SOVON/README.txt
index 8d3d92ea..202645a6 100644
--- a/data/SOVON/README.txt
+++ b/data/SOVON/README.txt
@@ -3,6 +3,9 @@ This folder contains data files related Sovon (https://www.sovon.nl/).
 - avimap_601_0_MAS_Vlaanderen_telpunten_xy.shp
   - Existing MAS locations before monitoring network
   - Used for recycling in sampling design
+- Interpretatie_Criteria_Broedvogels_v2.csv
+  - Breeding dates birds
+  - Used for data selection distance sampling
 
 - Any files related to Sovon for input or output in scripts 
 

From 5c4a15d1f757b5ee2468b2282543ae0a8744a2c2 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Wed, 13 Nov 2024 17:11:11 +0100
Subject: [PATCH 03/13] add data exploration

---
 .../markdown/test_densiteitsmodellering.Rmd   | 149 +++++++++++++++++-
 1 file changed, 146 insertions(+), 3 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index 5fc4b07e..cfc31a70 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -32,13 +32,13 @@ library(tidyverse)
 library(sf)
 library(targets)
 
+library(INBOtheme)
+theme_set(theme_inbo(transparent = TRUE))
+
 # Conflicts
 conflicted::conflicts_prefer(dplyr::filter)
 conflicted::conflicts_prefer(dplyr::select)
 
-# Source
-source(here::here("source", "R", "tar_read_sf.R"))
-
 # Paths
 mbag_dir <- here::here()
 targets_store <- here::here("source", "targets", "data_preparation", "_targets")
@@ -60,3 +60,146 @@ We willen de volgende onderzoeksvragen beantwoorden.
 - Hoe is de densiteit van doelsoorten verspreid in Vlaanderen?
 - Kunnen we kerngebieden afbakenen voor doelsoorten binnen Vlaanderen?
 
+We willen deze analyses automatiseren via een [targets pipeline](https://books.ropensci.org/targets/) die we branchen per jaar en per soort.
+Hiervoor moeten we eerst enkele dingen uittesten:
+
+- Kunnen we de detectiecurve fitten op alle data en de densiteitsschatting op een subset van de data toepassen (maxima)?
+- Hoe creëren we een spatiaal model voor densiteiten? Hoe nemen we de detectiekans en onzekerheid mee?
+
+# Test data
+
+We doen de analyses voor 1 branch van de targets pipeline. Namelijk voor de Veldleeuwerik in 2024.
+We selecteren alle waarnemingen met broedcode > 0.
+Alle aantallen met broedcode > 0 zijn 1.
+
+> Broedcodes ok?
+
+```{r}
+# Load mas data
+mas_data_clean <- tar_read("mas_data_clean", store = targets_store)
+
+# Select data Veldleeuwerik 2024
+veldleeuwerik_2024_df <- mas_data_clean %>%
+  filter(naam == "Veldleeuwerik",
+         jaar == 2024,
+         wrntype > 0) %>%
+  mutate(aantal = 1,
+         regio = ifelse(grepl("\\sleemstreek$", regio), "Leemstreek", regio),
+         stratum = paste(openheid_klasse, sbp, sep = " - "))
+```
+
+We kijken naar het aantal broedparen per regio.
+
+```{r}
+# Visualise
+veldleeuwerik_2024_df %>%
+  ggplot(aes(x = periode_in_jaar, fill = wrntype)) +
+    geom_bar() +
+    facet_wrap(~regio, scales = "free") +
+    labs(x = "Telperiode", y = "Aantal broedparen", fill = "broedcode")
+```
+
+En de per stratum.
+
+```{r}
+# Visualise
+veldleeuwerik_2024_df %>%
+  ggplot(aes(x = periode_in_jaar, fill = stratum)) +
+    geom_bar() +
+    facet_wrap(~regio, scales = "free") +
+    labs(x = "Telperiode", y = "Aantal broedparen", fill = "Stratum")
+```
+
+Alle telperiodes bevinden zich binnen de datumgrenzen.
+We nemen alle telperiode dus mee voor de densiteitsschattingen.
+
+```{r}
+# Load breeding dates file
+datumgrenzen_df <- read_csv2(
+  file.path(mbag_dir, "data", "SOVON",
+            "Interpretatie_Criteria_Broedvogels_v2.csv"))
+
+# Tidy dataframe
+trans_vec <- c("jan", "feb", "mrt", "apr", "mei", "jun", "jul", "aug")
+
+datumgrenzen_df2 <- datumgrenzen_df %>%
+  mutate(across(Datum_begin:`Datum eind_oud`, ~ gsub("-", "/", .x))) %>%
+  separate(Datum_begin, into = c("Datum_begin_dag", "Datum_begin_maand"), 
+           sep = "/") %>%
+  separate(Datum_eind, into = c("Datum_eind_dag", "Datum_eind_maand"), 
+           sep = "/") %>%
+  separate(`Datum begin_oud`, into = c("Datum_begin_dag_oud", 
+                                       "Datum_begin_maand_oud"), 
+           sep = "/") %>%
+  separate(`Datum eind_oud`, into = c("Datum_eind_dag_oud", 
+                                      "Datum_eind_maand_oud"), 
+           sep = "/") %>%
+  mutate(Datum_begin_maand = match(Datum_begin_maand, trans_vec),
+         datum_begin = paste(Datum_begin_maand, Datum_begin_dag, sep = "-")) %>%
+  mutate(Datum_eind_maand = match(Datum_eind_maand, trans_vec),
+         datum_eind = paste(Datum_eind_maand, Datum_eind_dag, sep = "-")) %>%
+  mutate(Datum_begin_maand_oud = match(Datum_begin_maand_oud, trans_vec),
+         datum_begin_oud = paste(Datum_begin_maand_oud, Datum_begin_dag_oud, 
+                                 sep = "-")) %>%
+  mutate(Datum_eind_maand_oud = match(Datum_eind_maand_oud, trans_vec),
+         datum_eind_oud = paste(Datum_eind_maand_oud, Datum_eind_dag_oud,
+                                sep = "-")) %>%
+  mutate(across(datum_begin:datum_eind_oud, 
+                ~ ifelse(.x == "NA-NA", NA, .x))) %>%
+  select(id = Id, soort = Soort, datum_begin, datum_eind, datum_begin_oud, 
+         datum_eind_oud, broedcode, broedcode_oud = bc_oud, 
+         wnm_vereist = `Waarnemingen verreist`, Opmerking)
+
+# MAS dates
+r1_start <- "04-01"
+r1_stop <- "04-20"
+r2_start <- "04-21"
+r2_stop <- "05-10"
+r3_start <- "05-11"
+r3_stop <- "06-10"
+r4_start <- "06-21"
+r4_stop <- "07-15"
+
+# Classify
+datumgrenzen_mas <- datumgrenzen_df2 %>%
+  select(id, soort, datum_begin, datum_eind) %>%
+  mutate(
+    jaar = 2024,
+    datum_begin2 = ymd(paste(jaar, datum_begin, sep = "-")),
+    datum_eind2 = ymd(paste(jaar, datum_eind, sep = "-")),
+    R1 = ifelse(ymd(paste(jaar, r1_start, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2) &
+                ymd(paste(jaar, r1_stop, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2), TRUE, FALSE),
+    R2 = ifelse(ymd(paste(jaar, r2_start, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2) &
+                ymd(paste(jaar, r2_stop, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2), TRUE, FALSE),
+    R3 = ifelse(ymd(paste(jaar, r3_start, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2) &
+                ymd(paste(jaar, r3_stop, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2), TRUE, FALSE),
+    R4 = ifelse(ymd(paste(jaar, r4_start, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2) &
+                ymd(paste(jaar, r4_stop, sep = "-")) %within% 
+                  interval(datum_begin2, datum_eind2), TRUE, FALSE)
+    
+    ) %>%
+  select(id, soort, datum_begin, datum_eind, R1, R2, R3, R4)
+
+# Show results Veldleeuwerik
+datumgrenzen_mas %>%
+  filter(soort == "Veldleeuwerik") %>%
+  kable()
+```
+
+Ten slotte kijken we naar de afstanden.
+
+```{r}
+# Visualise
+veldleeuwerik_2024_df %>%
+  ggplot(aes(x = distance2plot, fill = stratum)) +
+    geom_histogram() +
+    facet_wrap(~regio, scales = "free") +
+    labs(x = "Afstand (m)", y = "Aantal broedparen", fill = "Stratum")
+```

From ad96f1e6bfae3a49c54a51ee8fbe4edba1be5d03 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Wed, 13 Nov 2024 17:46:49 +0100
Subject: [PATCH 04/13] data preparation

---
 .../markdown/test_densiteitsmodellering.Rmd   | 53 ++++++++++++++++++-
 1 file changed, 51 insertions(+), 2 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index cfc31a70..256bf671 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -29,8 +29,9 @@ opts_knit$set(root.dir = here::here())
 ```{r}
 # packages
 library(tidyverse)
-library(sf)
+library(Distance)
 library(targets)
+library(sf)
 
 library(INBOtheme)
 theme_set(theme_inbo(transparent = TRUE))
@@ -85,7 +86,8 @@ veldleeuwerik_2024_df <- mas_data_clean %>%
          wrntype > 0) %>%
   mutate(aantal = 1,
          regio = ifelse(grepl("\\sleemstreek$", regio), "Leemstreek", regio),
-         stratum = paste(openheid_klasse, sbp, sep = " - "))
+         stratum = paste(openheid_klasse, sbp, sep = " - ")) %>%
+  st_drop_geometry()
 ```
 
 We kijken naar het aantal broedparen per regio.
@@ -203,3 +205,50 @@ veldleeuwerik_2024_df %>%
     facet_wrap(~regio, scales = "free") +
     labs(x = "Afstand (m)", y = "Aantal broedparen", fill = "Stratum")
 ```
+
+# Distance sampling
+
+We schatten de detectiekans en densiteiten via distance sampling.
+We vergelijken een tweestapsmethode met **mrds** met de directe methode van **Distance**.
+We beschouwen voor de formulatie van de detectiefunctie de half-normal en hazard-rate sleutelfuncties die we laten afhangen van `openheid`.
+
+## Data preparatie
+
+We hebben een kolom `object` nodig: volgnummer voor elke waarneming; een kolom `size`: aantal broedparen per waarneming (overal 1); een kolom `distance`: afstand tot elke waarneming; kolommen met covariaten (hier 1 variabele): `openheid`.
+
+```{r}
+veldleeuwerik_distance <- veldleeuwerik_2024_df %>%
+  select(object = oid,
+         size = aantal,
+         distance = distance2plot,
+         openheid = openheid_klasse)
+```
+
+Om abundanties te schatten, moeten we de oppervlakte van de strata berekenen voor extrapolatie. We trekken een buffer van 300 m rond het steekproefkader.
+
+## Model specificatie
+### MRDS
+
+```{r}
+test <- ddf(
+  data = veldleeuwerik_distance,
+  dsmodel = ~mcds(key="hr", formula = ~openheid),
+  method = "ds",
+  meta.data = list(width = 300, point = TRUE))
+summary(test)
+
+ddf.gof(test)
+plot(test, pdf = TRUE)
+```
+
+### Distance
+
+## Model selectie
+### MRDS
+
+```{r}
+ddf.gof(test)
+plot(test, pdf = TRUE)
+```
+
+### Distance

From 27caf320a3fe4ca14c9bce5d0589389244e82b45 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Thu, 14 Nov 2024 16:12:04 +0100
Subject: [PATCH 05/13] create dataframes for density estimation

---
 .../markdown/test_densiteitsmodellering.Rmd   | 143 ++++++++++++++++++
 1 file changed, 143 insertions(+)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index 256bf671..afe6e55d 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -226,6 +226,149 @@ veldleeuwerik_distance <- veldleeuwerik_2024_df %>%
 
 Om abundanties te schatten, moeten we de oppervlakte van de strata berekenen voor extrapolatie. We trekken een buffer van 300 m rond het steekproefkader.
 
+```{r}
+# Lees steekproefkader in
+steekproefkader <- st_read(
+  file.path(mbag_dir, "data", "steekproefkaders",
+            "steekproefkader_mbag_mas.gpkg"))
+
+# Buffer and area
+## Calculate area per region
+region_sf <- steekproefkader %>%
+  mutate(regio = ifelse(grepl("\\sleemstreek$", regio),
+                        "Leemstreek", regio)) %>%
+  st_buffer(dist = 300) %>%
+  group_by(regio) %>%
+  summarise(geom = st_union(geom)) %>%
+  ungroup() %>%
+  mutate(Area = as.numeric(st_area(geom)) / 1e6) %>%
+  select(regio, Area, everything())
+
+## Calculate area for Flanders
+flanders_sf <- steekproefkader %>%
+  st_buffer(dist = 300) %>%
+  summarise(geom = st_union(geom)) %>%
+  ungroup() %>%
+  mutate(regio = "Flanders",
+         Area = as.numeric(st_area(geom)) / 1e6) %>%
+  select(regio, Area, everything())
+
+## Calculate area per stratum
+strata_sf <- steekproefkader %>%
+  mutate(regio = ifelse(grepl("\\sleemstreek$", regio),
+                        "Leemstreek", regio)) %>%
+  st_buffer(dist = 300) %>%
+  group_by(regio, "openheid" = openheid_klasse, sbp) %>%
+  summarise(geom = st_union(geom)) %>%
+  ungroup() %>%
+  mutate(Area = as.numeric(st_area(geom)) / 1e6) %>%
+  select(regio, openheid, sbp, Area, everything())
+```
+
+Voor Vlaanderen:
+
+```{r}
+flanders_sf %>%
+  st_drop_geometry() %>%
+  rename("opp. (km²)" = "Area") %>%
+  kable(digits = 3)
+```
+
+```{r}
+# Visualisation
+flanders_sf %>%
+  ggplot() +
+    geom_sf(aes(fill = regio)) +
+    labs(fill = "") +
+    theme(legend.position = "bottom")
+```
+
+Per landbouwregio:
+
+```{r}
+region_sf %>%
+  st_drop_geometry() %>%
+  mutate(`tot. opp. (km²)` = sum(Area)) %>%
+  rename("opp. (km²)" = "Area") %>%
+  kable(digits = 3)
+```
+
+Als we de oppervlakte optellen komen we dezelfde oppervlakte uit als volledig Vlaanderen.
+
+```{r}
+# Visualisation
+region_sf %>%
+  ggplot() +
+    geom_sf(aes(fill = regio)) +
+    labs(fill = "") +
+    theme(legend.position = "bottom")
+```
+
+Per stratum:
+
+```{r}
+strata_sf %>%
+  st_drop_geometry() %>%
+  group_by(regio) %>%
+  mutate(`tot. opp. regio (km²)` = sum(Area)) %>%
+  ungroup() %>%
+  mutate(`totale oppervlakte (km²)` = sum(Area)) %>%
+  rename("opp. (km²)" = "Area") %>%
+  kable(digits = 3)
+```
+
+Als we de oppervlakte optellen komen we niet dezelfde oppervlakte uit als voor de regio's en volledig Vlaanderen!
+
+```{r, fig.width = 10}
+# Visualisation
+for (region in sort(unique(strata_sf$regio))) {
+  p <- strata_sf %>%
+    filter(regio == region) %>%
+    mutate(stratum = paste(openheid, sbp, sep = " - ")) %>%
+    ggplot() +
+      geom_sf(aes(fill = stratum)) +
+      facet_wrap(~regio, ncol = 2) +
+      labs(fill = "") +
+      theme(legend.position = "bottom")
+  print(p)
+}
+```
+
+We willen een densiteitsschatting voor elk stratum. We maken daarom dataframes om de abundanties en densiteiten te berekenen. `Region_table` met de oppervlaktes per stratum (regio-openheid-sbp), `sample_table` geeft aan welke plots in welke strata voorkomen en welke effort is geleverd (effort is gelijk aan het aantal beschouwde telperiodes per jaar) en `obs_table` die aangeeft welke waarnemingen in welke plots en strata zitten. We selecteren enkel de data per telpunt in de periode met het meeste waarnemingen.
+
+```{r}
+design <- read_csv(
+  file.path(mbag_dir, "data", "steekproefkaders",
+            "steekproef_avimap_mbag_mas.csv"))
+
+# Oppervlakte per stratum
+region_table <- strata_sf %>%
+  st_drop_geometry() %>%
+  mutate(stratum = paste(openheid, sbp, sep = " - ")) %>%
+  mutate(Region.Label = paste(regio, openheid, sbp, sep = " - ")) %>%
+  select(Region.Label, Area)
+
+# Effort per sample en plot per stratum
+effort <- 1
+
+sample_table <- design %>%
+  distinct(pointid, regio, openheid = openheid_klasse, sbp) %>%
+  mutate(Region.Label = paste(regio, openheid, sbp, sep = " - "),
+         Effort = effort) %>%
+  select(Sample.Label = pointid, Region.Label, Effort)
+
+# Waarnemingen per plots en strata
+obs_table <- veldleeuwerik_2024_df %>%
+  group_by(periode_in_jaar,  plotnaam) %>%
+  mutate(n = n()) %>%
+  slice_max(order_by = n, n = 1) %>%
+  ungroup() %>%
+  mutate(Region.Label = paste(regio, openheid = openheid_klasse, sbp,
+                              sep = " - ")) %>%
+  select(object = oid, Region.Label, Sample.Label = plotnaam)
+```
+
+
 ## Model specificatie
 ### MRDS
 

From 864d9cb6d40d49ad6187bedb0770ae01c46c0c00 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Thu, 14 Nov 2024 16:44:43 +0100
Subject: [PATCH 06/13] fit distance model

---
 .../markdown/test_densiteitsmodellering.Rmd   | 83 +++++++++++++++----
 1 file changed, 65 insertions(+), 18 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index afe6e55d..85375c54 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -209,7 +209,6 @@ veldleeuwerik_2024_df %>%
 # Distance sampling
 
 We schatten de detectiekans en densiteiten via distance sampling.
-We vergelijken een tweestapsmethode met **mrds** met de directe methode van **Distance**.
 We beschouwen voor de formulatie van de detectiefunctie de half-normal en hazard-rate sleutelfuncties die we laten afhangen van `openheid`.
 
 ## Data preparatie
@@ -318,6 +317,8 @@ strata_sf %>%
 ```
 
 Als we de oppervlakte optellen komen we niet dezelfde oppervlakte uit als voor de regio's en volledig Vlaanderen!
+Dit wil zeggen dat we een apart model moeten fitten voor regio's en specifieke strata.
+Voor Vlaanderen krijgen we een schatting via het regio's model.
 
 ```{r, fig.width = 10}
 # Visualisation
@@ -339,7 +340,8 @@ We willen een densiteitsschatting voor elk stratum. We maken daarom dataframes o
 ```{r}
 design <- read_csv(
   file.path(mbag_dir, "data", "steekproefkaders",
-            "steekproef_avimap_mbag_mas.csv"))
+            "steekproef_avimap_mbag_mas.csv"))  %>%
+  mutate(regio = ifelse(grepl("\\sleemstreek$", regio), "Leemstreek", regio))
 
 # Oppervlakte per stratum
 region_table <- strata_sf %>%
@@ -368,30 +370,75 @@ obs_table <- veldleeuwerik_2024_df %>%
   select(object = oid, Region.Label, Sample.Label = plotnaam)
 ```
 
+Ten slotte zet de conversiefactor de afstanden om van meter naar vierkante kilometer (= 100 ha).
+
+```{r}
+# Conversion factor from meters to 100 hectares
+conversion_factor <- convert_units("meter", NULL, "Square kilometer")
+```
 
 ## Model specificatie
-### MRDS
+
+We beschouwen voor de formulatie van de detectiefunctie de half-normal en hazard-rate sleutelfuncties die we laten afhangen van `openheid`.
 
 ```{r}
-test <- ddf(
-  data = veldleeuwerik_distance,
-  dsmodel = ~mcds(key="hr", formula = ~openheid),
-  method = "ds",
-  meta.data = list(width = 300, point = TRUE))
-summary(test)
-
-ddf.gof(test)
-plot(test, pdf = TRUE)
+veldleeuwerik_dist_hn <- ds(data = veldleeuwerik_distance, key = "hn", 
+  formula = ~openheid, adjustment = NULL, truncation = 300, 
+  transect = "point", dht_group = FALSE, 
+  convert_units = conversion_factor, region_table = region_table, 
+  sample_table = sample_table, obs_table = obs_table)
+
+veldleeuwerik_dist_hr <- ds(data = veldleeuwerik_distance, key = "hr", 
+  formula = ~openheid, adjustment = NULL, truncation = 300, 
+  transect = "point", dht_group = FALSE, 
+  convert_units = conversion_factor, region_table = region_table, 
+  sample_table = sample_table, obs_table = obs_table)
 ```
 
-### Distance
-
 ## Model selectie
-### MRDS
+
+We vergelijken de AIC van de verschillende modellen.
+Als het verschil tussen AIC’s kleiner is dan 2, kiezen we het eenvoudigste van deze modellen (model met minste parameters).
+Het Hazard-rate model scoort beter.
+
+```{r}
+summarize_ds_models(
+    veldleeuwerik_dist_hn,
+    veldleeuwerik_dist_hr
+  ) %>%
+  kable()
+```
+
+## Model fit
+
+Om te controleren of ons model goed past bij de data (“goodness of fit”) maken we een Q-Q plot waarbij de cumulatieve distributiefunctie van de gefitte detectiefunctie (CDF) wordt vergeleken met de distributie van de data (EDF). De Cramér-von Mises test kwantificeert de informatie van de Q-Q-plot door te testen of punten van de EDF en CDF uit dezelfde verdeling komen. Voor studies met punttransecten geven kansdichtheidsfunctieplots ook een goed idee over model fit.
+
+De fit is niet perfect, maar goed genoeg voor deze oefening.
+
+```{r, fig.width = 10}
+par(mfrow = c(2, 2))
+von_mises <- gof_ds(veldleeuwerik_dist_hr)
+von_mises
+plot(veldleeuwerik_dist_hr, pdf = TRUE, showpoints = TRUE)
+
+# Plot separately for covariate categories
+plot(veldleeuwerik_dist_hr, pdf = TRUE, showpoints = FALSE,
+   subset = openheid == "HOL",
+   main = "HOL", pl.col = alpha("green", 0.5))
+add.df.covar.line(veldleeuwerik_dist_hr, lwd = 3, lty = 1, col = "green",
+  data = data.frame(openheid = "HOL"), pdf = TRUE)
+
+plot(veldleeuwerik_dist_hr, pdf = TRUE, showpoints = FALSE,
+   subset = openheid == "OL",
+   main = "OL", pl.col = alpha("red", 0.5))
+add.df.covar.line(veldleeuwerik_dist_hr, lwd = 3, lty = 1, col = "red",
+  data = data.frame(openheid = "OL"), pdf = TRUE)
+par(mfrow = c(1, 1))
+```
+
+## Resultaten
 
 ```{r}
-ddf.gof(test)
-plot(test, pdf = TRUE)
+summary(veldleeuwerik_dist_hr)
 ```
 
-### Distance

From 168473b39efb4b4901e67e753da148db1d943c67 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Thu, 14 Nov 2024 17:18:27 +0100
Subject: [PATCH 07/13] results ds models

---
 .../markdown/test_densiteitsmodellering.Rmd   | 55 ++++++++++++++++++-
 1 file changed, 53 insertions(+), 2 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index 85375c54..632ce220 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -415,7 +415,7 @@ Om te controleren of ons model goed past bij de data (“goodness of fit”) mak
 
 De fit is niet perfect, maar goed genoeg voor deze oefening.
 
-```{r, fig.width = 10}
+```{r, fig.heigth = 20}
 par(mfrow = c(2, 2))
 von_mises <- gof_ds(veldleeuwerik_dist_hr)
 von_mises
@@ -439,6 +439,57 @@ par(mfrow = c(1, 1))
 ## Resultaten
 
 ```{r}
-summary(veldleeuwerik_dist_hr)
+sum_dist_veldleeuwerik <- summary(veldleeuwerik_dist_hr)
+```
+
+### Detectiekans
+
+```{r}
+cbind(veldleeuwerik_dist_hr$ddf$data, 
+      "p(z)" = predict(veldleeuwerik_dist_hr, 
+                       se.fit = TRUE)$fitted, 
+      "standard error" = predict(veldleeuwerik_dist_hr, 
+                                 se.fit = TRUE)$se) %>%
+  select("openheid", "p(z)", "standard error") %>%
+  distinct() %>%
+  arrange(openheid) %>%
+  kable()
+```
+
+### Densiteiten
+
+```{r}
+abundance_veldleeuwerik <- sum_dist_veldleeuwerik$dht$individuals$N
+density_veldleeuwerik <- sum_dist_veldleeuwerik$dht$individuals$D
+```
+
+Abundantie:
+
+```{r}
+abundance_veldleeuwerik %>%
+  select(stratum = Label, abundantie = Estimate, ll = lcl, ul = ucl) %>%
+  kable(digits = 3)
 ```
 
+Densiteit:
+
+```{r, warning=FALSE}
+density_veldleeuwerik %>%
+  select(stratum = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
+  separate(stratum, into = c("regio", "openheid", "sbp"), sep = " - ") %>%
+  kable(digits = 3)
+```
+
+```{r}
+density_veldleeuwerik %>%
+  filter(Label != "Total") %>%
+  select(stratum = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
+  separate(stratum, into = c("regio", "openheid", "sbp"), sep = " - ") %>%
+  mutate(stratum = paste(openheid, sbp, sep = " - ")) %>%
+  ggplot(aes(x = openheid, y = densiteit, colour = sbp)) +
+    geom_point(position = position_dodge(width = 0.5)) +
+    geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2,
+                  position = position_dodge(width = 0.5)) +
+    labs(y = "aantal broedparen per 100 ha", x = "") +
+    facet_wrap(~regio, ncol = 2, scales = "free")
+```

From c9a30546830c311c228830b6e93d0b9452a0dd82 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Thu, 14 Nov 2024 17:28:21 +0100
Subject: [PATCH 08/13] add conclusion

---
 .../markdown/test_densiteitsmodellering.Rmd   | 29 ++++++++++++-------
 1 file changed, 18 insertions(+), 11 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index 632ce220..c2828acf 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -207,11 +207,12 @@ veldleeuwerik_2024_df %>%
 ```
 
 # Distance sampling
+## Test workflow
 
 We schatten de detectiekans en densiteiten via distance sampling.
 We beschouwen voor de formulatie van de detectiefunctie de half-normal en hazard-rate sleutelfuncties die we laten afhangen van `openheid`.
 
-## Data preparatie
+### Data preparatie
 
 We hebben een kolom `object` nodig: volgnummer voor elke waarneming; een kolom `size`: aantal broedparen per waarneming (overal 1); een kolom `distance`: afstand tot elke waarneming; kolommen met covariaten (hier 1 variabele): `openheid`.
 
@@ -232,7 +233,7 @@ steekproefkader <- st_read(
             "steekproefkader_mbag_mas.gpkg"))
 
 # Buffer and area
-## Calculate area per region
+### Calculate area per region
 region_sf <- steekproefkader %>%
   mutate(regio = ifelse(grepl("\\sleemstreek$", regio),
                         "Leemstreek", regio)) %>%
@@ -243,7 +244,7 @@ region_sf <- steekproefkader %>%
   mutate(Area = as.numeric(st_area(geom)) / 1e6) %>%
   select(regio, Area, everything())
 
-## Calculate area for Flanders
+### Calculate area for Flanders
 flanders_sf <- steekproefkader %>%
   st_buffer(dist = 300) %>%
   summarise(geom = st_union(geom)) %>%
@@ -252,7 +253,7 @@ flanders_sf <- steekproefkader %>%
          Area = as.numeric(st_area(geom)) / 1e6) %>%
   select(regio, Area, everything())
 
-## Calculate area per stratum
+### Calculate area per stratum
 strata_sf <- steekproefkader %>%
   mutate(regio = ifelse(grepl("\\sleemstreek$", regio),
                         "Leemstreek", regio)) %>%
@@ -377,7 +378,7 @@ Ten slotte zet de conversiefactor de afstanden om van meter naar vierkante kilom
 conversion_factor <- convert_units("meter", NULL, "Square kilometer")
 ```
 
-## Model specificatie
+### Model specificatie
 
 We beschouwen voor de formulatie van de detectiefunctie de half-normal en hazard-rate sleutelfuncties die we laten afhangen van `openheid`.
 
@@ -395,7 +396,7 @@ veldleeuwerik_dist_hr <- ds(data = veldleeuwerik_distance, key = "hr",
   sample_table = sample_table, obs_table = obs_table)
 ```
 
-## Model selectie
+### Model selectie
 
 We vergelijken de AIC van de verschillende modellen.
 Als het verschil tussen AIC’s kleiner is dan 2, kiezen we het eenvoudigste van deze modellen (model met minste parameters).
@@ -409,7 +410,7 @@ summarize_ds_models(
   kable()
 ```
 
-## Model fit
+### Model fit
 
 Om te controleren of ons model goed past bij de data (“goodness of fit”) maken we een Q-Q plot waarbij de cumulatieve distributiefunctie van de gefitte detectiefunctie (CDF) wordt vergeleken met de distributie van de data (EDF). De Cramér-von Mises test kwantificeert de informatie van de Q-Q-plot door te testen of punten van de EDF en CDF uit dezelfde verdeling komen. Voor studies met punttransecten geven kansdichtheidsfunctieplots ook een goed idee over model fit.
 
@@ -436,13 +437,13 @@ add.df.covar.line(veldleeuwerik_dist_hr, lwd = 3, lty = 1, col = "red",
 par(mfrow = c(1, 1))
 ```
 
-## Resultaten
+### Resultaten
 
 ```{r}
 sum_dist_veldleeuwerik <- summary(veldleeuwerik_dist_hr)
 ```
 
-### Detectiekans
+Detectiekans:
 
 ```{r}
 cbind(veldleeuwerik_dist_hr$ddf$data, 
@@ -456,8 +457,6 @@ cbind(veldleeuwerik_dist_hr$ddf$data,
   kable()
 ```
 
-### Densiteiten
-
 ```{r}
 abundance_veldleeuwerik <- sum_dist_veldleeuwerik$dht$individuals$N
 density_veldleeuwerik <- sum_dist_veldleeuwerik$dht$individuals$D
@@ -493,3 +492,11 @@ density_veldleeuwerik %>%
     labs(y = "aantal broedparen per 100 ha", x = "") +
     facet_wrap(~regio, ncol = 2, scales = "free")
 ```
+
+### Conclusie en vragen
+
+De workflow werkt zoals verwacht.
+
+- Hoe krijgen we confidence limits voor de detectiekans?
+- We moeten het model ook fitten voor densiteitsschattingen per regio. Is de detectiekans hetzelfde? We verwachten van wel.
+- Hebben we correct met de maxima kunnen werken?

From 37ca5f45bce6ed727ae5200d51d1a4ab61a36581 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Thu, 14 Nov 2024 18:05:11 +0100
Subject: [PATCH 09/13] structue

---
 source/markdown/test_densiteitsmodellering.Rmd | 14 +++++++++++++-
 1 file changed, 13 insertions(+), 1 deletion(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index c2828acf..51026911 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -497,6 +497,18 @@ density_veldleeuwerik %>%
 
 De workflow werkt zoals verwacht.
 
-- Hoe krijgen we confidence limits voor de detectiekans?
+- Hoe krijgen we confidence limits voor de detectiekans? Wat is de distributie van detectiekans?
 - We moeten het model ook fitten voor densiteitsschattingen per regio. Is de detectiekans hetzelfde? We verwachten van wel.
 - Hebben we correct met de maxima kunnen werken?
+
+## Distributie detectiekans
+
+> zie theorie
+
+## Densiteitsschattingen op regio en Vlaams niveau
+
+zelfde model maar met andere region sample en obs tabellen
+
+## Densiteitsschattingen o.b.v. alle telperiodes
+
+effort = 4

From e203dc6cf084eca9e6b4b7762588fa5826904165 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Mon, 18 Nov 2024 12:38:30 +0100
Subject: [PATCH 10/13] explore methods

---
 .../markdown/test_densiteitsmodellering.Rmd   | 181 +++++++++++++++++-
 1 file changed, 178 insertions(+), 3 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index 51026911..4515acc0 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -336,7 +336,7 @@ for (region in sort(unique(strata_sf$regio))) {
 }
 ```
 
-We willen een densiteitsschatting voor elk stratum. We maken daarom dataframes om de abundanties en densiteiten te berekenen. `Region_table` met de oppervlaktes per stratum (regio-openheid-sbp), `sample_table` geeft aan welke plots in welke strata voorkomen en welke effort is geleverd (effort is gelijk aan het aantal beschouwde telperiodes per jaar) en `obs_table` die aangeeft welke waarnemingen in welke plots en strata zitten. We selecteren enkel de data per telpunt in de periode met het meeste waarnemingen.
+We willen een densiteitsschatting voor elk stratum. We maken daarom dataframes om de abundanties en densiteiten te berekenen. `region_table` met de oppervlaktes per stratum (regio-openheid-sbp), `sample_table` geeft aan welke plots in welke strata voorkomen en welke effort is geleverd (effort is gelijk aan het aantal beschouwde telperiodes per jaar) en `obs_table` die aangeeft welke waarnemingen in welke plots en strata zitten. We selecteren enkel de data per telpunt in de periode met het meeste waarnemingen.
 
 ```{r}
 design <- read_csv(
@@ -473,6 +473,10 @@ abundance_veldleeuwerik %>%
 Densiteit:
 
 ```{r, warning=FALSE}
+density_veldleeuwerik_tot <- density_veldleeuwerik[
+  density_veldleeuwerik$Label == "Total", "Estimate"
+  ]
+
 density_veldleeuwerik %>%
   select(stratum = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
   separate(stratum, into = c("regio", "openheid", "sbp"), sep = " - ") %>%
@@ -493,6 +497,10 @@ density_veldleeuwerik %>%
     facet_wrap(~regio, ncol = 2, scales = "free")
 ```
 
+Voor de Weidestreek zijn we niet van plan densiteiten te rapporteren op stratumniveau.
+Voor deze oefening zijn deze hier toch getoond.
+Voor de Weidestreek zullen we enkel op niveau van de landbouwstreek rapporteren (zie verder).
+
 ### Conclusie en vragen
 
 De workflow werkt zoals verwacht.
@@ -506,9 +514,176 @@ De workflow werkt zoals verwacht.
 > zie theorie
 
 ## Densiteitsschattingen op regio en Vlaams niveau
+### Model specificatie
+
+We fitten het finale tabel als voordien, maar we maken de `region_table`, `sample_table` en `obs_table` aan op regio-niveau in plaats van op stratum-niveau.
+
+```{r}
+# Area per region
+region_table2 <- region_sf %>%
+  st_drop_geometry() %>%
+  select(Region.Label = regio, Area)
+
+# Effort per sample en plot per region
+effort <- 1
+
+sample_table2 <- design %>%
+  distinct(pointid, regio, openheid = openheid_klasse, sbp) %>%
+  mutate(Effort = effort) %>%
+  select(Sample.Label = pointid, Region.Label = regio, Effort)
+
+# Observations per plots en regions
+obs_table2 <- veldleeuwerik_2024_df %>%
+  group_by(periode_in_jaar,  plotnaam) %>%
+  mutate(n = n()) %>%
+  slice_max(order_by = n, n = 1) %>%
+  ungroup() %>%
+  select(object = oid, Region.Label = regio, Sample.Label = plotnaam)
+
+# Distance sampling
+veldleeuwerik_dist_hr2 <- ds(data = veldleeuwerik_distance, key = "hr", 
+  formula = ~openheid, adjustment = NULL, truncation = 300, 
+  transect = "point", dht_group = FALSE, 
+  convert_units = conversion_factor, region_table = region_table2, 
+  sample_table = sample_table2, obs_table = obs_table2)
+```
+
+### Resultaten
+
+```{r}
+sum_dist_veldleeuwerik2 <- summary(veldleeuwerik_dist_hr2)
+```
+
+Detectiekans:
+
+```{r}
+cbind(veldleeuwerik_dist_hr2$ddf$data, 
+      "p(z)" = predict(veldleeuwerik_dist_hr2, 
+                       se.fit = TRUE)$fitted, 
+      "standard error" = predict(veldleeuwerik_dist_hr2, 
+                                 se.fit = TRUE)$se) %>%
+  select("openheid", "p(z)", "standard error") %>%
+  distinct() %>%
+  arrange(openheid) %>%
+  kable()
+```
+
+```{r}
+density_veldleeuwerik2 <- sum_dist_veldleeuwerik2$dht$individuals$D
+```
+
+Densiteit:
 
-zelfde model maar met andere region sample en obs tabellen
+```{r, warning=FALSE}
+density_veldleeuwerik2_tot <- density_veldleeuwerik2[
+  density_veldleeuwerik2$Label == "Total", "Estimate"
+  ]
+
+density_veldleeuwerik2 %>%
+  select(regio = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
+  kable(digits = 3)
+```
+
+```{r}
+density_veldleeuwerik2 %>%
+  filter(Label != "Total") %>%
+  select(regio = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
+  ggplot(aes(x = regio, y = densiteit)) +
+    geom_point(size = 3) +
+    geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2, linewidth = 1) +
+    labs(y = "aantal broedparen per 100 ha", x = "")
+```
+
+De detectiekans is dezelfde als voordien.
+Voordien (bij het stratum-model) hadden we een onderschatting van de totale densiteit binnen het steekproefkader:`r density_veldleeuwerik_tot` t.o.v. `r density_veldleeuwerik2_tot` gemiddeld aantal broedparen per 100 ha.
+Dit is omdat we bij het stratum-model geografisch overlap hebben.
+
+Om correcte schattingen te krijgen op regio- en steekproefkader-niveau, kunnen we dus het finale model herfitten op aangepaste datasets.
 
 ## Densiteitsschattingen o.b.v. alle telperiodes
+### Model specificatie
+
+We fitten het regio-model opnieuw, maar deze keer gebruiken we alle data voor de densiteitsschattingen met `effort = 4`.
+Dit om uit te zoeken of het gebruik van de maxima t.o.v. alle data de verwachte uitkomst biedt.
+We verwachten:
+
+1. identieke detectiekans
+   - dezelfde data wordt gebruikt om de detectiekans te schatten als voordien
+2. lagere densiteitsschattingen
+   - door uitmiddeling sampling effort
+3. smallere betrouwbaarheidsintervallen densiteiten
+   - meer data wordt gebruikt
+
+```{r}
+# Effort per sample en plot per region
+effort <- length(unique(veldleeuwerik_2024_df$periode_in_jaar))
+
+sample_table3 <- design %>%
+  distinct(pointid, regio, openheid = openheid_klasse, sbp) %>%
+  mutate(Effort = effort) %>%
+  select(Sample.Label = pointid, Region.Label = regio, Effort)
+
+# Observations per plots en regions
+obs_table3 <- veldleeuwerik_2024_df %>%
+  select(object = oid, Region.Label = regio, Sample.Label = plotnaam)
+
+# Distance sampling
+veldleeuwerik_dist_hr3 <- ds(data = veldleeuwerik_distance, key = "hr", 
+  formula = ~openheid, adjustment = NULL, truncation = 300, 
+  transect = "point", dht_group = FALSE, 
+  convert_units = conversion_factor, region_table = region_table2, 
+  sample_table = sample_table3, obs_table = obs_table3)
+```
+
+### Resultaten
+
+```{r}
+sum_dist_veldleeuwerik3 <- summary(veldleeuwerik_dist_hr3)
+```
+
+Detectiekans:
+
+```{r}
+cbind(veldleeuwerik_dist_hr3$ddf$data, 
+      "p(z)" = predict(veldleeuwerik_dist_hr3, 
+                       se.fit = TRUE)$fitted, 
+      "standard error" = predict(veldleeuwerik_dist_hr3, 
+                                 se.fit = TRUE)$se) %>%
+  select("openheid", "p(z)", "standard error") %>%
+  distinct() %>%
+  arrange(openheid) %>%
+  kable()
+```
+
+```{r}
+density_veldleeuwerik3 <- sum_dist_veldleeuwerik3$dht$individuals$D
+```
+
+Densiteit:
+
+```{r, warning=FALSE}
+density_veldleeuwerik3 %>%
+  select(regio = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
+  kable(digits = 3)
+```
+
+```{r}
+density_veldleeuwerik3 %>%
+  filter(Label != "Total") %>%
+  mutate(Dataset = "totaal (effort = 4)") %>%
+  bind_rows(
+    density_veldleeuwerik2 %>%
+      filter(Label != "Total") %>%
+      mutate(Dataset = "maxima (effort = 1)")) %>%
+  select(regio = Label, Dataset, densiteit = Estimate, ll = lcl, ul = ucl) %>%
+  ggplot(aes(x = regio, y = densiteit, colour = Dataset)) +
+    geom_point(size = 3, position = position_dodge(width = 0.5)) +
+    geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2, linewidth = 1,
+                  position = position_dodge(width = 0.5)) +
+    labs(y = "aantal broedparen per 100 ha", x = "")
+```
 
-effort = 4
+Aan alle verwachtingen is voldaan.
+Dit geeft aan dat we op correcte wijze de detectiekans kunnen schatten o.b.v. de volledige dataset en de densiteit o.b.v. een subset (maxima) of de volledige dataset.
+We zien echter wel dat de schattingen en onzekerheden tussen beide methodes sterk verschillend zijn.
+De keuze van ene of de andere methode moet onderbouwd worden o.b.v. ecologische kennis.

From 36314ecf43660ec68105c6d09910bc3d1ae8093c Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Mon, 18 Nov 2024 12:58:25 +0100
Subject: [PATCH 11/13] coding style

---
 .../markdown/test_densiteitsmodellering.Rmd   | 383 +++++++++++-------
 1 file changed, 231 insertions(+), 152 deletions(-)

diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index 4515acc0..e4e62831 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -81,12 +81,16 @@ mas_data_clean <- tar_read("mas_data_clean", store = targets_store)
 
 # Select data Veldleeuwerik 2024
 veldleeuwerik_2024_df <- mas_data_clean %>%
-  filter(naam == "Veldleeuwerik",
-         jaar == 2024,
-         wrntype > 0) %>%
-  mutate(aantal = 1,
-         regio = ifelse(grepl("\\sleemstreek$", regio), "Leemstreek", regio),
-         stratum = paste(openheid_klasse, sbp, sep = " - ")) %>%
+  filter(
+    naam == "Veldleeuwerik",
+    jaar == 2024,
+    wrntype > 0
+  ) %>%
+  mutate(
+    aantal = 1,
+    regio = ifelse(grepl("\\sleemstreek$", regio), "Leemstreek", regio),
+    stratum = paste(openheid_klasse, sbp, sep = " - ")
+  ) %>%
   st_drop_geometry()
 ```
 
@@ -96,9 +100,9 @@ We kijken naar het aantal broedparen per regio.
 # Visualise
 veldleeuwerik_2024_df %>%
   ggplot(aes(x = periode_in_jaar, fill = wrntype)) +
-    geom_bar() +
-    facet_wrap(~regio, scales = "free") +
-    labs(x = "Telperiode", y = "Aantal broedparen", fill = "broedcode")
+  geom_bar() +
+  facet_wrap(~regio, scales = "free") +
+  labs(x = "Telperiode", y = "Aantal broedparen", fill = "broedcode")
 ```
 
 En de per stratum.
@@ -107,9 +111,9 @@ En de per stratum.
 # Visualise
 veldleeuwerik_2024_df %>%
   ggplot(aes(x = periode_in_jaar, fill = stratum)) +
-    geom_bar() +
-    facet_wrap(~regio, scales = "free") +
-    labs(x = "Telperiode", y = "Aantal broedparen", fill = "Stratum")
+  geom_bar() +
+  facet_wrap(~regio, scales = "free") +
+  labs(x = "Telperiode", y = "Aantal broedparen", fill = "Stratum")
 ```
 
 Alle telperiodes bevinden zich binnen de datumgrenzen.
@@ -118,39 +122,68 @@ We nemen alle telperiode dus mee voor de densiteitsschattingen.
 ```{r}
 # Load breeding dates file
 datumgrenzen_df <- read_csv2(
-  file.path(mbag_dir, "data", "SOVON",
-            "Interpretatie_Criteria_Broedvogels_v2.csv"))
+  file.path(
+    mbag_dir, "data", "SOVON",
+    "Interpretatie_Criteria_Broedvogels_v2.csv"
+  )
+)
 
 # Tidy dataframe
 trans_vec <- c("jan", "feb", "mrt", "apr", "mei", "jun", "jul", "aug")
 
 datumgrenzen_df2 <- datumgrenzen_df %>%
   mutate(across(Datum_begin:`Datum eind_oud`, ~ gsub("-", "/", .x))) %>%
-  separate(Datum_begin, into = c("Datum_begin_dag", "Datum_begin_maand"), 
-           sep = "/") %>%
-  separate(Datum_eind, into = c("Datum_eind_dag", "Datum_eind_maand"), 
-           sep = "/") %>%
-  separate(`Datum begin_oud`, into = c("Datum_begin_dag_oud", 
-                                       "Datum_begin_maand_oud"), 
-           sep = "/") %>%
-  separate(`Datum eind_oud`, into = c("Datum_eind_dag_oud", 
-                                      "Datum_eind_maand_oud"), 
-           sep = "/") %>%
-  mutate(Datum_begin_maand = match(Datum_begin_maand, trans_vec),
-         datum_begin = paste(Datum_begin_maand, Datum_begin_dag, sep = "-")) %>%
-  mutate(Datum_eind_maand = match(Datum_eind_maand, trans_vec),
-         datum_eind = paste(Datum_eind_maand, Datum_eind_dag, sep = "-")) %>%
-  mutate(Datum_begin_maand_oud = match(Datum_begin_maand_oud, trans_vec),
-         datum_begin_oud = paste(Datum_begin_maand_oud, Datum_begin_dag_oud, 
-                                 sep = "-")) %>%
-  mutate(Datum_eind_maand_oud = match(Datum_eind_maand_oud, trans_vec),
-         datum_eind_oud = paste(Datum_eind_maand_oud, Datum_eind_dag_oud,
-                                sep = "-")) %>%
-  mutate(across(datum_begin:datum_eind_oud, 
-                ~ ifelse(.x == "NA-NA", NA, .x))) %>%
-  select(id = Id, soort = Soort, datum_begin, datum_eind, datum_begin_oud, 
-         datum_eind_oud, broedcode, broedcode_oud = bc_oud, 
-         wnm_vereist = `Waarnemingen verreist`, Opmerking)
+  separate(Datum_begin,
+    into = c("Datum_begin_dag", "Datum_begin_maand"),
+    sep = "/"
+  ) %>%
+  separate(Datum_eind,
+    into = c("Datum_eind_dag", "Datum_eind_maand"),
+    sep = "/"
+  ) %>%
+  separate(`Datum begin_oud`,
+    into = c(
+      "Datum_begin_dag_oud",
+      "Datum_begin_maand_oud"
+    ),
+    sep = "/"
+  ) %>%
+  separate(`Datum eind_oud`,
+    into = c(
+      "Datum_eind_dag_oud",
+      "Datum_eind_maand_oud"
+    ),
+    sep = "/"
+  ) %>%
+  mutate(
+    Datum_begin_maand = match(Datum_begin_maand, trans_vec),
+    datum_begin = paste(Datum_begin_maand, Datum_begin_dag, sep = "-")
+  ) %>%
+  mutate(
+    Datum_eind_maand = match(Datum_eind_maand, trans_vec),
+    datum_eind = paste(Datum_eind_maand, Datum_eind_dag, sep = "-")
+  ) %>%
+  mutate(
+    Datum_begin_maand_oud = match(Datum_begin_maand_oud, trans_vec),
+    datum_begin_oud = paste(Datum_begin_maand_oud, Datum_begin_dag_oud,
+      sep = "-"
+    )
+  ) %>%
+  mutate(
+    Datum_eind_maand_oud = match(Datum_eind_maand_oud, trans_vec),
+    datum_eind_oud = paste(Datum_eind_maand_oud, Datum_eind_dag_oud,
+      sep = "-"
+    )
+  ) %>%
+  mutate(across(
+    datum_begin:datum_eind_oud,
+    ~ ifelse(.x == "NA-NA", NA, .x)
+  )) %>%
+  select(
+    id = Id, soort = Soort, datum_begin, datum_eind, datum_begin_oud,
+    datum_eind_oud, broedcode, broedcode_oud = bc_oud,
+    wnm_vereist = `Waarnemingen verreist`, Opmerking
+  )
 
 # MAS dates
 r1_start <- "04-01"
@@ -169,24 +202,23 @@ datumgrenzen_mas <- datumgrenzen_df2 %>%
     jaar = 2024,
     datum_begin2 = ymd(paste(jaar, datum_begin, sep = "-")),
     datum_eind2 = ymd(paste(jaar, datum_eind, sep = "-")),
-    R1 = ifelse(ymd(paste(jaar, r1_start, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2) &
-                ymd(paste(jaar, r1_stop, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2), TRUE, FALSE),
-    R2 = ifelse(ymd(paste(jaar, r2_start, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2) &
-                ymd(paste(jaar, r2_stop, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2), TRUE, FALSE),
-    R3 = ifelse(ymd(paste(jaar, r3_start, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2) &
-                ymd(paste(jaar, r3_stop, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2), TRUE, FALSE),
-    R4 = ifelse(ymd(paste(jaar, r4_start, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2) &
-                ymd(paste(jaar, r4_stop, sep = "-")) %within% 
-                  interval(datum_begin2, datum_eind2), TRUE, FALSE)
-    
-    ) %>%
+    R1 = ifelse(ymd(paste(jaar, r1_start, sep = "-")) %within%
+      interval(datum_begin2, datum_eind2) &
+      ymd(paste(jaar, r1_stop, sep = "-")) %within%
+        interval(datum_begin2, datum_eind2), TRUE, FALSE),
+    R2 = ifelse(ymd(paste(jaar, r2_start, sep = "-")) %within%
+      interval(datum_begin2, datum_eind2) &
+      ymd(paste(jaar, r2_stop, sep = "-")) %within%
+        interval(datum_begin2, datum_eind2), TRUE, FALSE),
+    R3 = ifelse(ymd(paste(jaar, r3_start, sep = "-")) %within%
+      interval(datum_begin2, datum_eind2) &
+      ymd(paste(jaar, r3_stop, sep = "-")) %within%
+        interval(datum_begin2, datum_eind2), TRUE, FALSE),
+    R4 = ifelse(ymd(paste(jaar, r4_start, sep = "-")) %within%
+      interval(datum_begin2, datum_eind2) &
+      ymd(paste(jaar, r4_stop, sep = "-")) %within%
+        interval(datum_begin2, datum_eind2), TRUE, FALSE)
+  ) %>%
   select(id, soort, datum_begin, datum_eind, R1, R2, R3, R4)
 
 # Show results Veldleeuwerik
@@ -201,9 +233,9 @@ Ten slotte kijken we naar de afstanden.
 # Visualise
 veldleeuwerik_2024_df %>%
   ggplot(aes(x = distance2plot, fill = stratum)) +
-    geom_histogram() +
-    facet_wrap(~regio, scales = "free") +
-    labs(x = "Afstand (m)", y = "Aantal broedparen", fill = "Stratum")
+  geom_histogram() +
+  facet_wrap(~regio, scales = "free") +
+  labs(x = "Afstand (m)", y = "Aantal broedparen", fill = "Stratum")
 ```
 
 # Distance sampling
@@ -218,10 +250,12 @@ We hebben een kolom `object` nodig: volgnummer voor elke waarneming; een kolom `
 
 ```{r}
 veldleeuwerik_distance <- veldleeuwerik_2024_df %>%
-  select(object = oid,
-         size = aantal,
-         distance = distance2plot,
-         openheid = openheid_klasse)
+  select(
+    object = oid,
+    size = aantal,
+    distance = distance2plot,
+    openheid = openheid_klasse
+  )
 ```
 
 Om abundanties te schatten, moeten we de oppervlakte van de strata berekenen voor extrapolatie. We trekken een buffer van 300 m rond het steekproefkader.
@@ -229,14 +263,18 @@ Om abundanties te schatten, moeten we de oppervlakte van de strata berekenen voo
 ```{r}
 # Lees steekproefkader in
 steekproefkader <- st_read(
-  file.path(mbag_dir, "data", "steekproefkaders",
-            "steekproefkader_mbag_mas.gpkg"))
+  file.path(
+    mbag_dir, "data", "steekproefkaders",
+    "steekproefkader_mbag_mas.gpkg"
+  )
+)
 
 # Buffer and area
 ### Calculate area per region
 region_sf <- steekproefkader %>%
   mutate(regio = ifelse(grepl("\\sleemstreek$", regio),
-                        "Leemstreek", regio)) %>%
+    "Leemstreek", regio
+  )) %>%
   st_buffer(dist = 300) %>%
   group_by(regio) %>%
   summarise(geom = st_union(geom)) %>%
@@ -249,14 +287,17 @@ flanders_sf <- steekproefkader %>%
   st_buffer(dist = 300) %>%
   summarise(geom = st_union(geom)) %>%
   ungroup() %>%
-  mutate(regio = "Flanders",
-         Area = as.numeric(st_area(geom)) / 1e6) %>%
+  mutate(
+    regio = "Flanders",
+    Area = as.numeric(st_area(geom)) / 1e6
+  ) %>%
   select(regio, Area, everything())
 
 ### Calculate area per stratum
 strata_sf <- steekproefkader %>%
   mutate(regio = ifelse(grepl("\\sleemstreek$", regio),
-                        "Leemstreek", regio)) %>%
+    "Leemstreek", regio
+  )) %>%
   st_buffer(dist = 300) %>%
   group_by(regio, "openheid" = openheid_klasse, sbp) %>%
   summarise(geom = st_union(geom)) %>%
@@ -278,9 +319,9 @@ flanders_sf %>%
 # Visualisation
 flanders_sf %>%
   ggplot() +
-    geom_sf(aes(fill = regio)) +
-    labs(fill = "") +
-    theme(legend.position = "bottom")
+  geom_sf(aes(fill = regio)) +
+  labs(fill = "") +
+  theme(legend.position = "bottom")
 ```
 
 Per landbouwregio:
@@ -299,9 +340,9 @@ Als we de oppervlakte optellen komen we dezelfde oppervlakte uit als volledig Vl
 # Visualisation
 region_sf %>%
   ggplot() +
-    geom_sf(aes(fill = regio)) +
-    labs(fill = "") +
-    theme(legend.position = "bottom")
+  geom_sf(aes(fill = regio)) +
+  labs(fill = "") +
+  theme(legend.position = "bottom")
 ```
 
 Per stratum:
@@ -328,10 +369,10 @@ for (region in sort(unique(strata_sf$regio))) {
     filter(regio == region) %>%
     mutate(stratum = paste(openheid, sbp, sep = " - ")) %>%
     ggplot() +
-      geom_sf(aes(fill = stratum)) +
-      facet_wrap(~regio, ncol = 2) +
-      labs(fill = "") +
-      theme(legend.position = "bottom")
+    geom_sf(aes(fill = stratum)) +
+    facet_wrap(~regio, ncol = 2) +
+    labs(fill = "") +
+    theme(legend.position = "bottom")
   print(p)
 }
 ```
@@ -340,8 +381,11 @@ We willen een densiteitsschatting voor elk stratum. We maken daarom dataframes o
 
 ```{r}
 design <- read_csv(
-  file.path(mbag_dir, "data", "steekproefkaders",
-            "steekproef_avimap_mbag_mas.csv"))  %>%
+  file.path(
+    mbag_dir, "data", "steekproefkaders",
+    "steekproef_avimap_mbag_mas.csv"
+  )
+) %>%
   mutate(regio = ifelse(grepl("\\sleemstreek$", regio), "Leemstreek", regio))
 
 # Oppervlakte per stratum
@@ -356,18 +400,22 @@ effort <- 1
 
 sample_table <- design %>%
   distinct(pointid, regio, openheid = openheid_klasse, sbp) %>%
-  mutate(Region.Label = paste(regio, openheid, sbp, sep = " - "),
-         Effort = effort) %>%
+  mutate(
+    Region.Label = paste(regio, openheid, sbp, sep = " - "),
+    Effort = effort
+  ) %>%
   select(Sample.Label = pointid, Region.Label, Effort)
 
 # Waarnemingen per plots en strata
 obs_table <- veldleeuwerik_2024_df %>%
-  group_by(periode_in_jaar,  plotnaam) %>%
+  group_by(periode_in_jaar, plotnaam) %>%
   mutate(n = n()) %>%
   slice_max(order_by = n, n = 1) %>%
   ungroup() %>%
-  mutate(Region.Label = paste(regio, openheid = openheid_klasse, sbp,
-                              sep = " - ")) %>%
+  mutate(Region.Label = paste(regio,
+    openheid = openheid_klasse, sbp,
+    sep = " - "
+  )) %>%
   select(object = oid, Region.Label, Sample.Label = plotnaam)
 ```
 
@@ -383,17 +431,21 @@ conversion_factor <- convert_units("meter", NULL, "Square kilometer")
 We beschouwen voor de formulatie van de detectiefunctie de half-normal en hazard-rate sleutelfuncties die we laten afhangen van `openheid`.
 
 ```{r}
-veldleeuwerik_dist_hn <- ds(data = veldleeuwerik_distance, key = "hn", 
-  formula = ~openheid, adjustment = NULL, truncation = 300, 
-  transect = "point", dht_group = FALSE, 
-  convert_units = conversion_factor, region_table = region_table, 
-  sample_table = sample_table, obs_table = obs_table)
+veldleeuwerik_dist_hn <- ds(
+  data = veldleeuwerik_distance, key = "hn",
+  formula = ~openheid, adjustment = NULL, truncation = 300,
+  transect = "point", dht_group = FALSE,
+  convert_units = conversion_factor, region_table = region_table,
+  sample_table = sample_table, obs_table = obs_table
+)
 
-veldleeuwerik_dist_hr <- ds(data = veldleeuwerik_distance, key = "hr", 
-  formula = ~openheid, adjustment = NULL, truncation = 300, 
-  transect = "point", dht_group = FALSE, 
-  convert_units = conversion_factor, region_table = region_table, 
-  sample_table = sample_table, obs_table = obs_table)
+veldleeuwerik_dist_hr <- ds(
+  data = veldleeuwerik_distance, key = "hr",
+  formula = ~openheid, adjustment = NULL, truncation = 300,
+  transect = "point", dht_group = FALSE,
+  convert_units = conversion_factor, region_table = region_table,
+  sample_table = sample_table, obs_table = obs_table
+)
 ```
 
 ### Model selectie
@@ -404,9 +456,9 @@ Het Hazard-rate model scoort beter.
 
 ```{r}
 summarize_ds_models(
-    veldleeuwerik_dist_hn,
-    veldleeuwerik_dist_hr
-  ) %>%
+  veldleeuwerik_dist_hn,
+  veldleeuwerik_dist_hr
+) %>%
   kable()
 ```
 
@@ -423,17 +475,25 @@ von_mises
 plot(veldleeuwerik_dist_hr, pdf = TRUE, showpoints = TRUE)
 
 # Plot separately for covariate categories
-plot(veldleeuwerik_dist_hr, pdf = TRUE, showpoints = FALSE,
-   subset = openheid == "HOL",
-   main = "HOL", pl.col = alpha("green", 0.5))
-add.df.covar.line(veldleeuwerik_dist_hr, lwd = 3, lty = 1, col = "green",
-  data = data.frame(openheid = "HOL"), pdf = TRUE)
-
-plot(veldleeuwerik_dist_hr, pdf = TRUE, showpoints = FALSE,
-   subset = openheid == "OL",
-   main = "OL", pl.col = alpha("red", 0.5))
-add.df.covar.line(veldleeuwerik_dist_hr, lwd = 3, lty = 1, col = "red",
-  data = data.frame(openheid = "OL"), pdf = TRUE)
+plot(veldleeuwerik_dist_hr,
+  pdf = TRUE, showpoints = FALSE,
+  subset = openheid == "HOL",
+  main = "HOL", pl.col = alpha("green", 0.5)
+)
+add.df.covar.line(veldleeuwerik_dist_hr,
+  lwd = 3, lty = 1, col = "green",
+  data = data.frame(openheid = "HOL"), pdf = TRUE
+)
+
+plot(veldleeuwerik_dist_hr,
+  pdf = TRUE, showpoints = FALSE,
+  subset = openheid == "OL",
+  main = "OL", pl.col = alpha("red", 0.5)
+)
+add.df.covar.line(veldleeuwerik_dist_hr,
+  lwd = 3, lty = 1, col = "red",
+  data = data.frame(openheid = "OL"), pdf = TRUE
+)
 par(mfrow = c(1, 1))
 ```
 
@@ -446,11 +506,14 @@ sum_dist_veldleeuwerik <- summary(veldleeuwerik_dist_hr)
 Detectiekans:
 
 ```{r}
-cbind(veldleeuwerik_dist_hr$ddf$data, 
-      "p(z)" = predict(veldleeuwerik_dist_hr, 
-                       se.fit = TRUE)$fitted, 
-      "standard error" = predict(veldleeuwerik_dist_hr, 
-                                 se.fit = TRUE)$se) %>%
+cbind(veldleeuwerik_dist_hr$ddf$data,
+  "p(z)" = predict(veldleeuwerik_dist_hr,
+    se.fit = TRUE
+  )$fitted,
+  "standard error" = predict(veldleeuwerik_dist_hr,
+    se.fit = TRUE
+  )$se
+) %>%
   select("openheid", "p(z)", "standard error") %>%
   distinct() %>%
   arrange(openheid) %>%
@@ -475,7 +538,7 @@ Densiteit:
 ```{r, warning=FALSE}
 density_veldleeuwerik_tot <- density_veldleeuwerik[
   density_veldleeuwerik$Label == "Total", "Estimate"
-  ]
+]
 
 density_veldleeuwerik %>%
   select(stratum = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
@@ -490,11 +553,13 @@ density_veldleeuwerik %>%
   separate(stratum, into = c("regio", "openheid", "sbp"), sep = " - ") %>%
   mutate(stratum = paste(openheid, sbp, sep = " - ")) %>%
   ggplot(aes(x = openheid, y = densiteit, colour = sbp)) +
-    geom_point(position = position_dodge(width = 0.5)) +
-    geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2,
-                  position = position_dodge(width = 0.5)) +
-    labs(y = "aantal broedparen per 100 ha", x = "") +
-    facet_wrap(~regio, ncol = 2, scales = "free")
+  geom_point(position = position_dodge(width = 0.5)) +
+  geom_errorbar(aes(ymin = ll, ymax = ul),
+    width = 0.2,
+    position = position_dodge(width = 0.5)
+  ) +
+  labs(y = "aantal broedparen per 100 ha", x = "") +
+  facet_wrap(~regio, ncol = 2, scales = "free")
 ```
 
 Voor de Weidestreek zijn we niet van plan densiteiten te rapporteren op stratumniveau.
@@ -534,18 +599,20 @@ sample_table2 <- design %>%
 
 # Observations per plots en regions
 obs_table2 <- veldleeuwerik_2024_df %>%
-  group_by(periode_in_jaar,  plotnaam) %>%
+  group_by(periode_in_jaar, plotnaam) %>%
   mutate(n = n()) %>%
   slice_max(order_by = n, n = 1) %>%
   ungroup() %>%
   select(object = oid, Region.Label = regio, Sample.Label = plotnaam)
 
 # Distance sampling
-veldleeuwerik_dist_hr2 <- ds(data = veldleeuwerik_distance, key = "hr", 
-  formula = ~openheid, adjustment = NULL, truncation = 300, 
-  transect = "point", dht_group = FALSE, 
-  convert_units = conversion_factor, region_table = region_table2, 
-  sample_table = sample_table2, obs_table = obs_table2)
+veldleeuwerik_dist_hr2 <- ds(
+  data = veldleeuwerik_distance, key = "hr",
+  formula = ~openheid, adjustment = NULL, truncation = 300,
+  transect = "point", dht_group = FALSE,
+  convert_units = conversion_factor, region_table = region_table2,
+  sample_table = sample_table2, obs_table = obs_table2
+)
 ```
 
 ### Resultaten
@@ -557,11 +624,14 @@ sum_dist_veldleeuwerik2 <- summary(veldleeuwerik_dist_hr2)
 Detectiekans:
 
 ```{r}
-cbind(veldleeuwerik_dist_hr2$ddf$data, 
-      "p(z)" = predict(veldleeuwerik_dist_hr2, 
-                       se.fit = TRUE)$fitted, 
-      "standard error" = predict(veldleeuwerik_dist_hr2, 
-                                 se.fit = TRUE)$se) %>%
+cbind(veldleeuwerik_dist_hr2$ddf$data,
+  "p(z)" = predict(veldleeuwerik_dist_hr2,
+    se.fit = TRUE
+  )$fitted,
+  "standard error" = predict(veldleeuwerik_dist_hr2,
+    se.fit = TRUE
+  )$se
+) %>%
   select("openheid", "p(z)", "standard error") %>%
   distinct() %>%
   arrange(openheid) %>%
@@ -577,7 +647,7 @@ Densiteit:
 ```{r, warning=FALSE}
 density_veldleeuwerik2_tot <- density_veldleeuwerik2[
   density_veldleeuwerik2$Label == "Total", "Estimate"
-  ]
+]
 
 density_veldleeuwerik2 %>%
   select(regio = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
@@ -589,9 +659,9 @@ density_veldleeuwerik2 %>%
   filter(Label != "Total") %>%
   select(regio = Label, densiteit = Estimate, ll = lcl, ul = ucl) %>%
   ggplot(aes(x = regio, y = densiteit)) +
-    geom_point(size = 3) +
-    geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2, linewidth = 1) +
-    labs(y = "aantal broedparen per 100 ha", x = "")
+  geom_point(size = 3) +
+  geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2, linewidth = 1) +
+  labs(y = "aantal broedparen per 100 ha", x = "")
 ```
 
 De detectiekans is dezelfde als voordien.
@@ -628,11 +698,13 @@ obs_table3 <- veldleeuwerik_2024_df %>%
   select(object = oid, Region.Label = regio, Sample.Label = plotnaam)
 
 # Distance sampling
-veldleeuwerik_dist_hr3 <- ds(data = veldleeuwerik_distance, key = "hr", 
-  formula = ~openheid, adjustment = NULL, truncation = 300, 
-  transect = "point", dht_group = FALSE, 
-  convert_units = conversion_factor, region_table = region_table2, 
-  sample_table = sample_table3, obs_table = obs_table3)
+veldleeuwerik_dist_hr3 <- ds(
+  data = veldleeuwerik_distance, key = "hr",
+  formula = ~openheid, adjustment = NULL, truncation = 300,
+  transect = "point", dht_group = FALSE,
+  convert_units = conversion_factor, region_table = region_table2,
+  sample_table = sample_table3, obs_table = obs_table3
+)
 ```
 
 ### Resultaten
@@ -644,11 +716,14 @@ sum_dist_veldleeuwerik3 <- summary(veldleeuwerik_dist_hr3)
 Detectiekans:
 
 ```{r}
-cbind(veldleeuwerik_dist_hr3$ddf$data, 
-      "p(z)" = predict(veldleeuwerik_dist_hr3, 
-                       se.fit = TRUE)$fitted, 
-      "standard error" = predict(veldleeuwerik_dist_hr3, 
-                                 se.fit = TRUE)$se) %>%
+cbind(veldleeuwerik_dist_hr3$ddf$data,
+  "p(z)" = predict(veldleeuwerik_dist_hr3,
+    se.fit = TRUE
+  )$fitted,
+  "standard error" = predict(veldleeuwerik_dist_hr3,
+    se.fit = TRUE
+  )$se
+) %>%
   select("openheid", "p(z)", "standard error") %>%
   distinct() %>%
   arrange(openheid) %>%
@@ -674,13 +749,17 @@ density_veldleeuwerik3 %>%
   bind_rows(
     density_veldleeuwerik2 %>%
       filter(Label != "Total") %>%
-      mutate(Dataset = "maxima (effort = 1)")) %>%
+      mutate(Dataset = "maxima (effort = 1)")
+  ) %>%
   select(regio = Label, Dataset, densiteit = Estimate, ll = lcl, ul = ucl) %>%
   ggplot(aes(x = regio, y = densiteit, colour = Dataset)) +
-    geom_point(size = 3, position = position_dodge(width = 0.5)) +
-    geom_errorbar(aes(ymin = ll, ymax = ul), width = 0.2, linewidth = 1,
-                  position = position_dodge(width = 0.5)) +
-    labs(y = "aantal broedparen per 100 ha", x = "")
+  geom_point(size = 3, position = position_dodge(width = 0.5)) +
+  geom_errorbar(aes(ymin = ll, ymax = ul),
+    width = 0.2, linewidth = 1,
+    position = position_dodge(width = 0.5)
+  ) +
+  labs(y = "aantal broedparen per 100 ha", x = "") +
+  theme(legend.position = "bottom")
 ```
 
 Aan alle verwachtingen is voldaan.

From c4d22b4d8d1a9626888fd9ba57da4d8ffa767ac3 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Mon, 18 Nov 2024 13:01:43 +0100
Subject: [PATCH 12/13] spelling check

---
 inst/nl_be.dic                                | 29 +++++++++++++++++++
 .../markdown/test_densiteitsmodellering.Rmd   |  6 ++--
 2 files changed, 32 insertions(+), 3 deletions(-)

diff --git a/inst/nl_be.dic b/inst/nl_be.dic
index e9bd95b7..55576bbe 100644
--- a/inst/nl_be.dic
+++ b/inst/nl_be.dic
@@ -45,6 +45,7 @@ Casteleynstraat
 Chi-kwadraattoetsingsgrootheid
 Core
 Core-termen
+Cramér-von
 Curruca
 DHMV
 DP
@@ -82,6 +83,7 @@ MBAG
 Maaibos
 Matkop
 Meetnetontwerp
+Mises
 NOVN00
 Natuur-
 OGC
@@ -92,6 +94,8 @@ Ol
 Parus
 Pl
 Poecile
+Q-Q
+Q-Q-plot
 QGIS
 R1
 R1-R4
@@ -136,25 +140,31 @@ WFS-verbinding
 WGS
 Westhoekreservaat
 \
+abundanties
 accepted
 add
 aeroway
 agro-ecologisch
 allocatieproporties
+and
 avimap
 avimapcodes
 batching
 beheren'
+branch
+branchen
 branching
 calculate
 canopy
 cf
 coden
 coding
+confidence
 confounders
 confounding
 control
 counts
+covariaten
 crs
 cumulative
 curruca
@@ -165,26 +175,37 @@ distance
 double
 dubbelgeteld
 dynamic
+effort
+estimation
 exclusiekaart
 exclusielaag
 final
+formulatie
+gefitte
 gelocaliseerd
 geojson-bestand
 gespecifieerd
 git2rdata
 globals
+goodness
 groups
+half-normal
 half-open
+hazard-rate
 height
+herfitten
 highway
+indensity
 indicatorfunctie
 inner
+kansdichtheidsfunctieplots
 kp
 landgebruiksamenstelling
 landuse
 layers
 leemstreek.Rmd
 lijn-
+limits
 ll'en
 lm.Rmd
 lm.Rproj
@@ -194,6 +215,7 @@ manual
 mapping
 markdown
 mas
+modelling
 monocotylen
 montanus
 multipolygon
@@ -220,8 +242,10 @@ probleemtelpunt
 proffessionele
 publicatieklaar
 puntenset
+punttransecten
 rastercel
 regex
+regio-openheid-sbp
 remove
 replicaten
 repo
@@ -233,11 +257,15 @@ size
 soortgrp
 soortnamen/auteurs
 sovon
+spatiaal
 spatiale
 spec
 statistical
 ste
+steekproefkader-niveau
 stratificatievariabelen
+stratum-model
+stratum-niveau
 stratumoppervlakte
 subspecies
 synonym
@@ -256,6 +284,7 @@ tertiare
 thinning
 toevoegen'
 topographic
+uitmiddeling
 uncut
 union
 variable
diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index e4e62831..e1d577be 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -47,7 +47,7 @@ targets_store <- here::here("source", "targets", "data_preparation", "_targets")
 
 # Achtergrond
 
-In kader van het meetnet agrarische soorten (MAS), worden vogels en Haas gemonitord volgens een vast protocol. Hierbij geven tellers waarnemingen in vanaf een telpunt tot 300 meter ver. Dit doen ze 4x per jaar. Doordat we de afstand van de teller tot de waarneming kennen, kunnen we de detectiekans per soort schatten en vervolgens de werkelijke densiteit berekenen (incl. onzekerheid).
+In kader van het meetnet agrarische soorten (MAS), worden vogels en Haas gemonitord volgens een vast protocol. Hierbij geven tellers waarnemingen in vanaf een telpunt tot 300 meter ver. Dit doen ze 4 x per jaar. Doordat we de afstand van de teller tot de waarneming kennen, kunnen we de detectiekans per soort schatten en vervolgens de werkelijke densiteit berekenen (incl. onzekerheid).
 
 Het meetnet dekt verschillende strata: landbouwregio’s, soortbeschermingsplan (SBP; binnen of buiten) en openheid landschap (OL: open landschap, HOL: half-open landschap).
 
@@ -452,7 +452,7 @@ veldleeuwerik_dist_hr <- ds(
 
 We vergelijken de AIC van de verschillende modellen.
 Als het verschil tussen AIC’s kleiner is dan 2, kiezen we het eenvoudigste van deze modellen (model met minste parameters).
-Het Hazard-rate model scoort beter.
+Het hazard-rate model scoort beter.
 
 ```{r}
 summarize_ds_models(
@@ -562,7 +562,7 @@ density_veldleeuwerik %>%
   facet_wrap(~regio, ncol = 2, scales = "free")
 ```
 
-Voor de Weidestreek zijn we niet van plan densiteiten te rapporteren op stratumniveau.
+Voor de Weidestreek zijn we niet van plan densiteiten te rapporteren op stratum-niveau.
 Voor deze oefening zijn deze hier toch getoond.
 Voor de Weidestreek zullen we enkel op niveau van de landbouwstreek rapporteren (zie verder).
 

From b59ce91618b22a0fa25c27a5d1453c198894f0b7 Mon Sep 17 00:00:00 2001
From: Ward Langeraert <ward.langeraert@inbo.be>
Date: Mon, 18 Nov 2024 14:14:09 +0100
Subject: [PATCH 13/13] add conclusion

---
 inst/nl_be.dic                                |  1 +
 .../markdown/test_densiteitsmodellering.Rmd   | 34 +++++++++++++++++++
 2 files changed, 35 insertions(+)

diff --git a/inst/nl_be.dic b/inst/nl_be.dic
index 55576bbe..ad2040e3 100644
--- a/inst/nl_be.dic
+++ b/inst/nl_be.dic
@@ -170,6 +170,7 @@ cumulative
 curruca
 cvvi
 day
+density
 design'
 distance
 double
diff --git a/source/markdown/test_densiteitsmodellering.Rmd b/source/markdown/test_densiteitsmodellering.Rmd
index e1d577be..34a38f39 100644
--- a/source/markdown/test_densiteitsmodellering.Rmd
+++ b/source/markdown/test_densiteitsmodellering.Rmd
@@ -766,3 +766,37 @@ Aan alle verwachtingen is voldaan.
 Dit geeft aan dat we op correcte wijze de detectiekans kunnen schatten o.b.v. de volledige dataset en de densiteit o.b.v. een subset (maxima) of de volledige dataset.
 We zien echter wel dat de schattingen en onzekerheden tussen beide methodes sterk verschillend zijn.
 De keuze van ene of de andere methode moet onderbouwd worden o.b.v. ecologische kennis.
+
+## Conclusie
+
+- We kunnen de detectiekans schatten o.b.v. de complete dataset en tegelijk de densiteit schatten o.b.v. een subset van de data.
+  - De keuze van de subset is afhankelijk van de ecologische context en moet met experts worden besproken.
+- We moeten twee modellen fitten om schattingen te krijgen om stratum-niveau, en op regio- en steekproefkader-niveau
+  - We fitten meerdere modellen (bv. op stratum-niveau) waarbij we het beste model voor de detectiefunctie selecteren (sleutelfunctie, covariaten ...).
+  - Het finale model wordt opnieuw gefit waarbij de dataset voor densiteitsschatting wordt aangepast (bv. op regio-niveau).
+
+Technisch gezien lijkt het dus dat we de methodiek van de pilootstudie eenvoudig kunnen toepassen voor de volledige MAS-data.
+In praktijk verwachten we nog altijd enkele uitdagingen:
+
+- Hoe kunnen we de analyses automatiseren als er in sommige strata weinig data aanwezig zijn?
+  - Door de inclusie van de Weidestreek (kleine steekproef) zullen interactietermen met regio waarschijnlijk niet vaak kunnen gefit worden.
+  - We kunnen eventueel strata met "te weinig data" weglaten uit de analyses en hiervan enkel de absolute aantallen rapporteren. Op regio- of steekproefniveau kunnen ze gewoon opgeteld worden bij de schattingen.
+    - Nood aan regels (bv. $n > 15$ per stratum).
+    - Automatisatie?
+  - Alternatief is enkel simpele detectiefuncties te fitten (zonder interactietermen). Of we kunnen regels verzinnen om simpele detectiefuncties te fitten als er te veel problemen bij de andere optreden (bv. voor zeldzame soorten is het $a priori$ al niet nuttig om interactietermen toe te voegen).
+
+```{r}
+design %>%
+  filter(regio == "Weidestreek") %>%
+  count(regio, openheid = openheid_klasse, sbp)
+```
+
+- Hoe automatiseren we detectiefunctie modelselectie in de targets pipeline?
+  - We kunnen bepaalde automatisatie inbouwen: kleinste AIC en als verschil kleiner is dan 2, kies dan het model met de minste parameters.
+  - De praktijk zal moeten uitwijzen of dit effectief correct verloopt. We kunnen de pipeline modelselectie plots en statistieken laten uitschrijven die we kunnen controleren in een apart bestand.
+  - Afhankelijk van de regels zullen we allicht meer of minder problemen krijgen bij model selectie. Enkel strata selecteren met veel data zal niet snel voor problemen zorgen, maar aan de andere kant wil je ook niet te veel data weglaten.
+
+In dit document onderzoeken we de technische mogelijkheden.
+Daarom kozen we voor Veldleeuwerik wat een heel algemene soort is.
+Bovenstaande uitdagingen zullen we op zeldzamere soorten moeten aftoetsen.
+Dit zal in een andere analyse gebeuren (mogelijks tijdens de implementatie van de pipeline) en niet hier.