Doorontwikkeling provinciale natuurindicatoren

In opdracht van BIJ12 en de Werkgroep Monitoring Natuur heeft het CBS enkele jaren geleden natuurindicatoren op provinciaal (en lager) schaalniveau ontwikkeld. Deze natuurindicatoren geven een gemiddeld beeld van de trend in de populaties van hele soortgroepen, bijvoorbeeld alle broedvogels of alle dagvlinders per provincie, of van kenmerkende soorten per ecosysteemtype, zoals bos, heide of open duin. 

De gegevens voor de provinciale natuurindicatoren zijn in meerderheid afkomstig uit de verschillende meetprogramma’s van het Netwerk Ecologische Monitoring. Deze meetprogramma’s zijn destijds opgezet om met name landelijke trends per soort op te leveren. Representativiteit en statistische robuustheid zijn daarom op landelijk niveau grotendeels geborgd. Op provinciaal niveau is dat voor een gedeelte van de soorten niet het geval. Om grip te houden op de kwaliteit van de provinciale indicatoren gebaseerd op provinciale trends is nadere inspectie van de onderliggende trends per soort noodzakelijk. 

Aanvankelijk werd hiertoe een plausibiliteitstoets ontwikkeld, waarbij het aantal meetpunten, de statistische significantie van de trend en de richting van de trend in vergelijking met het landelijke beeld en dat van aanpalende provincies, meewogen in het besluit de trend al dan niet te gebruiken voor de betreffende indicator. Het aantal meetpunten (de “teldekking”) woog dus wel mee, maar was niet diskwalificerend op zichzelf; als de andere onderdelen van de plausibiliteitstoets ‘de goede kant’ op wezen, dan kon een trend die bepaald was met zeer magere teldekking gewoon meedoen.

Die teldekking is recent echter dusdanig cruciaal gebleken om een stijging of daling van een trend vast te kúnnen stellen, dat de plausibiliteitstoets alleen niet volstaat. In 2025 is in opdracht van BIJ12 daarom ingezet op verdere ontwikkeling van een methode om vast te stellen of een trend op voldoende meetpunten is gebaseerd om opname in de provinciale natuurindicatoren te rechtvaardigen. Het onderzoek heeft geresulteerd in vuistregels voor teldekking. Toepassing van deze vuistregels bij de actualisatie van de indicatoren heeft geleid tot uitsluiting van een beperkt aantal soorten per indicator. 

Al eerder was gebleken dat de indicatoren door de bank genomen vrij robuust zijn tegen het afvallen van een beperkte set aan soorten. De meeste indicatoren zijn als gevolg van het afvallen van een paar soorten dan ook niet heel erg veranderd wat trendrichting en trendbeoordeling betreft t.o.v. de indicatorenset van vorige ronde, al zijn er zeker uitzonderingen. 

Dit rapport beschrijft de aanpak van het onderzoek naar de teldekking en het formuleren van de vuistregels voor opname in de provinciale natuurindicatoren.

In dit onderzoek is een simulatiestudie uitgevoerd om te bepalen vanaf welke teldekking (aantal jaren, plots en observaties) het mogelijk is om een trend van ten minste 5% toename of afname betrouwbaar te detecteren. Aangezien de kwaliteit van trends sterk afhankelijk is van de teldekking, is het van belang om te begrijpen hoe deze factoren de statistische uitkomsten beïnvloeden, vooral bij provinciale of andere kleinschalige trends waar de dekking vaak beperkter is.

Genereren van synthetische data

Synthetische datasets werden gegenereerd door variaties in het aantal jaren, plots en observaties. De jaren werden altijd vanaf 1980 genomen, waarbij een minimale teldekking wordt meegegeven van tenminste 10% bij minder dan 30 plots, of tenminste 5% bij 30 plots of meer. Waarbij er minstens twee observaties per plot moeten zijn (anders heeft het plot geen invloed op de trend), en minstens één observatie per jaar. De observaties worden vervolgens aangevuld tot het minimum bereikt is en volgens de andere meegegeven parameters (o.a. skewness, schaarste).

De aantallen per observatie worden getrokken uit een Poisson distributie (benodigde verdeling voor het TRIM model waarmee trends berekend worden, zie volgende sectie) op een gesimuleerde toenemende of afnemende trend van 5% per jaar. Per dataset wordt ook nog gevarieerd in het percentage hoge of lage aantallen over de plots. Daarnaast wordt ook een percentage variatie per plot meegenomen, waarbij het startpunt van de getrokken aantallen kan variëren, bovenop de hoge of lage aantallen. 

Tot slot is er nog rekening gehouden met mogelijke verschillen in de verdeling van de observaties over de jaren en de plots. Zo kunnen er voornamelijk tellingen zijn geweest in het begin of aan het eind van de monitoringsdata. Dit noemen we hier ‘skewness’, waarbij per dataset is meegegeven of de jaren die geteld zijn per plot voornamelijk in het begin van de tijdreeks (linkse skewness), in het midden van de tijdreeks (symmetrische skewness) aan het eind van de tijdreeks (rechtse skewness) of aan de randen van de tijdreeks (u-shaped skewness) zitten. 

