In inbo/abvanalysis: Analysis of the Common Breeding Bird Survey
(PART) Methodologie {-}
Inzameling van de tellingen in het veld
Hier beschrijven we beknopt het veldwerk, voor de details verwijzen we naar @ABVmethode.
Steekproef {#s:steekproef}
De UTM 1 x 1 km hokken vormen de basis waaruit de steekproef getrokken is.
Bij de start van het meetnet hebben we de hokken opgedeeld in een aantal strata.
De regels zijn gebaseerd op het oppervlakteaandeel van een bepaald landgebruik op basis van de Biologische Waarderingkaart [@BWK].
- Landbouw: minstens 80% landbouw.
r strata$n[strata$description == "Landbouw"] hokken.
- Urbaan: minstens 80% urbaan.
r strata$n[strata$description == "Urbaan"] hokken.
- Bos: minstens 80% bos.
r strata$n[strata$description == "Bos"] hokken.
- Suburbaan: minstens 80% suburbaan.
r strata$n[strata$description == "Suburbaan"] hokken.
- Heide en duin: minstens 20% heide of duin.
r strata$n[strata$description == "Heide en duin"] hokken.
- Moeras en water: minstens 20% moeras en water.
r strata$n[strata$description == "Moeras en water"] hokken.
Uit deze set trekken we een aselecte, gestratificeerde steekproef van 1200 hokken waarbij zeldzamere habitats overbemonsterd worden.
We streven er naar om jaarlijks 300 hokken te bemonsteren in een driejarige rotatie.
De waarnemers mochten in het eerste jaar 300 hokken kiezen uit de set van 1200.
Deze set van hokken komen in principe opnieuw aan bod in jaren 4, 7, 10, ...
In jaar 2 kiezen ze 300 hokken uit de overgebleven 900 hokken.
Deze set komt opnieuw aan bod in jaren 5, 8, 11, ...
Tenslotte kiezen de waarnemers in het derde jaar een laatste set van 300 hokken uit de laatste 900 hokken.
Deze set hokken bemonsteren we in de jaren 3, 6, 9, 12, ...
Figuur \@ref(fig:inspanning) geeft de effectieve monitoringsinspanning weer.
In deze figuur hebben we de hokken gesorteerd volgens 1) het eerste jaar met gegevens, 2) het laatste jaar met gegevens, 3) het tweede jaar met gegevens, 4) het derde jaar met gegevens, ...
Hierdoor staan hokken met een meer gelijkende onderzoeksgeschiedenis dicht bij elkaar.
Merk op dat de driejarige cyclus voor de meest hokken wordt gerespecteerd.
Voor sommige hokken is de inspanning variabel, soms frequenter dan om de drie jaar, soms zit er meer tijd tussen.
Sommige hokken werden slechts in een of twee jaar onderzocht (fig. \@ref(fig:inspanning-hist)).
Voor het onderscheid tussen mogelijk bruikbaar en voorlopig niet bruikbaar verwijzen we naar §\@ref(s:relevant).
effort %>%
arrange(.data$year) %>%
group_by(.data$location) %>%
transmute(
.data$year, end = max(.data$year), year2 = lead(.data$year),
year3 = lead(.data$year2), year4 = lead(.data$year3),
year5 = lead(.data$year4), year6 = lead(.data$year5),
year7 = lead(.data$year6)
) %>%
slice(1) %>%
ungroup() %>%
arrange(
.data$year, .data$end, .data$year2, .data$year3, .data$year4, .data$year5,
.data$year6, .data$year7
) %>%
rowid_to_column("hok") %>%
transmute(
.data$location, hok = factor(.data$hok),
bruikbaar = ((.data$year - 2007) %/% 3) < ((end - 2007) %/% 3),
bruikbaar = factor(
.data$bruikbaar, levels = c(TRUE, FALSE),
labels = c("mogelijk", "voorlopig\nniet")
)
) %>%
inner_join(effort, by = "location") -> inspanning
ggplot(inspanning, aes(x = year, y = hok, group = hok, colour = bruikbaar)) +
geom_line(alpha = 0.2) +
geom_point(alpha = 0.5, size = ifelse(output_format == "html", 4, 1)) +
facet_wrap(~stratum, scales = "free_y") +
scale_x_continuous(breaks = pretty) +
theme(
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90),
axis.text.y = element_blank()
)
inspanning %>%
count(stratum, location, bruikbaar) %>%
ggplot(aes(x = n, fill = bruikbaar)) +
geom_histogram(binwidth = 1) +
facet_wrap(~stratum) +
scale_x_continuous(
"Aantal jaren waarin een hok onderzocht is", breaks = pretty,
limits = c(0, NA)
) +
ylab("Aantal hokken")
Steekproefeenheid
De steekproefeenheid bestaat uit een hok van 1 x 1 km.
Binnen dit hok worden zes punten vastgelegd in een vaste patroon (fig. \@ref(fig:patroon)).
Indien de punten in de praktijk niet bereikbaar zijn, mag de waarnemer ze verplaatsen naar het dichtstbijzijnde bereikbare punt.
De waarnemer documenteert deze wijziging zodat we in de toekomst steeds op dezelfde punten blijven waarnemen.
expand.grid(
x = seq(1 / 4, by = 1 / 2, length = 2),
y = seq(1 / 6, by = 1 / 3, length = 3)
) %>%
ggplot(aes(x = x, y = y)) +
geom_rect(
data = tribble(~x, ~y, ~xend, ~yend, 0, 0, 1, 1),
aes(xmin = x, ymin = y, xmax = xend, ymax = yend),
fill = NA, colour = inbo_steun_blauw
) +
geom_point() +
coord_fixed() +
theme_void()
Tellingen
In het jaar dat we een hok bemonsteren zal de waarnemer het hok in drie periodes bezoeken: 1 maart - 15 april, 16 april - 31 mei, 1 juni - 15 juli.
Tussen twee opeenvolgende bezoeken moet er minstens twee weken liggen.
Alle meetpunten van een hok worden op dezelfde dag onderzocht tussen zonsopgang en 4 uur na zonsopgang.
Op elk meetpunt telt de waarnemer gedurende 5 minuten het aantal volwassen vogels per soort.
Overvliegende groepen vogels worden hierbij niet meegeteld.
Verwerking van de gegevens
Selectie van relevante gegevens per soort {#s:relevant}
Hoewel we alle waargenomen vogelsoorten op een meetpunt noteren, toch zijn niet alle hokken even relevant voor een bepaalde soort.
Om de redenering op te bouwen, beschouwen we een fictieve soort die enkel op een afgelegen eiland voorkomt.
Alle vogeltellingen buiten dit eiland zijn uiteraard niet informatief voor wijzigingen in de aantallen van deze fictieve soort.
Buiten het eiland zullen de getelde aantallen immers per definitie nul zijn.
Wanneer we deze tellingen toch zouden in rekening brengen, dan zullen ze de globale trend afvlakken.
De globale trend is min of meer equivalent met een "gemiddelde" trend over de verschillende locaties.
Als een locatie steeds nul geeft, dan zijn lokaal de aantallen wiskundig "stabiel" (geen lokale trend).
Daarom is het noodzakelijk om de locaties die niet relevant zijn voor een soort buiten beschouwing te houden.
In Vlaanderen is de situatie niet zo zwart-wit als bij het voorbeeld met het afgelegen eiland.
Daarom hebben nood aan een set regels die we toepassen op de beschikbare gegevens.
- We beschouwen een hok als relevant voor een soort wanneer we deze soort tijdens minstens twee verschillende driejarige cycli waarnemen.
Op deze manier sluiten we hokken uit waar de soort nooit of slechts toevallig waargenomen is.
Hokken die slechts binnen één driejarige cyclus onderzocht zijn, kunnen we hierdoor voorlopig niet gebruiken.
Figuur \@ref(fig:potentieel) geeft een overzicht van het aantal hokken die gedurende minstens twee cycli onderzocht zijn.
- Binnen een hok is een meetpunt pas relevant wanneer we de soort er minstens eenmaal waarnemen.
- Een stratum is relevant voor een soort wanneer het minstens drie relevante hokken bevat.
- Op basis van de overblijvende gegevens berekenen we het (meetkundig[^geometrisch]) gemiddelde van de aantallen in elk van de drie perioden.
De periode met het hoogste gemiddelde gebruiken we als referentieperiode.
We houden enkel de periodes waarbij het gemiddelde minstens 15% van de referentieperiode bedraagt.
Op deze manier sluiten we periodes uit waarbij de soort zelden gezien wordt.
Bijvoorbeeld een trekvogel die pas in de loop van de tweede periode toekomt en waarvoor de eerste periode bijgevolg niet relevant is.
- We passen de regels voor relevante hokken en strata opnieuw toe.
Het negeren van een of twee irrelevante periodes kan er voor zorgen dat een hok niet langer relevant is.
[^geometrisch]: Het meetkundig gemiddelde of geometrisch gemiddelde van $n$ getallen wordt verkregen door de getallen met elkaar te vermenigvuldigen en vervolgens van het product de n^de^ machtswortel te nemen. $$m=\left(\prod {i=1}^{n}a{i}\right)^{1/n}$$
inspanning %>%
mutate(cyclus = (year - 2007) %/% 3) %>%
distinct(location, cyclus) %>%
count(location) %>%
filter(n >= 2) %>%
semi_join(x = inspanning, by = "location") %>%
count(year, stratum, name = "hokken") %>%
mutate(stratum = reorder(stratum, -hokken)) %>%
ggplot(aes(x = year, y = hokken, colour = stratum)) +
geom_line() +
scale_x_continuous(breaks = pretty) +
scale_y_continuous("Aantal potentieel relevante hokken", limits = c(0, NA)) +
theme(axis.title.x = element_blank())
Na het toepassen van deze regels gaan we na of we minstens 100 waarnemingen van de soort hebben.
Hierbij is een waarneming de telling van minstens één dier in een combinatie van jaar, periode, hok en meetpunt.
Stel dat we een soort hebben waarbij we de drie periodes weerhouden en de waarnemer ziet de soort op elk van de zes meetpunten tijdens alle periodes van een bepaald jaar, dan hebben we $3 \times 6 = 18$ waarnemingen voor dat hok in dat jaar.
Als we minder dan 100 waarnemingen voorhouden voor een soort, dan stopt voorlopig het verhaal voor deze soort.
We hebben dan te weinig waarnemingen voor een betrouwbare trendberekening.
Dit wil niet zeggen dat we deze soort permanent afschrijven!
Wanneer we deze regels in de toekomst opnieuw toepassen hebben we mogelijk wel voldoende waarnemingen.
Ook hokken of strata in de toekomst gekoloniseerd worden, worden dan opgepikt.
Trendberekening {#s:trendberekening}
Niet-lineaire trends
Om de trends te berekenen, voorspellen we de waargenomen aantallen aan de hand van een statistisch model.
$Y_{jpshm}$ is de voorspelling voor jaar $j$, periode $p$, stratum $s$, hok $h$ en meetpunt $m$.
Veronderstellen dat deze aantallen uit een zero-inflated negative binomial^[Een negatief binomiale verdeling met een overmaat aan nullen.] verdeling komen met parameters $\mu_{jpshm}$, $\nu$ en $\pi$.
$$Y_{jpshm} \sim ZINB(\mu_{jpshm}, \nu, \pi)$$
Laat ons eerst focussen op $\mu_{jpshm}$.
Deze parameter is via een $\log$ link gerelateerd aan de lineaire predictor $\eta_{jpshm}$.
Door deze $\log$ link kunnen we de individuele termen van de lineare predictor interpreteren als relatieve effecten.
$$\log\mu_{jpshm} = \eta_{jpshm}$$
En nu wordt het interessant want deze lineaire predictor hangt af van vijf termen:
- $\beta_s$: het effect van stratum $s$.
- $\beta_p$: het verschil tussen periode $p$ en de referentieperiode.
- $b_h$: het effect van hok $h$.
Deze effecten komen uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_h$.
- $b_m$: het effect van meetpunt $m$.
Deze effecten komen uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_m$.
- $b_{j,s}$: het effect van jaar $j$ in stratum $s$.
Dit effect modelleert een eerste orde random walk per stratum.
Het verschil tussen twee opeenvolgende jaren komt uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_j$.
$$\eta_{jpshm} = \beta_s + \beta_p + b_h + b_m + b_{j,s}$$
$$b_h \sim \mathcal{N}(0, \sigma^2_h)$$
$$b_m \sim \mathcal{N}(0, \sigma^2_m)$$
$$b_{j,s} - b_{j-1, s} = \Delta b_{j,s} \sim \mathcal{N}(0, \sigma^2_j)$$
De random walk $b_{j,s}$ is de term die ons het meest interesseert aangezien deze modelleert hoe de aantallen in de tijd wijzigen.
Het komt er op neer dat de aantallen in jaar $j$ en stratum $s$ een factor $e^{\Delta b_{j,s}}$ verschillen van het voorgaande jaar in datzelfde stratum.
De variantie $\sigma^2_j$ bepaalt hoe groot de sterkste schommelingen tussen de opeenvolgende jaren kunnen zijn.
Het model laat toe dat elk stratum zijn eigen trend heeft.
Statistische verdelingen
Een standaard distributie voor aantallen is de Poisson verdeling.
Deze heeft als kenmerk dat de variantie gelijk is aan het gemiddelde ($\sigma^2 = \mu$).
Bij veel ecologisch gegevens zien we dat de variantie in de praktijk groter is dan het gemiddelde, een fenomeen dat we overdispersie noemen.
In dat geval kunnen we een negatief binomiale verdeling gebruiken.
Deze heeft als variantie $\sigma^2 = \mu + \mu^2/\nu$, waarbij de parameter $\nu$ is een maat voor overdispersie.
Merk op dat als de overdispersieparameter zeer groot wordt ($\nu = \infty$), dan wordt de term $\mu/\nu = 0$ en bijgevolg reduceert de negatief binomiale verdeling dan tot een Poisson verdeling.
Sommige ecologische gegevens vertonen een "overmaat" aan nullen, dat zijn meer nulwaarnemingen dan de Poisson of negatief binomiale verdeling kan modelleren.
In dergelijke gevallen kunnen we overschakelen naar de zero-inflated versie van deze verdelingen.
Deze hebben een parameter $\pi$ die een maat is voor de kans op een overmatige nul.
Bij het modelleren zullen we in eerste instantie deze twee parameters instellen op $\nu = \infty$ en $\pi = 0$, m.a.w. geen overdispersie en geen overmaat aan nullen zodat we een Poisson verdeling krijgen.
Vervolgens gaan we na of er voldoende aanwijzingen zijn voor overdispersie of een overmaat aan nullen.
In het geval van overdispersie laten we $\nu$ door het model schatten waardoor we overgaan van een Poisson naar een negatief binomiaal.
In het geval van een overmaat aan nullen laten we $\pi$ door het model schatten waardoor we overgaan van een Poisson naar een zero-inflated Poisson.
Soms hebben we zowel aanwijzingen voor overdispersie als een overmaat aan nullen.
In die gevallen kiezen we de negatief binomiaal of zero-inflated Poisson naargelang welke de sterkste aanwijzingen heeft.
Vervolgens gaan we na of we de andere parameter ook nog een probleem vormt.
Zo ja, gaan we over naar een zero-inflated negatief binomiaal waarbij het model zowel $\nu$ als $\pi$ zal schatten.
Lineaire trends {#s:lineairetrend}
Lineaire trends veronderstellen dat er een constante wijziging is over de volledige looptijd.
Het model dat we hiervoor gebruiken is nagenoeg identiek aan het niet-lineaire model.
Het enige verschil zit in de lineaire predictor waar we de first order random walk ($b_{j,s}$) vervangen door een lineaire trend per stratum ($\beta_{s1} j$).
$$\eta_{jpshm} = \beta_s + \beta_p + b_h + b_m + \beta_{s1} j$$
Het lineaire model is een vereenvoudiging van het niet-lineaire model.
We kunnen deze modellen met elkaar vergelijken op basis van het Wantanabe-Akaike Information Criterion (WAIC) [@gelmanUnderstandingPredictiveInformation2014].
De WAIC waarde daalt naarmate het model de gegevens beter kan beschrijven en stijgt wanneer het model complexer wordt.
Als we modellen met elkaar vergelijken op basis van WAIC, zal het model met de laagste WAIC de beste mix zijn tussen een goede beschrijving van de gegevens en een zo eenvoudig mogelijk model.
We beschouwen de trend als lineair wanneer het lineaire model de laagste WAIC heeft.
Wanneer de WAIC van het niet-lineaire model minder dan 2 eenheden lager is dan deze van het lineaire model, beschouwen de trend als mogelijk niet-lineair.
Pas wanneer de WAIC van het niet-lineaire model duidelijk lager is dan het lineaire model, stellen we dat de trend niet-lineair is.
In dat geval de gebruiker moet de gerapporteerde lineaire trend met de nodige voorzichtigheid interpreteren.
Kijk zeker naar de bijhorende niet-lineaire trend vooraleer de cijfers te interpreteren.
Resultaten per driejarige cyclus
Zoals aangegeven in §\@ref(s:steekproef) zullen we een bepaald hok in principe om de drie jaar herbezoeken.
Hierdoor krijgen we drie sets van hokken.
Er bestaat een kans dat de schatting voor een bepaald jaar beïnvloed is door de set van hokken die in dat jaar onderzocht worden.
Om dit effect uit te schakelen, analyseren we de gegevens tevens op basis van de driejarige cyclus $c$ i.p.v. jaar $j$.
De eerste cyclus omvat de eerste drie jaar sinds de start van het meetnet (2007-2009).
De volgende cyclus omvat telkens de volgende drie jaar aansluitend op de vorige cyclus.
Aangezien we alle beschikbare gegevens gebruiken bij de analyse bevat de laatste cyclus mogelijk minder dan drie jaar.
Dit is duidelijk zichtbaar doordat het laatste jaar van de cyclus op dat ogenblik in de toekomst ligt.
Als bijvoorbeeld 2019 het meest recente jaar met gegevens is, dan is de laatste cyclus 2019-2021.
Hieronder geven we aan op welke manier we de eerste beschreven modellen aanpassen.
$$Y_{cpshm} \sim ZINB(\mu_{cpshm}, \nu, \pi)$$
$$\log\mu_{cpshm} = \eta_{cpshm}$$
Niet-lineaire trend
$$\eta_{cpshm} = \beta_s + \beta_p + b_h + b_m + b_{c,s}$$
$$b_{c,s} - b_{c-1, s} = \Delta b_{c,s} \sim \mathcal{N}(0, \sigma^2_c)$$
$b_{c,s}$: het effect van driejarige cyclus $c$ in stratum $s$.
Dit effect modelleert een eerste orde random walk per stratum.
Het verschil tussen twee opeenvolgende driejarige cycli komt uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_c$.
Lineaire trend
$$\eta_{cpshm} = \beta_s + \beta_p + b_h + b_m + \beta_{s1} c$$
Modellen fitten
We fitten de statische modellen in R [@R] met het INLA package [@rueINLAFunctionsWhich2009].
INLA gebruikt een Bayesiaanse benadering om de modellen te fitten.
Daarom moeten we priors specificeren voor de parameters en hyperparameters.
- $\beta_s$, $\beta_p$ en $\beta_{s1}$ krijgen een Gaussiaanse prior $\mathcal{N}(0, 5)$.
- $\sigma^2_h$ krijgt een
penalised complexity prior [@simpson.rue.ea2017PenalisingModelComponent] zodat $Prob(\sigma_h > 1) = 0.01$.
- $\sigma^2_m$ krijgt een
penalised complexity prior zodat $Prob(\sigma_m > 1) = 0.01$.
- $\sigma^2_j$ krijgt een
penalised complexity prior zodat $Prob(\sigma_j > 0.5) = 0.01$.
- $\sigma^2_c$ krijgt een
penalised complexity prior zodat $Prob(\sigma_c > 0.5) = 0.01$.
- $\nu$ krijgt een Gamma prior zodat $\nu \sim \Gamma(1/7, 1/7)$.
- $\pi$ wordt intern geparametriseerd als $\theta = \log\frac{\pi}{1 - \pi}$ met Gaussiaanse prior $\theta ~ \sim \mathcal{N}(-1, 5)$.
Gemiddeld aantal dieren per meetpunt
Een gemiddeld hok heeft als effect $b_h = 0$ en een gemiddeld meetpunt $b_m = 0$.
In de referentieperiode is $\beta_p = 0$.
In deze gevallen vereenvoudigt de lineaire predictor tot het effect van stratum $s$ en de trend in dat stratum:
$$\eta_{js} = \beta_s + b_{j,s}$$
Door de effecten van de strata te vermenigvuldigen met hun stratumgewicht ($\gamma_s$) krijgen we een schatting voor Vlaanderen.
$$\eta_{j} = \sum_s(\gamma_s\beta_s + \gamma_sb_{j,s})$$
Om de schatting van het gemiddelde aantal in Vlaanderen te krijgen, moeten we de lineaire predictor terug omzetten van de $\log$ schaal naar de natuurlijke schaal en corrigeren voor de eventuele overmaat aan nullen.
We krijgen dan:
$$E[Y_j] = (1 - \pi) e^{\eta_j} \prod_s(e^{\gamma_s \beta_s} e ^ {\gamma_s b_{j,s}})$$
Vergelijken van jaren
We kunnen twee jaren $a$ en $b$ met elkaar vergelijken door hun lineaire predictoren voor het gemiddeld aantal dieren van elkaar af te trekken.
$$\eta_a - \eta_b = \sum_s(\gamma_s\beta_s + \gamma_sb_{a,s}) - \sum_s(\gamma_s\beta_s + \gamma_sb_{b,s})$$
Aangezien het globale stratumeffecten en de stratumgewichten niet wijzigen in de tijd, kunnen we dit vereenvoudigen tot
$$\eta_a - \eta_b = \sum_s\gamma_s(b_{a,s} -b_{b,s})$$
De wijziging in Vlaanderen is het gewogen gemiddelde van de wijzigingen in de strata.
Na omzetting van de $log$ schaal naar de natuurlijke schaal krijgen we het relatieve aantal $I_{a|b}$ van jaar $a$ t.o.v. jaar $b$
$$I_{a|b}=\frac{e^{\nu_a}}{e^{\nu_b}} = \prod_s\left(\frac{e^{\gamma_sb_{a,s}}}{e^{\gamma_sb_{b,s}}}\right)$$
Wanneer we een jaar met zichzelf vergelijken krijgen we per definitie $I_{a|a} = 1 = 100\%$.
Wanneer we meerdere jaren met eenzelfde referentiejaar vergelijken krijgen we een indexwaarde: het relatieve verschil van elk jaar t.o.v. van een bepaald referentiejaar.
Gewicht van de strata {#s:stratumgewicht}
In het vorige onderdeel hebben gebruikt gemaakt van stratumgewichten $\gamma_s$ zonder deze te definiëren.
We bepalen het gewicht van een stratum op basis van drie kenmerken: het totaal aantal hokken van het stratum in Vlaanderen ($N_s$), het aantal onderzochte hokken in het stratum ($T_s$) ongeacht of ze al dan niet relevant zijn voor de soort en het aantal relevante hokken voor de soort in het stratum ($R_s$).
Het aandeel relevante hokken per stratum ($R_s/T_s$) is een goede maat van de frequentie waarmee een soort aanwezig is binnen een stratum.
Wanneer een soort in nagenoeg alle onderzochte hokken van het stratum waargenomen wordt is $R_s/T_s \simeq 1$.
Is de soort zeer zeldzaam (komt slechts in een paar van de onderzochte hokken voor) dan is $R_s/T_s \simeq 0$.
Het basisgewicht $\gamma_{sb}$ voor stratum $s$ is het aandeel relevante hokken vermenigvuldigd met het totaal aantal hokken in het stratum.
Dit is een schatting van het totaal aantal hokken in het stratum waar de soort voldoende frequent voorkomt.
$$\gamma_{sb} = \frac{R_s}{T_s}N_s$$
Om makkelijker te kunnen rekenen delen we de basisgewichten door hun som zodat de stratumgewichten sommeren tot 1.
$$\gamma_s = \frac{\gamma_{sb}}{\sum_s\gamma_{sb}}$$
Merk op dat het aantal relevante hokken per stratum ($R_s$) soortafhankelijk is.
Bijgevolg zijn de stratumgewichten eveneens soortafhankelijk.
Samengestelde indices {#s:samengesteld}
We berekenen een aantal samengestelde indices die de trends voor een groep van soorten aggregeren.
Deze indices zijn het meetkundig gemiddelde van de verschillen tussen jaren voor alle soorten van de groep.
Bij de berekening maken we gebruik van het trucje dat we een product kunnen schrijven als een som van logaritmes die we nadien terug exponentiëren.
Een meetkundig gemiddelde in de natuurlijke schaal kunnen we dat schrijven als een rekenkundig gemiddelde in de log-schaal.
$$\sqrt[n]{\prod_{i+1}^na_i} = \exp\left(\frac{\sum_{i = 1}^n\log a_i}{n}\right)$$
Het voordeel aan deze techniek is dat de schattingen van de paarsgewijze verschillen tussen de jaren reeds beschikbaar zijn in de log-schaal.
Bovendien beschikken we tevens over hun variantie in de log-schaal.
Dit laat ons toe om makkelijk het betrouwbaarheidsinterval te berekenen aan de hand van een paar vuistregels.
- Bij onafhankelijke variabelen geldt dat de variantie van hun som gelijk is aan de som van de varianties. $\sigma^2_{\sum X_i} = \sum \sigma^2_{X_i}$
- De variantie van het product van een variabele en een constante is gelijk aan de variantie van de variabele vermenigvuldigd met die constante. $\sigma^2_{aX} = a \sigma^2_X$
- Het model veronderstelt dat de individuele modelparameters een Gaussiaanse verdeling volgen.
Hierdoor volgen de paarsgewijze verschillen tussen de jaren en hun gemiddelde (telkens in de log-schaal) een Gaussiaanse verdeling.
Aangezien we zowel het gemiddeld als de variantie van deze verdeling kennen, kunnen we hieruit de gewenste kwantielen voor de betrouwbaarheidsintervallen berekenen (nog steeds in de log-schaal).
- De log-transformatie is een monotoon stijgende functie aangezien $\log(x)$ steeds groter wordt wanneer $x$ groter wordt.
Een kenmerk van een monotoon stijgende transformatie is dat ze de volgorde behouden: het 5% kleinste element zal ook na de transformatie het 5% kleinste element zijn.
Waardoor we de kwantielen van de betrouwbaarheidsintervallen zonder probleem kunnen terugrekenen naar de natuurlijke schaal.
Voorstelling van de gegevens
Onzekerheid {#s:onzekerheid}
Alle resultaten zijn gebaseerd op een steekproef en op de waarnemingen zit onvermijdelijk een zekere meetfout.
Vandaar dat we bij de puntschattingen tevens een betrouwbaarheidsinterval weergeven.
In de tekst gebruiken we het 90% (5%; 95%) interval waarbij er 5% kans is dat de werkelijke waarde kleiner is dan de ondergrens en 5% dat ze groter is dan de bovengrens.
Het 90% interval is iets smaller dan het traditionele 95% (2.5%; 97.5%) interval.
Door een smaller interval te kiezen zullen we sneller uitspraken kunnen doen, waardoor de kans kleiner wordt dat we ten onrechte stellen dat er geen effect is.
De prijs die we hiervoor betalen is dat de kans dat we ten onrechte stellen dat er een significant effect is, stijgt van 5% naar 10%.
De doelstelling van de algemene broedvogelmonitoring is zo spoedig mogelijk detecteren wanneer er iets aan de hand is de broedvogels.
Vanuit dat oogpunt is het vermijden van vals negatieve signalen (ten onrechte stellen dat er niets aan de hand is) belangrijker dan het vermijden van vals positieve signalen (ten onrechte stellen dat er iets aan de hand is).
Wanneer we het interval op een figuur (bijvoorbeeld fig. \@ref(fig:trendklasse)) weergeven, vullen we het 90% interval aan met een 60% (20%; 80%) interval en een 30% (35%; 65%) interval.
@britton.fisher.ea1998InflationReportProjections waren de inspiratie voor deze manier van weergeven.
Het 30% interval vormt het donkerste deel van het interval en geeft de meeste waarschijnlijke locatie van de werkelijke waarde weer.
Naarmate het interval lichter wordt, neemt de onzekerheid toe.
Opdeling van de effecten in een aantal klassen {#s:trendklasse}
Bij het niet-lineaire model berekenen we alle paarsgewijze relatieve verschillen tussen de verschillende jaren.
Bij het lineaire model hebben we de gemiddelde jaarlijkse relatieve verandering.
Deze laatste rekenen we tevens om naar de totale wijziging over de looptijd van het meetnet omdat dit eenvoudiger te interpreteren is.
Vergelijk een daling met -5% per jaar of een daling met -50% over 15 jaar.
Deze laatste klinkt dramatischer door het grotere cijfer, terwijl -5% per jaar overeenkomt met -54% over 15 jaar.
De totale wijziging over de looptijd van het lineaire model is tevens vergelijkbaar met de wijziging tussen het eerste en laatste jaar van het niet-lineaire model.
Om de interpretatie makkelijker te maken, delen we de wijzigingen op in tien klassen door hun 90% interval te vergelijken met een referentie, ondergrens en bovengrens.
Een effect is significant wanneer de referentie buiten het 90% interval ligt.
We spreken over een toename (afname) als het interval volledig boven (onder) de referentie ligt.
Niet significante effect is ook informatief wanneer het bijhorende interval voldoende smal is.
Bijvoorbeeld als het interval volledig tussen de ondergrens en bovengrens ligt.
In dat geval kunnen we stellen dat het effect niet significant en klein is, het immers zeker minder sterk dan de ondergrens en minder sterk dan de bovengrens.
Dergelijk effect krijgt de naam stabiel.
Heeft het effect een breed interval dat zowel de boven- als ondergrens bevat, spreken we over een onduidelijk effect.
Daarnaast is er nog de mogelijkheid dat het interval zowel de bovengrens (ondergrens) als de referentie bevat maar niet de ondergrens (bovengrens).
Dan spreken we over een mogelijke toename (mogelijke afname).
We kunnen de boven- en ondergrens eveneens gebruiken om een verder onderscheid te maken binnen de significante effecten.
Een interval volledig boven (onder) de bovengrens (ondergrens) wordt dan een sterke toename (sterke afname).
Een interval volledig tussen de referentie en de de bovengrens (ondergrens) wordt dan een matige toename (matige afname).
Een interval dat de referentie niet bevat maar wel de bovengrens (ondergrens) blijft een toename (afname).
Merk op dat de indeling volledig gebaseerd is op de onzekerheid rond het effect en niet op de puntschatting van het effect zelf.
We vatten de opdeling met bijhorende afkortingen en regels samen in tabel \@ref(tab:regels).
Figuur \@ref(fig:trendklasse) geeft een grafische voorstelling waarbij we de afkortingen in combinatie met aangepaste symbolen gebruiken.
De afkortingen zelf zijn te fijn om als symbool te gebruiken.
(ref:regels) Overzicht van de benamingen van de tien effectklassen met hun afkorting en de regels. $R$: referentie, $L$: ondergrens, $B$: bovengrens, $l$: ondergrens van het 90% interval, $b$: bovengrens van het 90% interval. $L < R < B$ en $l < b$.
tribble(
~ benaming, ~afkorting, ~regels,
"sterke toename", "`++`", "$B < l$",
"toename", "`+`", "$R < l < B$ en $B < b$",
"matige toename", "`+~`", "$R < l < B$ en $b < B$",
"stabiel", "`~`", "$L < l < R$ en $R < b < B$",
"matige afname", "`-~`", "$L < l < R$ en $b < R$",
"afname", "`-`", "$l < L$ en $L < b < R$",
"sterke afname", "`--`", "$l < L$",
"mogelijke toename", "`?+`", "$L < l < R$ en $B < b$",
"mogelijke afname", "`?-`", "$l < L$ en $R < b < B$",
"onduidelijk", "`?`", "$l < L$ en $B < b$"
) %>%
kable(
caption = "(ref:regels)", format = "pandoc", escape = FALSE, align = "lcc"
)
(ref:trendklasse) Voorbeeld van de tien mogelijke interpretaties van een effect door het 90% interval te vergelijken met een referentie, ondergrens en bovengrens.
ggplot(effect, aes(x = x_int, y = y, ymin = lcl, ymax = ucl, link_sd = s)) +
geom_hline(yintercept = 0, linetype = 2) +
geom_hline(yintercept = c(-1, 1), linetype = 3) +
geom_rect(
aes(xmin = x_int - 0.4, xmax = x_int + 0.4, fill = klasse),
show.legend = FALSE
) +
geom_point(size = 6, colour = inbo_hoofd) +
geom_text(
aes(label = klasse), hjust = 0.5, vjust = 0.5, colour = "white"
) +
scale_x_continuous(breaks = effect$x_int, labels = effect$x) +
scale_y_continuous(
"effect", breaks = c(-1, 0, 1),
labels = c("ondergrens", "referentie", "bovengrens")
) +
scale_fill_manual(values = kleurgradient) +
coord_flip() +
theme(axis.title.y = element_blank())
Uiteraard hangt de opdeling sterk af van de keuze van de boven- en ondergrens.
De soortenmeetnetten voor de Natura 2000 monitoring streven er naar om een daling in populatiegrootte met -25% over 25 jaar tijd vast te kunnen stellen.
Hierbij wordt een wijziging in de populatiegrootte van -25% als belangrijk aanzien.
Daarom hanteren we voor de algemene broedvogelmonitoring ook -25% als ondergrens, zowel bij de paarsgewijze verschillen tussen de jaren als de lineaire wijziging over de looptijd van het meetnet.
Een van daling van -25% komt overeen met aantallen die nog drie kwart van de uitgangssituatie bedragen.
Een even sterke wijziging in de omgekeerde richting zorgt er voor dat de aantallen toenemen tot vier derde van de uitgangssituatie, of een toename met +33% wat we als bovengrens gebruiken.
Bij de samengestelde indices hebben we aangepaste grenswaarden nodig.
Veronderstel een samengestelde index op basis van $n$ soorten.
We berekenen het rekenkundig gemiddelde in de log-schaal, dan is de variantie van dit gemiddelde een factor $n$ kleiner dan de som van de varianties.
De breedte van een betrouwbaarheidsinterval hangt samen met de standaard afwijking, wat de vierkantswortel van de variantie is.
Hierdoor zullen de breedtes van de betrouwbaarheidsintervallen van de samengestelde index een factor $\sqrt{n}$ kleiner zijn.
Vandaar dat we de grenswaarden van de indices tevens aanpassen door ze met een factor $\sqrt{n}$ te verkleinen.
Elke samengestelde index heeft zijn eigen soortenlijst met een variabele aantal soorten.
Daarom zullen we bij elke samengestelde index zijn aangepaste grenswaarden vermelden.
Overzicht van de lineaire trends
Tabellen \@ref(tab:driejaarlijkse-trends) en \@ref(tab:jaarlijkse-trends) geven een overzicht van de lineaire trends voor elke soort.
We hebben de trends gesorteerd volgens opdeling van de klassen en binnen de klasse volgens puntschatting van de trend.
Hierdoor start de tabel met de soorten met de sterkste positieve trends.
Bij de volgende soorten zal de trend minder sterk worden tot we aan de soorten met een stabiele trend komen.
Daarna volgen de soorten met een negatieve trend waarbij de trend steeds sterker negatief wordt.
We sluiten de tabel af met de soorten met een mogelijke of onduidelijke trend.
De tabel bevat de trend uitgedrukt als een jaarlijkse wijziging en als een wijziging over de volledige looptijd van het meetnet.
Verder bevat de tabel de opdeling van de trend in klassen en een indicatie of de trend al dan niet lineair is (zie §\@ref(s:lineairetrend)).
De naam van de soort is een snelkoppeling naar de detail van de soort zelf.
We raden aan om hiervan gebruikt te maken bij de interpretatie van niet-lineaire trends.
Evolutie van gemiddelde aantallen per soort
Deze figuur geeft de evolutie van de gemiddelde aantallen per meetpunt weer volgens het niet-lineair model.
Indien het model lineair is, zal het patroon van de niet-lineaire trend dicht bij een lineaire trend liggen.
De lijn bevat de puntschatting van het gemiddelde aantal in elk jaar.
Rond de lijn tonen we het 30%, 60% en 90% interval (zie §\@ref(s:onzekerheid)).
In de online versie van het rapport is de figuur interactief.
Wanneer de gebruiker met de muis over de figuur gaat verschijnt er in de buurt van de lijn een pop-up met de exacte schatting van dat jaar inclusief het 90% betrouwbaarheidsinterval en het jaartal.
Paarsgewijze vergelijking van jaren
Deze informatie hebben we telkens in twee figuren samengevat.
De eerste figuur geeft het verschil van elk jaar t.o.v. 2007, het jaar waarin de metingen gestart zijn.
Voor 2007 geven we geen cijfer omdat het per definitie 1 is en geen informatie bevat.
Het symbool geeft de puntschatting van het relatieve verschil t.o.v. 2007 weer.
De vorm van het symbool geeft de opdeling van de sterkte van het effect weer (zie fig. \@ref(fig:trendklasse)).
Rond elke puntschatting tonen we het 30%, 60% en 90% interval (zie §\@ref(s:onzekerheid)).
De horizontale streepjeslijn geeft de referentie van 0% verschil weer.
De horizontale puntlijnen geven de bovengrens (+33%) en ondergrens (-25%) weer.
Deze laten toe om de vlot zelf de betrouwbaarheidsintervallen te vergelijken met de referentie, boven- en ondergrens.
De figuur heeft twee y-assen.
Beide assen geven dezelfde informatie weer, enkel de formattering van de labels is anders.
De linkeras toont procentuele verschillen terwijl de rechteras de relatieve verschillen toont.
Een procentueel verschil van +50% is equivalent met een relatieve verschil (verhouding) met een factor 1,5.
We stellen vast dat veel gebruikers ook andere jaren met elkaar willen vergelijken.
Een correcte vergelijking is enkel mogelijk indien we een van deze jaren als referentie gebruiken.
Voor elk jaar een afzonderlijke figuur maken, zou het rapport onoverzichtelijk groot maken.
Om de vergelijkingen toch mogelijk te maken hebben we alle paarsgewijze verschillen tussen de jaren in een raster weergeven.
Elke rij in het raster staat voor een ander referentiejaar.
De kolommen geven de verschillen van een bepaald jaar weer t.o.v. de verschillende referentiejaren.
Op de diagonaal staan geen waarden omdat we daar een jaar met zichzelf vergelijken.
De kleur van de symbolen geeft de sterkte van het verschil (gebaseerd op de puntschatting).
Zwakke verschillen zijn grijs, sterke positieve verschillen rood, sterke negatieve verschillen blauw.
De vorm van de symbolen geeft zicht op de sterkte en onzekerheid van het effect (zie tab. \@ref(tab:regels)).
Wanneer een rij volledig rood (blauw) is, zijn alle verschillen met dit referentiejaar positief (negatief) m.a.w. dit is het referentiejaar met de laagste (hoogste) aantallen.
Wanneer een kolom volledig rood (blauw) is, zijn alle verschillen van dit jaar t.o.v. alle referentiejaren negatief (positief) m.a.w. dit is het jaar met de hoogste (laagste) aantallen.
Clusters van punten met een gelijkaardige kleur geven periodes aan waarin de aantallen geleidelijk wijzigen.
In de online versie van dit rapport is deze figuur interactief.
Wanneer de gebruiker met de muis over de figuur gaat verschijnt er in de buurt van de lijn een pop-up met de exacte schatting van dat punt inclusief het 90% betrouwbaarheidsinterval en het jaartal en referentiejaar.
Reproduceerbaarheid en traceerbaarheid van de gegevensverwerking
Versiebeheer
De waarnemers voeren hun waarnemingen via een webapplicatie toe aan de databank.
De analyse start met het importeren van de relevante gegevens uit de databank aan de hand van een R-script.
De bekomen dataset bewaren we als een collectie van tekstbestanden met behulp van het git2rdata package [@git2rdata].
Zowel deze code als de tekstbestanden bewaren we onze versiebeheer met git.
Git is een gratis en open source gedistribueerd versiebeheersysteem dat is ontworpen om alles van kleine tot zeer grote projecten snel en efficiëntie te beheren.
De code zit vervat in het abvanalysis package [@abvanalysis], dat vrij online beschikbaar is.
De tekstbestanden met de gegevens zijn momenteel niet publiek toegankelijk.
De ruwe gegevens worden geaggregreerd per hok en met drie jaar vertraging publiek ontsloten via GBIF [@vermeersch.anselin.ea2018ABVCommonBreeding].
Reproduceerbaarheid
Om de reproduceerbaarheid te garanderen starten we de analyses steeds vanaf de tekstbestanden.
Aan de hand van code in het abvanalysis package definiëren we hoe we de tekstbestanden omzetten naar de verschillende analyses (§\@ref(s:trendberekening)).
Dit resulteert in een resem op zichzelf staande analyseobjecten van het n2kanalysis package [@n2kanalysis].
Dergelijk analyseobject bevat naast de nodige gegevens en de modeldefinitie tevens de nodige metadata zoals alle gebruikte R packages, inclusief hun versie, en een verwijzing naar de voorafgaande analyses waarvan deze analyse afhangt.
Zo verwijst de analyse van een samengestelde index (§\@ref(s:samengesteld)) naar de trendberekening van de bijhorende soorten (§\@ref(s:trendberekening)).
De individuele trendberekeningen verwijzen naar de "analyse" die de gegevens importeerde uit de databank.
Deze laatste bevat de nodige links naar de tekstbestanden en hun versie.
Traceerbaarheid {#s:traceerbaarheid}
De n2kanalysis objecten hebben elk twee data hashes.
Elk data hash is een reeks van 40 hexadecimale^[Hexadecimaal betekent letterlijk zestientallig. Het is een talstelsel waarbij niet, zoals gebruikelijk, met tien cijfers wordt gewerkt, maar met zestien cijfers. De cijfers 0 t.e.m. 9 worden daarom uitgebreid met a (=10) t.e.m. f (=15).] cijfers die het resultaat zijn van de cryptografische hashfunctie SHA-1.
Deze hashfunctie heeft een aantal belangrijke eigenschappen:
- Ze zetten elke invoer om naar een uitvoer met vaste lengte (40 hexadecimale cijfers).
- De uitvoer is stabiel: als je de hash van een bepaalde invoer opnieuw berekent krijg je steeds dezelfde uitvoer.
- Het is niet mogelijk om de invoer te reconstrueren op basis van de uitvoer.
- Eender welke kleine wijziging aan de invoer resulteert in een sterke wijziging van de uitvoer.
- De kans dat twee verschillende invoeren dezelfde uitvoer opleveren is zeer klein.
De eerste data hash van het analyseobject is gebaseerd op alle informatie die gekend is op het moment dat we het analyseobject definiëren en die nooit zal wijzigen tijdens de analyse.
Denk hierbij aan de definitie van de analyse, de gegevens, de soort, ...
Gezien de eigenschappen van de hashfunctie kunnen we deze data hash gebruiken om ondubbelzinnig te verwijzen naar een specifieke analyse (inclusief de gebruikte gegevens.)
De tweede data hash van het analyseobject baseren we enerzijds op de eerste data hash en anderzijds van alle onderdelen van het analyseobject die wijzigen in de loop van de analyse.
Dit is o.a. het resultaat van de statistische analyse, de gebruikte software, ...
De analyse opnieuw uitrekenen met software van een andere versie zal de tweede data hash aanpassen.
De combinatie van deze twee data hashes laat enerzijds toe om naar een specifieke versie van de analyse te verwijzen.
Anderzijds bieden ze een garantie over de inhoud van het analyseobject.
In het geval van twijfel over een analyse kunnen we teruggrijpen naar de analyseobject in kwestie.
De data hashes bewijzen dan dat we het correcte analyseobject hebben.
inbo/abvanalysis documentation built on April 6, 2023, 5:14 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
(PART) Methodologie {-}
Inzameling van de tellingen in het veld
Hier beschrijven we beknopt het veldwerk, voor de details verwijzen we naar @ABVmethode.
Steekproef {#s:steekproef}
De UTM 1 x 1 km hokken vormen de basis waaruit de steekproef getrokken is. Bij de start van het meetnet hebben we de hokken opgedeeld in een aantal strata. De regels zijn gebaseerd op het oppervlakteaandeel van een bepaald landgebruik op basis van de Biologische Waarderingkaart [@BWK].
- Landbouw: minstens 80% landbouw.
r strata$n[strata$description == "Landbouw"]hokken. - Urbaan: minstens 80% urbaan.
r strata$n[strata$description == "Urbaan"]hokken. - Bos: minstens 80% bos.
r strata$n[strata$description == "Bos"]hokken. - Suburbaan: minstens 80% suburbaan.
r strata$n[strata$description == "Suburbaan"]hokken. - Heide en duin: minstens 20% heide of duin.
r strata$n[strata$description == "Heide en duin"]hokken. - Moeras en water: minstens 20% moeras en water.
r strata$n[strata$description == "Moeras en water"]hokken.
Uit deze set trekken we een aselecte, gestratificeerde steekproef van 1200 hokken waarbij zeldzamere habitats overbemonsterd worden. We streven er naar om jaarlijks 300 hokken te bemonsteren in een driejarige rotatie. De waarnemers mochten in het eerste jaar 300 hokken kiezen uit de set van 1200. Deze set van hokken komen in principe opnieuw aan bod in jaren 4, 7, 10, ... In jaar 2 kiezen ze 300 hokken uit de overgebleven 900 hokken. Deze set komt opnieuw aan bod in jaren 5, 8, 11, ... Tenslotte kiezen de waarnemers in het derde jaar een laatste set van 300 hokken uit de laatste 900 hokken. Deze set hokken bemonsteren we in de jaren 3, 6, 9, 12, ...
Figuur \@ref(fig:inspanning) geeft de effectieve monitoringsinspanning weer. In deze figuur hebben we de hokken gesorteerd volgens 1) het eerste jaar met gegevens, 2) het laatste jaar met gegevens, 3) het tweede jaar met gegevens, 4) het derde jaar met gegevens, ... Hierdoor staan hokken met een meer gelijkende onderzoeksgeschiedenis dicht bij elkaar. Merk op dat de driejarige cyclus voor de meest hokken wordt gerespecteerd. Voor sommige hokken is de inspanning variabel, soms frequenter dan om de drie jaar, soms zit er meer tijd tussen. Sommige hokken werden slechts in een of twee jaar onderzocht (fig. \@ref(fig:inspanning-hist)). Voor het onderscheid tussen mogelijk bruikbaar en voorlopig niet bruikbaar verwijzen we naar §\@ref(s:relevant).
effort %>% arrange(.data$year) %>% group_by(.data$location) %>% transmute( .data$year, end = max(.data$year), year2 = lead(.data$year), year3 = lead(.data$year2), year4 = lead(.data$year3), year5 = lead(.data$year4), year6 = lead(.data$year5), year7 = lead(.data$year6) ) %>% slice(1) %>% ungroup() %>% arrange( .data$year, .data$end, .data$year2, .data$year3, .data$year4, .data$year5, .data$year6, .data$year7 ) %>% rowid_to_column("hok") %>% transmute( .data$location, hok = factor(.data$hok), bruikbaar = ((.data$year - 2007) %/% 3) < ((end - 2007) %/% 3), bruikbaar = factor( .data$bruikbaar, levels = c(TRUE, FALSE), labels = c("mogelijk", "voorlopig\nniet") ) ) %>% inner_join(effort, by = "location") -> inspanning ggplot(inspanning, aes(x = year, y = hok, group = hok, colour = bruikbaar)) + geom_line(alpha = 0.2) + geom_point(alpha = 0.5, size = ifelse(output_format == "html", 4, 1)) + facet_wrap(~stratum, scales = "free_y") + scale_x_continuous(breaks = pretty) + theme( axis.title.x = element_blank(), axis.text.x = element_text(angle = 90), axis.text.y = element_blank() )
inspanning %>% count(stratum, location, bruikbaar) %>% ggplot(aes(x = n, fill = bruikbaar)) + geom_histogram(binwidth = 1) + facet_wrap(~stratum) + scale_x_continuous( "Aantal jaren waarin een hok onderzocht is", breaks = pretty, limits = c(0, NA) ) + ylab("Aantal hokken")
Steekproefeenheid
De steekproefeenheid bestaat uit een hok van 1 x 1 km. Binnen dit hok worden zes punten vastgelegd in een vaste patroon (fig. \@ref(fig:patroon)). Indien de punten in de praktijk niet bereikbaar zijn, mag de waarnemer ze verplaatsen naar het dichtstbijzijnde bereikbare punt. De waarnemer documenteert deze wijziging zodat we in de toekomst steeds op dezelfde punten blijven waarnemen.
expand.grid( x = seq(1 / 4, by = 1 / 2, length = 2), y = seq(1 / 6, by = 1 / 3, length = 3) ) %>% ggplot(aes(x = x, y = y)) + geom_rect( data = tribble(~x, ~y, ~xend, ~yend, 0, 0, 1, 1), aes(xmin = x, ymin = y, xmax = xend, ymax = yend), fill = NA, colour = inbo_steun_blauw ) + geom_point() + coord_fixed() + theme_void()
Tellingen
In het jaar dat we een hok bemonsteren zal de waarnemer het hok in drie periodes bezoeken: 1 maart - 15 april, 16 april - 31 mei, 1 juni - 15 juli. Tussen twee opeenvolgende bezoeken moet er minstens twee weken liggen. Alle meetpunten van een hok worden op dezelfde dag onderzocht tussen zonsopgang en 4 uur na zonsopgang. Op elk meetpunt telt de waarnemer gedurende 5 minuten het aantal volwassen vogels per soort. Overvliegende groepen vogels worden hierbij niet meegeteld.
Verwerking van de gegevens
Selectie van relevante gegevens per soort {#s:relevant}
Hoewel we alle waargenomen vogelsoorten op een meetpunt noteren, toch zijn niet alle hokken even relevant voor een bepaalde soort. Om de redenering op te bouwen, beschouwen we een fictieve soort die enkel op een afgelegen eiland voorkomt. Alle vogeltellingen buiten dit eiland zijn uiteraard niet informatief voor wijzigingen in de aantallen van deze fictieve soort. Buiten het eiland zullen de getelde aantallen immers per definitie nul zijn. Wanneer we deze tellingen toch zouden in rekening brengen, dan zullen ze de globale trend afvlakken. De globale trend is min of meer equivalent met een "gemiddelde" trend over de verschillende locaties. Als een locatie steeds nul geeft, dan zijn lokaal de aantallen wiskundig "stabiel" (geen lokale trend). Daarom is het noodzakelijk om de locaties die niet relevant zijn voor een soort buiten beschouwing te houden.
In Vlaanderen is de situatie niet zo zwart-wit als bij het voorbeeld met het afgelegen eiland. Daarom hebben nood aan een set regels die we toepassen op de beschikbare gegevens.
- We beschouwen een hok als relevant voor een soort wanneer we deze soort tijdens minstens twee verschillende driejarige cycli waarnemen. Op deze manier sluiten we hokken uit waar de soort nooit of slechts toevallig waargenomen is. Hokken die slechts binnen één driejarige cyclus onderzocht zijn, kunnen we hierdoor voorlopig niet gebruiken. Figuur \@ref(fig:potentieel) geeft een overzicht van het aantal hokken die gedurende minstens twee cycli onderzocht zijn.
- Binnen een hok is een meetpunt pas relevant wanneer we de soort er minstens eenmaal waarnemen.
- Een stratum is relevant voor een soort wanneer het minstens drie relevante hokken bevat.
- Op basis van de overblijvende gegevens berekenen we het (meetkundig[^geometrisch]) gemiddelde van de aantallen in elk van de drie perioden. De periode met het hoogste gemiddelde gebruiken we als referentieperiode. We houden enkel de periodes waarbij het gemiddelde minstens 15% van de referentieperiode bedraagt. Op deze manier sluiten we periodes uit waarbij de soort zelden gezien wordt. Bijvoorbeeld een trekvogel die pas in de loop van de tweede periode toekomt en waarvoor de eerste periode bijgevolg niet relevant is.
- We passen de regels voor relevante hokken en strata opnieuw toe. Het negeren van een of twee irrelevante periodes kan er voor zorgen dat een hok niet langer relevant is.
[^geometrisch]: Het meetkundig gemiddelde of geometrisch gemiddelde van $n$ getallen wordt verkregen door de getallen met elkaar te vermenigvuldigen en vervolgens van het product de n^de^ machtswortel te nemen. $$m=\left(\prod {i=1}^{n}a{i}\right)^{1/n}$$
inspanning %>% mutate(cyclus = (year - 2007) %/% 3) %>% distinct(location, cyclus) %>% count(location) %>% filter(n >= 2) %>% semi_join(x = inspanning, by = "location") %>% count(year, stratum, name = "hokken") %>% mutate(stratum = reorder(stratum, -hokken)) %>% ggplot(aes(x = year, y = hokken, colour = stratum)) + geom_line() + scale_x_continuous(breaks = pretty) + scale_y_continuous("Aantal potentieel relevante hokken", limits = c(0, NA)) + theme(axis.title.x = element_blank())
Na het toepassen van deze regels gaan we na of we minstens 100 waarnemingen van de soort hebben. Hierbij is een waarneming de telling van minstens één dier in een combinatie van jaar, periode, hok en meetpunt. Stel dat we een soort hebben waarbij we de drie periodes weerhouden en de waarnemer ziet de soort op elk van de zes meetpunten tijdens alle periodes van een bepaald jaar, dan hebben we $3 \times 6 = 18$ waarnemingen voor dat hok in dat jaar.
Als we minder dan 100 waarnemingen voorhouden voor een soort, dan stopt voorlopig het verhaal voor deze soort. We hebben dan te weinig waarnemingen voor een betrouwbare trendberekening. Dit wil niet zeggen dat we deze soort permanent afschrijven! Wanneer we deze regels in de toekomst opnieuw toepassen hebben we mogelijk wel voldoende waarnemingen. Ook hokken of strata in de toekomst gekoloniseerd worden, worden dan opgepikt.
Trendberekening {#s:trendberekening}
Niet-lineaire trends
Om de trends te berekenen, voorspellen we de waargenomen aantallen aan de hand van een statistisch model.
$Y_{jpshm}$ is de voorspelling voor jaar $j$, periode $p$, stratum $s$, hok $h$ en meetpunt $m$.
Veronderstellen dat deze aantallen uit een zero-inflated negative binomial^[Een negatief binomiale verdeling met een overmaat aan nullen.] verdeling komen met parameters $\mu_{jpshm}$, $\nu$ en $\pi$.
$$Y_{jpshm} \sim ZINB(\mu_{jpshm}, \nu, \pi)$$
Laat ons eerst focussen op $\mu_{jpshm}$.
Deze parameter is via een $\log$ link gerelateerd aan de lineaire predictor $\eta_{jpshm}$.
Door deze $\log$ link kunnen we de individuele termen van de lineare predictor interpreteren als relatieve effecten.
$$\log\mu_{jpshm} = \eta_{jpshm}$$
En nu wordt het interessant want deze lineaire predictor hangt af van vijf termen:
- $\beta_s$: het effect van stratum $s$.
- $\beta_p$: het verschil tussen periode $p$ en de referentieperiode.
- $b_h$: het effect van hok $h$. Deze effecten komen uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_h$.
- $b_m$: het effect van meetpunt $m$. Deze effecten komen uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_m$.
- $b_{j,s}$: het effect van jaar $j$ in stratum $s$. Dit effect modelleert een eerste orde random walk per stratum. Het verschil tussen twee opeenvolgende jaren komt uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_j$.
$$\eta_{jpshm} = \beta_s + \beta_p + b_h + b_m + b_{j,s}$$ $$b_h \sim \mathcal{N}(0, \sigma^2_h)$$ $$b_m \sim \mathcal{N}(0, \sigma^2_m)$$ $$b_{j,s} - b_{j-1, s} = \Delta b_{j,s} \sim \mathcal{N}(0, \sigma^2_j)$$
De random walk $b_{j,s}$ is de term die ons het meest interesseert aangezien deze modelleert hoe de aantallen in de tijd wijzigen. Het komt er op neer dat de aantallen in jaar $j$ en stratum $s$ een factor $e^{\Delta b_{j,s}}$ verschillen van het voorgaande jaar in datzelfde stratum. De variantie $\sigma^2_j$ bepaalt hoe groot de sterkste schommelingen tussen de opeenvolgende jaren kunnen zijn. Het model laat toe dat elk stratum zijn eigen trend heeft.
Statistische verdelingen
Een standaard distributie voor aantallen is de Poisson verdeling. Deze heeft als kenmerk dat de variantie gelijk is aan het gemiddelde ($\sigma^2 = \mu$). Bij veel ecologisch gegevens zien we dat de variantie in de praktijk groter is dan het gemiddelde, een fenomeen dat we overdispersie noemen. In dat geval kunnen we een negatief binomiale verdeling gebruiken. Deze heeft als variantie $\sigma^2 = \mu + \mu^2/\nu$, waarbij de parameter $\nu$ is een maat voor overdispersie. Merk op dat als de overdispersieparameter zeer groot wordt ($\nu = \infty$), dan wordt de term $\mu/\nu = 0$ en bijgevolg reduceert de negatief binomiale verdeling dan tot een Poisson verdeling.
Sommige ecologische gegevens vertonen een "overmaat" aan nullen, dat zijn meer nulwaarnemingen dan de Poisson of negatief binomiale verdeling kan modelleren.
In dergelijke gevallen kunnen we overschakelen naar de zero-inflated versie van deze verdelingen.
Deze hebben een parameter $\pi$ die een maat is voor de kans op een overmatige nul.
Bij het modelleren zullen we in eerste instantie deze twee parameters instellen op $\nu = \infty$ en $\pi = 0$, m.a.w. geen overdispersie en geen overmaat aan nullen zodat we een Poisson verdeling krijgen.
Vervolgens gaan we na of er voldoende aanwijzingen zijn voor overdispersie of een overmaat aan nullen.
In het geval van overdispersie laten we $\nu$ door het model schatten waardoor we overgaan van een Poisson naar een negatief binomiaal.
In het geval van een overmaat aan nullen laten we $\pi$ door het model schatten waardoor we overgaan van een Poisson naar een zero-inflated Poisson.
Soms hebben we zowel aanwijzingen voor overdispersie als een overmaat aan nullen.
In die gevallen kiezen we de negatief binomiaal of zero-inflated Poisson naargelang welke de sterkste aanwijzingen heeft.
Vervolgens gaan we na of we de andere parameter ook nog een probleem vormt.
Zo ja, gaan we over naar een zero-inflated negatief binomiaal waarbij het model zowel $\nu$ als $\pi$ zal schatten.
Lineaire trends {#s:lineairetrend}
Lineaire trends veronderstellen dat er een constante wijziging is over de volledige looptijd.
Het model dat we hiervoor gebruiken is nagenoeg identiek aan het niet-lineaire model.
Het enige verschil zit in de lineaire predictor waar we de first order random walk ($b_{j,s}$) vervangen door een lineaire trend per stratum ($\beta_{s1} j$).
$$\eta_{jpshm} = \beta_s + \beta_p + b_h + b_m + \beta_{s1} j$$
Het lineaire model is een vereenvoudiging van het niet-lineaire model.
We kunnen deze modellen met elkaar vergelijken op basis van het Wantanabe-Akaike Information Criterion (WAIC) [@gelmanUnderstandingPredictiveInformation2014].
De WAIC waarde daalt naarmate het model de gegevens beter kan beschrijven en stijgt wanneer het model complexer wordt.
Als we modellen met elkaar vergelijken op basis van WAIC, zal het model met de laagste WAIC de beste mix zijn tussen een goede beschrijving van de gegevens en een zo eenvoudig mogelijk model.
We beschouwen de trend als lineair wanneer het lineaire model de laagste WAIC heeft.
Wanneer de WAIC van het niet-lineaire model minder dan 2 eenheden lager is dan deze van het lineaire model, beschouwen de trend als mogelijk niet-lineair.
Pas wanneer de WAIC van het niet-lineaire model duidelijk lager is dan het lineaire model, stellen we dat de trend niet-lineair is.
In dat geval de gebruiker moet de gerapporteerde lineaire trend met de nodige voorzichtigheid interpreteren.
Kijk zeker naar de bijhorende niet-lineaire trend vooraleer de cijfers te interpreteren.
Resultaten per driejarige cyclus
Zoals aangegeven in §\@ref(s:steekproef) zullen we een bepaald hok in principe om de drie jaar herbezoeken. Hierdoor krijgen we drie sets van hokken. Er bestaat een kans dat de schatting voor een bepaald jaar beïnvloed is door de set van hokken die in dat jaar onderzocht worden. Om dit effect uit te schakelen, analyseren we de gegevens tevens op basis van de driejarige cyclus $c$ i.p.v. jaar $j$.
De eerste cyclus omvat de eerste drie jaar sinds de start van het meetnet (2007-2009). De volgende cyclus omvat telkens de volgende drie jaar aansluitend op de vorige cyclus. Aangezien we alle beschikbare gegevens gebruiken bij de analyse bevat de laatste cyclus mogelijk minder dan drie jaar. Dit is duidelijk zichtbaar doordat het laatste jaar van de cyclus op dat ogenblik in de toekomst ligt. Als bijvoorbeeld 2019 het meest recente jaar met gegevens is, dan is de laatste cyclus 2019-2021.
Hieronder geven we aan op welke manier we de eerste beschreven modellen aanpassen.
$$Y_{cpshm} \sim ZINB(\mu_{cpshm}, \nu, \pi)$$ $$\log\mu_{cpshm} = \eta_{cpshm}$$
Niet-lineaire trend
$$\eta_{cpshm} = \beta_s + \beta_p + b_h + b_m + b_{c,s}$$ $$b_{c,s} - b_{c-1, s} = \Delta b_{c,s} \sim \mathcal{N}(0, \sigma^2_c)$$
$b_{c,s}$: het effect van driejarige cyclus $c$ in stratum $s$. Dit effect modelleert een eerste orde random walk per stratum. Het verschil tussen twee opeenvolgende driejarige cycli komt uit een Gaussiaanse verdeling met gemiddelde 0 en variantie $\sigma^2_c$.
Lineaire trend
$$\eta_{cpshm} = \beta_s + \beta_p + b_h + b_m + \beta_{s1} c$$
Modellen fitten
We fitten de statische modellen in R [@R] met het INLA package [@rueINLAFunctionsWhich2009]. INLA gebruikt een Bayesiaanse benadering om de modellen te fitten. Daarom moeten we priors specificeren voor de parameters en hyperparameters.
- $\beta_s$, $\beta_p$ en $\beta_{s1}$ krijgen een Gaussiaanse prior $\mathcal{N}(0, 5)$.
- $\sigma^2_h$ krijgt een
penalised complexityprior [@simpson.rue.ea2017PenalisingModelComponent] zodat $Prob(\sigma_h > 1) = 0.01$. - $\sigma^2_m$ krijgt een
penalised complexityprior zodat $Prob(\sigma_m > 1) = 0.01$. - $\sigma^2_j$ krijgt een
penalised complexityprior zodat $Prob(\sigma_j > 0.5) = 0.01$. - $\sigma^2_c$ krijgt een
penalised complexityprior zodat $Prob(\sigma_c > 0.5) = 0.01$. - $\nu$ krijgt een Gamma prior zodat $\nu \sim \Gamma(1/7, 1/7)$.
- $\pi$ wordt intern geparametriseerd als $\theta = \log\frac{\pi}{1 - \pi}$ met Gaussiaanse prior $\theta ~ \sim \mathcal{N}(-1, 5)$.
Gemiddeld aantal dieren per meetpunt
Een gemiddeld hok heeft als effect $b_h = 0$ en een gemiddeld meetpunt $b_m = 0$.
In de referentieperiode is $\beta_p = 0$.
In deze gevallen vereenvoudigt de lineaire predictor tot het effect van stratum $s$ en de trend in dat stratum:
$$\eta_{js} = \beta_s + b_{j,s}$$
Door de effecten van de strata te vermenigvuldigen met hun stratumgewicht ($\gamma_s$) krijgen we een schatting voor Vlaanderen.
$$\eta_{j} = \sum_s(\gamma_s\beta_s + \gamma_sb_{j,s})$$
Om de schatting van het gemiddelde aantal in Vlaanderen te krijgen, moeten we de lineaire predictor terug omzetten van de $\log$ schaal naar de natuurlijke schaal en corrigeren voor de eventuele overmaat aan nullen.
We krijgen dan:
$$E[Y_j] = (1 - \pi) e^{\eta_j} \prod_s(e^{\gamma_s \beta_s} e ^ {\gamma_s b_{j,s}})$$
Vergelijken van jaren
We kunnen twee jaren $a$ en $b$ met elkaar vergelijken door hun lineaire predictoren voor het gemiddeld aantal dieren van elkaar af te trekken.
$$\eta_a - \eta_b = \sum_s(\gamma_s\beta_s + \gamma_sb_{a,s}) - \sum_s(\gamma_s\beta_s + \gamma_sb_{b,s})$$
Aangezien het globale stratumeffecten en de stratumgewichten niet wijzigen in de tijd, kunnen we dit vereenvoudigen tot
$$\eta_a - \eta_b = \sum_s\gamma_s(b_{a,s} -b_{b,s})$$
De wijziging in Vlaanderen is het gewogen gemiddelde van de wijzigingen in de strata. Na omzetting van de $log$ schaal naar de natuurlijke schaal krijgen we het relatieve aantal $I_{a|b}$ van jaar $a$ t.o.v. jaar $b$
$$I_{a|b}=\frac{e^{\nu_a}}{e^{\nu_b}} = \prod_s\left(\frac{e^{\gamma_sb_{a,s}}}{e^{\gamma_sb_{b,s}}}\right)$$
Wanneer we een jaar met zichzelf vergelijken krijgen we per definitie $I_{a|a} = 1 = 100\%$. Wanneer we meerdere jaren met eenzelfde referentiejaar vergelijken krijgen we een indexwaarde: het relatieve verschil van elk jaar t.o.v. van een bepaald referentiejaar.
Gewicht van de strata {#s:stratumgewicht}
In het vorige onderdeel hebben gebruikt gemaakt van stratumgewichten $\gamma_s$ zonder deze te definiëren. We bepalen het gewicht van een stratum op basis van drie kenmerken: het totaal aantal hokken van het stratum in Vlaanderen ($N_s$), het aantal onderzochte hokken in het stratum ($T_s$) ongeacht of ze al dan niet relevant zijn voor de soort en het aantal relevante hokken voor de soort in het stratum ($R_s$).
Het aandeel relevante hokken per stratum ($R_s/T_s$) is een goede maat van de frequentie waarmee een soort aanwezig is binnen een stratum. Wanneer een soort in nagenoeg alle onderzochte hokken van het stratum waargenomen wordt is $R_s/T_s \simeq 1$. Is de soort zeer zeldzaam (komt slechts in een paar van de onderzochte hokken voor) dan is $R_s/T_s \simeq 0$.
Het basisgewicht $\gamma_{sb}$ voor stratum $s$ is het aandeel relevante hokken vermenigvuldigd met het totaal aantal hokken in het stratum. Dit is een schatting van het totaal aantal hokken in het stratum waar de soort voldoende frequent voorkomt.
$$\gamma_{sb} = \frac{R_s}{T_s}N_s$$
Om makkelijker te kunnen rekenen delen we de basisgewichten door hun som zodat de stratumgewichten sommeren tot 1.
$$\gamma_s = \frac{\gamma_{sb}}{\sum_s\gamma_{sb}}$$
Merk op dat het aantal relevante hokken per stratum ($R_s$) soortafhankelijk is. Bijgevolg zijn de stratumgewichten eveneens soortafhankelijk.
Samengestelde indices {#s:samengesteld}
We berekenen een aantal samengestelde indices die de trends voor een groep van soorten aggregeren. Deze indices zijn het meetkundig gemiddelde van de verschillen tussen jaren voor alle soorten van de groep. Bij de berekening maken we gebruik van het trucje dat we een product kunnen schrijven als een som van logaritmes die we nadien terug exponentiëren. Een meetkundig gemiddelde in de natuurlijke schaal kunnen we dat schrijven als een rekenkundig gemiddelde in de log-schaal.
$$\sqrt[n]{\prod_{i+1}^na_i} = \exp\left(\frac{\sum_{i = 1}^n\log a_i}{n}\right)$$ Het voordeel aan deze techniek is dat de schattingen van de paarsgewijze verschillen tussen de jaren reeds beschikbaar zijn in de log-schaal. Bovendien beschikken we tevens over hun variantie in de log-schaal. Dit laat ons toe om makkelijk het betrouwbaarheidsinterval te berekenen aan de hand van een paar vuistregels.
- Bij onafhankelijke variabelen geldt dat de variantie van hun som gelijk is aan de som van de varianties. $\sigma^2_{\sum X_i} = \sum \sigma^2_{X_i}$
- De variantie van het product van een variabele en een constante is gelijk aan de variantie van de variabele vermenigvuldigd met die constante. $\sigma^2_{aX} = a \sigma^2_X$
- Het model veronderstelt dat de individuele modelparameters een Gaussiaanse verdeling volgen. Hierdoor volgen de paarsgewijze verschillen tussen de jaren en hun gemiddelde (telkens in de log-schaal) een Gaussiaanse verdeling. Aangezien we zowel het gemiddeld als de variantie van deze verdeling kennen, kunnen we hieruit de gewenste kwantielen voor de betrouwbaarheidsintervallen berekenen (nog steeds in de log-schaal).
- De log-transformatie is een monotoon stijgende functie aangezien $\log(x)$ steeds groter wordt wanneer $x$ groter wordt. Een kenmerk van een monotoon stijgende transformatie is dat ze de volgorde behouden: het 5% kleinste element zal ook na de transformatie het 5% kleinste element zijn. Waardoor we de kwantielen van de betrouwbaarheidsintervallen zonder probleem kunnen terugrekenen naar de natuurlijke schaal.
Voorstelling van de gegevens
Onzekerheid {#s:onzekerheid}
Alle resultaten zijn gebaseerd op een steekproef en op de waarnemingen zit onvermijdelijk een zekere meetfout. Vandaar dat we bij de puntschattingen tevens een betrouwbaarheidsinterval weergeven. In de tekst gebruiken we het 90% (5%; 95%) interval waarbij er 5% kans is dat de werkelijke waarde kleiner is dan de ondergrens en 5% dat ze groter is dan de bovengrens. Het 90% interval is iets smaller dan het traditionele 95% (2.5%; 97.5%) interval. Door een smaller interval te kiezen zullen we sneller uitspraken kunnen doen, waardoor de kans kleiner wordt dat we ten onrechte stellen dat er geen effect is. De prijs die we hiervoor betalen is dat de kans dat we ten onrechte stellen dat er een significant effect is, stijgt van 5% naar 10%. De doelstelling van de algemene broedvogelmonitoring is zo spoedig mogelijk detecteren wanneer er iets aan de hand is de broedvogels. Vanuit dat oogpunt is het vermijden van vals negatieve signalen (ten onrechte stellen dat er niets aan de hand is) belangrijker dan het vermijden van vals positieve signalen (ten onrechte stellen dat er iets aan de hand is).
Wanneer we het interval op een figuur (bijvoorbeeld fig. \@ref(fig:trendklasse)) weergeven, vullen we het 90% interval aan met een 60% (20%; 80%) interval en een 30% (35%; 65%) interval. @britton.fisher.ea1998InflationReportProjections waren de inspiratie voor deze manier van weergeven. Het 30% interval vormt het donkerste deel van het interval en geeft de meeste waarschijnlijke locatie van de werkelijke waarde weer. Naarmate het interval lichter wordt, neemt de onzekerheid toe.
Opdeling van de effecten in een aantal klassen {#s:trendklasse}
Bij het niet-lineaire model berekenen we alle paarsgewijze relatieve verschillen tussen de verschillende jaren. Bij het lineaire model hebben we de gemiddelde jaarlijkse relatieve verandering. Deze laatste rekenen we tevens om naar de totale wijziging over de looptijd van het meetnet omdat dit eenvoudiger te interpreteren is. Vergelijk een daling met -5% per jaar of een daling met -50% over 15 jaar. Deze laatste klinkt dramatischer door het grotere cijfer, terwijl -5% per jaar overeenkomt met -54% over 15 jaar. De totale wijziging over de looptijd van het lineaire model is tevens vergelijkbaar met de wijziging tussen het eerste en laatste jaar van het niet-lineaire model.
Om de interpretatie makkelijker te maken, delen we de wijzigingen op in tien klassen door hun 90% interval te vergelijken met een referentie, ondergrens en bovengrens. Een effect is significant wanneer de referentie buiten het 90% interval ligt. We spreken over een toename (afname) als het interval volledig boven (onder) de referentie ligt. Niet significante effect is ook informatief wanneer het bijhorende interval voldoende smal is. Bijvoorbeeld als het interval volledig tussen de ondergrens en bovengrens ligt. In dat geval kunnen we stellen dat het effect niet significant en klein is, het immers zeker minder sterk dan de ondergrens en minder sterk dan de bovengrens. Dergelijk effect krijgt de naam stabiel. Heeft het effect een breed interval dat zowel de boven- als ondergrens bevat, spreken we over een onduidelijk effect. Daarnaast is er nog de mogelijkheid dat het interval zowel de bovengrens (ondergrens) als de referentie bevat maar niet de ondergrens (bovengrens). Dan spreken we over een mogelijke toename (mogelijke afname). We kunnen de boven- en ondergrens eveneens gebruiken om een verder onderscheid te maken binnen de significante effecten. Een interval volledig boven (onder) de bovengrens (ondergrens) wordt dan een sterke toename (sterke afname). Een interval volledig tussen de referentie en de de bovengrens (ondergrens) wordt dan een matige toename (matige afname). Een interval dat de referentie niet bevat maar wel de bovengrens (ondergrens) blijft een toename (afname). Merk op dat de indeling volledig gebaseerd is op de onzekerheid rond het effect en niet op de puntschatting van het effect zelf. We vatten de opdeling met bijhorende afkortingen en regels samen in tabel \@ref(tab:regels). Figuur \@ref(fig:trendklasse) geeft een grafische voorstelling waarbij we de afkortingen in combinatie met aangepaste symbolen gebruiken. De afkortingen zelf zijn te fijn om als symbool te gebruiken.
(ref:regels) Overzicht van de benamingen van de tien effectklassen met hun afkorting en de regels. $R$: referentie, $L$: ondergrens, $B$: bovengrens, $l$: ondergrens van het 90% interval, $b$: bovengrens van het 90% interval. $L < R < B$ en $l < b$.
tribble( ~ benaming, ~afkorting, ~regels, "sterke toename", "`++`", "$B < l$", "toename", "`+`", "$R < l < B$ en $B < b$", "matige toename", "`+~`", "$R < l < B$ en $b < B$", "stabiel", "`~`", "$L < l < R$ en $R < b < B$", "matige afname", "`-~`", "$L < l < R$ en $b < R$", "afname", "`-`", "$l < L$ en $L < b < R$", "sterke afname", "`--`", "$l < L$", "mogelijke toename", "`?+`", "$L < l < R$ en $B < b$", "mogelijke afname", "`?-`", "$l < L$ en $R < b < B$", "onduidelijk", "`?`", "$l < L$ en $B < b$" ) %>% kable( caption = "(ref:regels)", format = "pandoc", escape = FALSE, align = "lcc" )
(ref:trendklasse) Voorbeeld van de tien mogelijke interpretaties van een effect door het 90% interval te vergelijken met een referentie, ondergrens en bovengrens.
ggplot(effect, aes(x = x_int, y = y, ymin = lcl, ymax = ucl, link_sd = s)) + geom_hline(yintercept = 0, linetype = 2) + geom_hline(yintercept = c(-1, 1), linetype = 3) + geom_rect( aes(xmin = x_int - 0.4, xmax = x_int + 0.4, fill = klasse), show.legend = FALSE ) + geom_point(size = 6, colour = inbo_hoofd) + geom_text( aes(label = klasse), hjust = 0.5, vjust = 0.5, colour = "white" ) + scale_x_continuous(breaks = effect$x_int, labels = effect$x) + scale_y_continuous( "effect", breaks = c(-1, 0, 1), labels = c("ondergrens", "referentie", "bovengrens") ) + scale_fill_manual(values = kleurgradient) + coord_flip() + theme(axis.title.y = element_blank())
Uiteraard hangt de opdeling sterk af van de keuze van de boven- en ondergrens. De soortenmeetnetten voor de Natura 2000 monitoring streven er naar om een daling in populatiegrootte met -25% over 25 jaar tijd vast te kunnen stellen. Hierbij wordt een wijziging in de populatiegrootte van -25% als belangrijk aanzien. Daarom hanteren we voor de algemene broedvogelmonitoring ook -25% als ondergrens, zowel bij de paarsgewijze verschillen tussen de jaren als de lineaire wijziging over de looptijd van het meetnet. Een van daling van -25% komt overeen met aantallen die nog drie kwart van de uitgangssituatie bedragen. Een even sterke wijziging in de omgekeerde richting zorgt er voor dat de aantallen toenemen tot vier derde van de uitgangssituatie, of een toename met +33% wat we als bovengrens gebruiken.
Bij de samengestelde indices hebben we aangepaste grenswaarden nodig. Veronderstel een samengestelde index op basis van $n$ soorten. We berekenen het rekenkundig gemiddelde in de log-schaal, dan is de variantie van dit gemiddelde een factor $n$ kleiner dan de som van de varianties. De breedte van een betrouwbaarheidsinterval hangt samen met de standaard afwijking, wat de vierkantswortel van de variantie is. Hierdoor zullen de breedtes van de betrouwbaarheidsintervallen van de samengestelde index een factor $\sqrt{n}$ kleiner zijn. Vandaar dat we de grenswaarden van de indices tevens aanpassen door ze met een factor $\sqrt{n}$ te verkleinen. Elke samengestelde index heeft zijn eigen soortenlijst met een variabele aantal soorten. Daarom zullen we bij elke samengestelde index zijn aangepaste grenswaarden vermelden.
Overzicht van de lineaire trends
Tabellen \@ref(tab:driejaarlijkse-trends) en \@ref(tab:jaarlijkse-trends) geven een overzicht van de lineaire trends voor elke soort. We hebben de trends gesorteerd volgens opdeling van de klassen en binnen de klasse volgens puntschatting van de trend. Hierdoor start de tabel met de soorten met de sterkste positieve trends. Bij de volgende soorten zal de trend minder sterk worden tot we aan de soorten met een stabiele trend komen. Daarna volgen de soorten met een negatieve trend waarbij de trend steeds sterker negatief wordt. We sluiten de tabel af met de soorten met een mogelijke of onduidelijke trend. De tabel bevat de trend uitgedrukt als een jaarlijkse wijziging en als een wijziging over de volledige looptijd van het meetnet. Verder bevat de tabel de opdeling van de trend in klassen en een indicatie of de trend al dan niet lineair is (zie §\@ref(s:lineairetrend)). De naam van de soort is een snelkoppeling naar de detail van de soort zelf. We raden aan om hiervan gebruikt te maken bij de interpretatie van niet-lineaire trends.
Evolutie van gemiddelde aantallen per soort
Deze figuur geeft de evolutie van de gemiddelde aantallen per meetpunt weer volgens het niet-lineair model. Indien het model lineair is, zal het patroon van de niet-lineaire trend dicht bij een lineaire trend liggen. De lijn bevat de puntschatting van het gemiddelde aantal in elk jaar. Rond de lijn tonen we het 30%, 60% en 90% interval (zie §\@ref(s:onzekerheid)). In de online versie van het rapport is de figuur interactief. Wanneer de gebruiker met de muis over de figuur gaat verschijnt er in de buurt van de lijn een pop-up met de exacte schatting van dat jaar inclusief het 90% betrouwbaarheidsinterval en het jaartal.
Paarsgewijze vergelijking van jaren
Deze informatie hebben we telkens in twee figuren samengevat. De eerste figuur geeft het verschil van elk jaar t.o.v. 2007, het jaar waarin de metingen gestart zijn. Voor 2007 geven we geen cijfer omdat het per definitie 1 is en geen informatie bevat. Het symbool geeft de puntschatting van het relatieve verschil t.o.v. 2007 weer. De vorm van het symbool geeft de opdeling van de sterkte van het effect weer (zie fig. \@ref(fig:trendklasse)). Rond elke puntschatting tonen we het 30%, 60% en 90% interval (zie §\@ref(s:onzekerheid)). De horizontale streepjeslijn geeft de referentie van 0% verschil weer. De horizontale puntlijnen geven de bovengrens (+33%) en ondergrens (-25%) weer. Deze laten toe om de vlot zelf de betrouwbaarheidsintervallen te vergelijken met de referentie, boven- en ondergrens. De figuur heeft twee y-assen. Beide assen geven dezelfde informatie weer, enkel de formattering van de labels is anders. De linkeras toont procentuele verschillen terwijl de rechteras de relatieve verschillen toont. Een procentueel verschil van +50% is equivalent met een relatieve verschil (verhouding) met een factor 1,5.
We stellen vast dat veel gebruikers ook andere jaren met elkaar willen vergelijken. Een correcte vergelijking is enkel mogelijk indien we een van deze jaren als referentie gebruiken. Voor elk jaar een afzonderlijke figuur maken, zou het rapport onoverzichtelijk groot maken. Om de vergelijkingen toch mogelijk te maken hebben we alle paarsgewijze verschillen tussen de jaren in een raster weergeven. Elke rij in het raster staat voor een ander referentiejaar. De kolommen geven de verschillen van een bepaald jaar weer t.o.v. de verschillende referentiejaren. Op de diagonaal staan geen waarden omdat we daar een jaar met zichzelf vergelijken. De kleur van de symbolen geeft de sterkte van het verschil (gebaseerd op de puntschatting). Zwakke verschillen zijn grijs, sterke positieve verschillen rood, sterke negatieve verschillen blauw. De vorm van de symbolen geeft zicht op de sterkte en onzekerheid van het effect (zie tab. \@ref(tab:regels)). Wanneer een rij volledig rood (blauw) is, zijn alle verschillen met dit referentiejaar positief (negatief) m.a.w. dit is het referentiejaar met de laagste (hoogste) aantallen. Wanneer een kolom volledig rood (blauw) is, zijn alle verschillen van dit jaar t.o.v. alle referentiejaren negatief (positief) m.a.w. dit is het jaar met de hoogste (laagste) aantallen. Clusters van punten met een gelijkaardige kleur geven periodes aan waarin de aantallen geleidelijk wijzigen. In de online versie van dit rapport is deze figuur interactief. Wanneer de gebruiker met de muis over de figuur gaat verschijnt er in de buurt van de lijn een pop-up met de exacte schatting van dat punt inclusief het 90% betrouwbaarheidsinterval en het jaartal en referentiejaar.
Reproduceerbaarheid en traceerbaarheid van de gegevensverwerking
Versiebeheer
De waarnemers voeren hun waarnemingen via een webapplicatie toe aan de databank.
De analyse start met het importeren van de relevante gegevens uit de databank aan de hand van een R-script.
De bekomen dataset bewaren we als een collectie van tekstbestanden met behulp van het git2rdata package [@git2rdata].
Zowel deze code als de tekstbestanden bewaren we onze versiebeheer met git.
Git is een gratis en open source gedistribueerd versiebeheersysteem dat is ontworpen om alles van kleine tot zeer grote projecten snel en efficiëntie te beheren.
De code zit vervat in het abvanalysis package [@abvanalysis], dat vrij online beschikbaar is.
De tekstbestanden met de gegevens zijn momenteel niet publiek toegankelijk.
De ruwe gegevens worden geaggregreerd per hok en met drie jaar vertraging publiek ontsloten via GBIF [@vermeersch.anselin.ea2018ABVCommonBreeding].
Reproduceerbaarheid
Om de reproduceerbaarheid te garanderen starten we de analyses steeds vanaf de tekstbestanden.
Aan de hand van code in het abvanalysis package definiëren we hoe we de tekstbestanden omzetten naar de verschillende analyses (§\@ref(s:trendberekening)).
Dit resulteert in een resem op zichzelf staande analyseobjecten van het n2kanalysis package [@n2kanalysis].
Dergelijk analyseobject bevat naast de nodige gegevens en de modeldefinitie tevens de nodige metadata zoals alle gebruikte R packages, inclusief hun versie, en een verwijzing naar de voorafgaande analyses waarvan deze analyse afhangt.
Zo verwijst de analyse van een samengestelde index (§\@ref(s:samengesteld)) naar de trendberekening van de bijhorende soorten (§\@ref(s:trendberekening)).
De individuele trendberekeningen verwijzen naar de "analyse" die de gegevens importeerde uit de databank.
Deze laatste bevat de nodige links naar de tekstbestanden en hun versie.
Traceerbaarheid {#s:traceerbaarheid}
De n2kanalysis objecten hebben elk twee data hashes.
Elk data hash is een reeks van 40 hexadecimale^[Hexadecimaal betekent letterlijk zestientallig. Het is een talstelsel waarbij niet, zoals gebruikelijk, met tien cijfers wordt gewerkt, maar met zestien cijfers. De cijfers 0 t.e.m. 9 worden daarom uitgebreid met a (=10) t.e.m. f (=15).] cijfers die het resultaat zijn van de cryptografische hashfunctie SHA-1.
Deze hashfunctie heeft een aantal belangrijke eigenschappen:
- Ze zetten elke invoer om naar een uitvoer met vaste lengte (40 hexadecimale cijfers).
- De uitvoer is stabiel: als je de hash van een bepaalde invoer opnieuw berekent krijg je steeds dezelfde uitvoer.
- Het is niet mogelijk om de invoer te reconstrueren op basis van de uitvoer.
- Eender welke kleine wijziging aan de invoer resulteert in een sterke wijziging van de uitvoer.
- De kans dat twee verschillende invoeren dezelfde uitvoer opleveren is zeer klein.
De eerste data hash van het analyseobject is gebaseerd op alle informatie die gekend is op het moment dat we het analyseobject definiëren en die nooit zal wijzigen tijdens de analyse. Denk hierbij aan de definitie van de analyse, de gegevens, de soort, ... Gezien de eigenschappen van de hashfunctie kunnen we deze data hash gebruiken om ondubbelzinnig te verwijzen naar een specifieke analyse (inclusief de gebruikte gegevens.)
De tweede data hash van het analyseobject baseren we enerzijds op de eerste data hash en anderzijds van alle onderdelen van het analyseobject die wijzigen in de loop van de analyse. Dit is o.a. het resultaat van de statistische analyse, de gebruikte software, ... De analyse opnieuw uitrekenen met software van een andere versie zal de tweede data hash aanpassen.
De combinatie van deze twee data hashes laat enerzijds toe om naar een specifieke versie van de analyse te verwijzen. Anderzijds bieden ze een garantie over de inhoud van het analyseobject. In het geval van twijfel over een analyse kunnen we teruggrijpen naar de analyseobject in kwestie. De data hashes bewijzen dan dat we het correcte analyseobject hebben.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.