Ook kan er variatie zitten in de teldekking per plot, met een aantal goed getelde plots en een aantal minder goed getelde plots over de jaren. De variatie in deze teldekking is ook nog meegegeven aan de gesimuleerde dataset in de vorm van een percentage voor het aantal minder goed getelde plots, en noemen we ‘schaarste’. Bij datasets met een laag aantal jaren en plots met een lage teldekking werkt deze berekening niet, en is daar dus niet meegenomen.

De variatie in de gebruikte parameters om synthetische datasets te genereren Tabel 1, in totaal zijn er meer dan 300,000 datasets genereerd. Een voorbeeld van de teldekking van een gesimuleerde dataset is te zien in Figuur 1. 

2.2 Voorbeeld van de teldekking van een gesimuleerde dataset van 30 jaar met 40 plots met een toenemende trend, met lage aantallen (met 50% variatie per plot), een rechtse skewness, en 75% schaarste. De blauwe vlakken zijn de genereerde observaties (met aantallen van 1 tot en met 50) en de grijze vlakken de missende data.

Resultaten simulatie data

Over alle simulatiedatasets hebben we populatietrends berekend met het programma TRIM (Trends & Indices for Monitoring data; Pannekoek & van Strien 2005; beschikbaar in R als package rtrim; Bogaart et al., 2024). Dit programma is speciaal door het CBS ontwikkeld voor de analyse van telgegevens op vaste meetpunten en wordt zowel in Nederland als het buitenland veel gebruikt in de natuurmonitoring. TRIM maakt gebruik van een Poisson-regressiemodel om veranderingen in aantallen door de tijd te modelleren. Daarbij wordt gecorrigeerd voor overdispersie (extra variatie tussen tellingen) en seriële autocorrelatie (afhankelijkheid tussen aantallen uit opeenvolgende jaren).

Hierbij hebben we gekeken of de berekende trendklassen met TRIM overeenkomen met de meegegeven trendklassen (5% toe- of afname). Op basis van de simulatiedatasets konden we vaststellen dat de juistheid van de berekende trendklasse sterk afhankelijk is van het aantal jaren, aantal plots, en aantal observaties waarop de data gebaseerd is. De andere parameters waarmee we hadden gevarieerd hadden minimale invloed. Maar bij een korte tijdsreeks is de trend wel beter te detecteren bij gemiddeld hogere aantallen per observatie en bij een u-shape verdeling (Figuur 2). Er is geen verschil in detectie van een stijgende of dalende trend.

Bepaling minimale teldekking 

Om te bepalen bij welke grenswaarden van aantal jaren, aantal plots, en aantal observaties de juiste trendklasse vastgesteld kan worden met TRIM, hebben we logistische regressie toegepast op de resultaten van de simulatie data. Hiervoor hebben we de simulatiedatasets met lage aantallen gebruikt, en zonder u-shape verdeling, om goed de laagste grenswaarden te kunnen bepalen. Verder is alleen een toenemende trend gebruikt om het aantal simulaties te beperken. 

De relatie tussen de binaire uitkomstvariabele ‘gelijk’ = 1 (voorspelde trendklasse met TRIM komt overeen berekende trendklasse) en de predictorvariabelen obs (aantal observaties), jaren (aantal jaren) en plots (aantal plots), evenals hun interacties, werd gemodelleerd met behulp van een generalized lineair model (GLM) met een binomiale verdeling. Het model kan als volgt worden uitgedrukt:

\begin{equation}
\text{logit}\big(P(\text{gelijk}=1)\big)
=
\beta_0
+ \beta_1\,\text{obs}
+ \beta_2\,\text{jaren}
+ \beta_3\,\text{plots}
+ \beta_4\,(\text{obs}\times\text{jaren})
+ \beta_5\,(\text{obs}\times\text{plots})
\end{equation}

Waarbij β-waarden in het model de geschatte effecten vertegenwoordigen van de verschillende predictorvariabelen.

Dit model is vervolgens gebruikt om het minimale aantal jaren, plots en observaties te voorspellen dat nodig is om met 95% zekerheid de juiste trendklasse weer te geven (Figuur 3).

2.3 Voorbeeld van het effect van skewness op het correct bepalen van de trendklasse in gesimuleerde datasets, waarbij rood aangeeft of de berekende trendklasse uit TRIM overeenkomt met de meegegeven trendklasse. Hierin is zowel links, rechts, midden, als u-shape skewness te zien (verschillende panelen) van het aantal observaties over het aantal jaren (a) en het aantal plots (b) met lage aantallen en het aantal observaties over het aantal jaren (c) en het aantal plots (d) met hoge aantallen.