Notitie bepaling onzekerheidsmarge ammoniakemissie bij bedrijfswijziging
Opgesteld door:
Nico Ogink, onderzoeker Veehouderij en Milieu, Wageningen UR
Paul Goedhart, statisticus Biometris, Wageningen UR
Julio Mosquera, onderzoeker Veehouderij en Milieu, Wageningen UR
November 2022
Introductie en doelstelling
De minister van IenW heeft emissiefactoren van stalsystemen vastgesteld in de Regeling ammoniak
en veehouderij (Rav) om te beoordelen of de veehouderij voldoet aan de grenswaarde voor
huisvesting op grond van het Besluit emissiearme huisvesting. De Rav-emissiefactoren bieden
volgens het oordeel van de rechter onvoldoende zekerheid over de ammoniakemissie en depositie
van een individuele veehouderij om te kunnen hanteren bij de vergunningsverlening in het kader van
de Wet Natuurbescherming (Wnb). Een mogelijke oplossingsrichting zou kunnen zijn om de
onzekerheid van het verschil tussen de Rav-factoren in de oude en nieuwe bedrijfssituatie in deze
beoordeling te betrekken. Daarvoor is een onderbouwde methode nodig.
In deze notitie wordt de vraag beantwoord of er een methode beschikbaar is waarmee een
onzekerheidsmarge voor het verschil tussen emissiefactoren bij bedrijfswijziging kan worden
vastgesteld. Daarmee moet, met een zekere betrouwbaarheid, kunnen worden beoordeeld of een
bedrijfswijziging met overgang naar een ander stalsysteem in het hiervoor genoemde Wnb kader al
of niet tot toename van emissie leidt in een individueel geval. Voor de beantwoording van deze vraag
is een statistische analyse verricht op een dataset met ammoniakemissie meetgegevens van
melkveebedrijven. De dataset bestond uit bedrijven met conventionele en verschillende emissiearme
staluitvoeringen. Op basis van de analyseresultaten is een methode uitgewerkt voor het vaststellen
van een onzekerheidsmarge bij wijziging van stalsysteem op een individueel bedrijf. Om de
methodiek te illustreren zijn de onzekerheidsmarges berekend voor de emissiearme stalsystemen in
de dataset. De volledige technische uitwerking van de statistische analyse is in de annex van deze
notitie opgenomen. In overleg met de opdrachtgever zijn eveneens een aantal hieraan verbonden
deelvragen beantwoord.
Methode voor het berekenen van de onzekerheidsmarge van het verschil tussen emissiefactoren
bij bedrijfsaanpassing met een gewijzigd stalsysteem
De afgelopen 25 jaar zijn in de melkveehouderij nieuwe stalsystemen ontwikkeld om de
ammoniakemissie per dier te verminderen t.o.v. de conventionele ligboxenstal met roostervloer en
mestopslag onder de vloer. Het meest gehanteerde reductieprincipe is de emissiearme vloer die in
verschillende uitvoeringsvarianten beschikbaar is. Zowel van het conventionele stalsysteem als de
nieuwere emissiearme systemen zijn de jaargemiddelde ammoniakemissies per dierplaats
vastgesteld conform meetprotocol1. Daarbij wordt op minimaal vier bedrijven de ammoniakemissie
gemeten op een aantal dagen verspreid over een jaar. Deze emissiefactoren worden gebruikt voor
vergunningverlening in het kader van de Wet ammoniak en veehouderij en staan in bijlage 1 van de
Rav.

1

N. Ogink, J. Mosquera & A. Hol, 2017. Protocol voor meting van ammoniakemissie uit
huisvestingssystemen in de veehouderij. WLR rapport 1032.

De ervaring uit metingen aan individuele stallen van hetzelfde stalsysteem is dat de emissie in de tijd
varieert als gevolg van de weersomstandigheden (zoals temperatuur en windsnelheid). Daarnaast is
er variatie tussen bedrijven als gevolg van verschillen in bedrijfsmanagement (zoals bijvoorbeeld
voersamenstelling, management mestafvoer vloeren, stalbezetting uitgesplitst naar melkkoeien,
droogstaande koeien en niet-drachtig jongvee). Er is daardoor een aanzienlijke spreiding in emissie
tussen individuele bedrijven met hetzelfde systeem, maar ook tussen meettijdstippen op hetzelfde
bedrijf. Die variatie tussen en binnen bedrijven is in de analyse van emissiedata van
melkveebedrijven berekend als variatiecoëfficiënt (%cv)2. De geschatte %cv tussen melkveebedrijven
met hetzelfde stalsysteem in deze analyse is 21%. De variatiecoëfficiënt tussen individuele metingen
op hetzelfde bedrijf is 37%; deze is hoger dan die tussen bedrijven.
De onzekerheid van emissiefactoren kan in de vorm van betrouwbaarheidsintervallen worden
uitgedrukt. Dit interval geeft de onderste en bovenste waarde waarbinnen de geschatte waarde van
de emissiefactor met een gegeven betrouwbaarheid ligt. Meer gewenste zekerheid (hogere
betrouwbaarheid) leidt altijd tot langere (grotere) intervallen. Als we spreken over
betrouwbaarheidsintervallen is in de eerste plaats van belang aan te geven waar deze precies
betrekking op heeft. In de Rav staan emissiefactoren opgenomen die zijn gebaseerd op de gemiddeld
gemeten waarde in een steekproef van bedrijven. Het betrouwbaarheidsinterval van de gemiddelde
emissiefactor hangt af van het aantal bedrijven en het aantal meettijdstippen per bedrijf in de
steekproef en de spreiding tussen en binnen bedrijven. Met meer bedrijven in de steekproef wordt
het betrouwbaarheidsinterval voor de gemiddelde emissiefactor korter. Het
betrouwbaarheidsinterval van diezelfde gemiddelde emissiefactor toegepast op een individueel
bedrijf is altijd langer omdat hierin zowel de onzekerheid van de schatting van de gemiddelde
emissiefactor als de onzekerheid ten gevolge van variatie tussen bedrijven moet worden
verdisconteerd. In figuur 4 van de bijlage zijn beide typen betrouwbaarheidsintervallen opgenomen
voor de emissie van stalsystemen in de beschikbare dataset.
Wanneer een veehouderijbedrijf overschakelt naar een ander stalsysteem is het, in het kader van de
Wnb-vergunningverlening, van belang hoe zeker het is dat de stalemissie en de daaraan gekoppelde
depositie daadwerkelijk afneemt. Vanuit de opdracht is meegegeven dit voor een individuele
bedrijfssituatie te beschouwen. Voor het beoordelen of emissie bij overschakeling afneemt kan
gekeken worden naar het verschil tussen beide Rav-emissiefactoren (oud en nieuw) en de hieraan
verbonden onzekerheid in de vorm van een betrouwbaarheidsinterval. Het gaat hier dus om het
verschil in emissie per dierplaats. In de bijlage wordt uitgelegd dat op statistische gronden dit verschil
tussen twee stalsystemen het best kan worden uitgedrukt als de ratio tussen beide factoren (100 x
emissiefactor nieuw/emissiefactor oud). Bijvoorbeeld, bij overschakeling van een oude stal met een
emissiefactor van 13,0 kg NH3/jaar per dierplaats naar een nieuwe emissiearme stal met 6,5 kg
NH3/jaar per dierplaats is de ratio gelijk aan 50% (100 x 6,5/13,0). Hoe effectiever de werking van het
nieuwe stalsysteem op het betreffende bedrijf, des te lager deze ratio, en omgekeerd hoe hoger deze
ratio des te minder effectief is de emissiearme werking. Ratio’s boven de 100% geven aan dat de
bedrijfsemissie van de nieuwe stal hoger is dan van de oude stal.
In paragraaf 2.4 van de bijlage wordt toegelicht hoe het betrouwbaarheidsinterval van deze ratio op
individueel bedrijfsniveau kan worden berekend uit de variaties3 van beide geschatte Rav-factoren en
de variatie tussen individuele bedrijven. Een belangrijk uitgangspunt hierbij is dat de variatie tussen
bedrijven en binnen bedrijven gebaseerd is op een gepoolde schatting van alle stalsystemen. De
2

) Variatie binnen een reeks getallen wordt in de statistiek vaak uitgedrukt als standaarddeviatie, een maat
voor de afwijking van getallen van het gemiddelde. De variatiecoëfficiënt is een relatieve maat voor spreiding.
Het drukt de standaarddeviatie uit als percentage van het gemiddelde.
3
In de annex wordt hiervoor de term variantie gebruikt, een maat voor variatie gelijk aan de geschatte
standaardafwijking in het kwadraat.

bovengrens van het betrouwbaarheidsinterval voor de ratio tussen het nieuwe staltype en oude
staltype voor een individueel bedrijf kan dan gebruikt worden als criterium in de
vergunningsverlening voor een individueel bedrijf. Als deze bovengrens kleiner is dan 100% kan er
besloten worden dat er voldoende zekerheid is dat de NH3-emissie bij overgang op het nieuwe
staltype niet toeneemt bij een gelijk aantal dieren. Daarbij wordt verondersteld dat het gemeten
praktijkgemiddelde van een emissiearme stal niet door ander of inadequaat gebruik systematisch
afwijkt van de Rav-factor. Een belangrijk deel van het betrouwbaarheidsinterval wordt bepaald door
de variatie tussen individuele bedrijven. Deze tussenbedrijf variatie speelt een veel kleinere rol
wanneer in plaats van toetsing op individueel bedrijfsniveau getoetst zou worden op een gemiddeld
te verwachten emissieafname/toename van de populatie bedrijven met een emissiearm stalsysteem
t.o.v. het gemiddelde van de referentiebedrijven.
Ter illustratie van de methode zijn in de onderstaande figuur (figuur 5 uit de bijlage), voor de
stalsystemen in de dataset, de 90%-betrouwbaarheidsintervallen voor ratio’s ten opzichte van het
referentie staltype A1.100 weergegeven. In de figuur geven de rode lijnen de intervallen voor
toetsing op individueel bedrijfsniveau weer. Ter vergelijking staan ook de
betrouwbaarheidsintervallen weergegeven als er getoetst zou worden op populatie niveau (zwarte
lijnen). De waarden van de onder- en bovengrenzen van de rode lijnen staan in de tabel volgend op
de figuur weergegeven (afgeleid uit tabel 5 van de bijlage). In de tabel is te zien dat voor 10 van de 15
emissiearme stalsystemen sprake is van een bovengrens lager dan 100%. Bij overschakeling naar
deze systemen kan besloten worden dat er voldoende zekerheid is dat de emissie bij een gelijk aantal
dieren niet toeneemt. Ook zou in deze gevallen berekend kunnen worden met hoeveel dierplaatsen
een bedrijf dan nog maximaal kan uitbreiden; het aantal dierplaatsen zou dan maximaal met een
factor (100/bovengrens) mogen toenemen.
In de hierna weergegeven figuur is te zien dat de lengte van de betrouwbaarheidsintervallen
toeneemt naarmate de ratio’s hoger worden. Dit wordt veroorzaakt doordat de spreiding tussen
bedrijven proportioneel toeneemt met het emissieniveau. Daarnaast is er ook een doorwerking in de
betrouwbaarheidsintervallen van het aantal stallen in elk van de emissiearme datasets (3, 4, 5 of 6
stallen). De lengte van de intervallen op populatieniveau, de zwarte lijnen, neemt af bij metingen op
meer stallen; daarbij geldt als vuistregel dat vier keer zoveel stallen het interval tot de helft
reduceert. De lengte van het betrouwbaarheidsinterval op individueel bedrijfsniveau, de rode lijnen,
wordt echter voor een groot deel bepaald door de variatie tussen individuele bedrijven en deze
component in het interval neemt niet af bij meer stallen in de dataset van het betreffende
stalsysteem. Toetsing op individueel bedrijfsniveau levert dus een aanzienlijk ruimere bovengrens op
dan wanneer getoetst zou worden op de gemiddelde toename op populatie niveau.
In de getoonde figuur is staltype A1.100 (conventionele stal met roostervloer) als referentiestaltype
gebruikt. In principe kan de berekening ook voor andere referentiestaltype worden uitgevoerd. Als
voorbeeld is hier een 90% tweezijdig betrouwbaarheidsinterval gegeven. Voor het beoordelen van
toe of afname van emissie is alleen het al of niet overschrijden van de bovengrens van belang. Voor
deze eenzijdige toetsing halveert de overgebleven kans op overschrijding van de grenswaarde (de
helft van 10%) en geldt daarom een betrouwbaarheid van 95%. De berekening kan naar wens voor
elke andere betrouwbaarheid worden uitgevoerd.
Samengevat is het antwoord op de gestelde vraag dat er een kwantitatieve methode beschikbaar is
waarmee een onzekerheidsmarge voor het verschil tussen emissiefactoren bij bedrijfswijziging kan
worden vastgesteld. De methode is gebaseerd op het vaststellen van de ratio tussen de Rav-factoren
van het nieuwe en de oude stalsysteem in de betreffende bedrijfswijziging, en het berekenen van
een betrouwbaarheidsinterval van die ratio met een gewenst betrouwbaarheidspercentage. De
bovengrens van dit interval kan vervolgens als criterium worden gebruikt voor de beoordeling of er
sprake is van voldoende zekerheid dat de emissie niet toeneemt bij bedrijfswijziging. Bij een gelijk

beschikt kan worden over alle meetdata die ten grondslag liggen aan de berekende Ravemissiefactoren. Bij toepassing in de praktijk kan in de loop van de tijd nieuwe data van nieuwe
stalsystemen beschikbaar komen. Voor deze nieuwe stalsystemen kunnen dan
betrouwbaarheidsintervallen berekend worden op basis van de eerder vastgestelde gepoolde
variatie tussen bedrijven. De verwachting is namelijk dat het aantal bedrijven (n=81) in de nu
beschikbare dataset een voldoende stabiele schatting van deze variatie oplevert, en dat het
toevoegen van nieuwe informatie daar niet veel verandering in zal brengen.

Beantwoording van de deelvragen

1) Hoe kunnen stalemissies in individuele gevallen worden bepaald in het kader van
vergunningverlening Wnb?
Voor de bepaling van een stalemissie in het kader van vergunningverlening is het niet mogelijk vooraf
te beschikken over de gemeten emissie van die betreffende stal op het betreffende bedrijf omdat
deze meestal nog gebouwd moet worden en omdat praktijkstallen niet beschikken over eigen
meetinstallatie waarmee de ammoniakemissie wordt gemeten. Om toch een inschatting te maken
kan gebruik gemaakt worden van metingen uitgevoerd op andere bedrijven met een gelijke
stalinrichting. In de beantwoording van de hoofdvraag is aangegeven op welke wijze met
inachtneming van onzekerheidsmarges dit kan worden vorm gegeven.
2) In hoeverre kunnen Rav-factoren daarbij een rol spelen?
Daarbij dient aandacht te worden besteed aan de volgende oplossingsrichtingen/subvragen:
a) Hoe nauwkeurig kan de inherente onzekerheid van stalemissies in individuele gevallen worden
bepaald, waarbij onderscheid gemaakt wordt in:
- Nauwkeurigheid van de metingen
- Hogere of lagere emissies als gevolg van voer- en managementmaatregen die niet als voorschriften
in de vergunning zijn opgenomen.
De Rav-factor van een stalsysteem is het gemiddelde van de gemeten emissie op een aantal
bedrijven volgens een vastgelegd meetprotocol. Zoals hierboven toegelicht is de nauwkeurigheid van
de geschatte Rav-factor afhankelijk van de steekproefomvang. Daarnaast blijkt uit meetseries dat er
aanzienlijke variatie tussen bedrijven bestaat met hetzelfde stalsysteem. De nauwkeurigheid
waarmee een jaaremissie van een individueel bedrijf kan worden geschat hangt af van de
nauwkeurigheid van de geschatte Rav-factor en de variantie tussen bedrijven. In de statistische
analyse die voor deze notitie is uitgevoerd zijn de 90%-betrouwbaarheidsintervallen berekend voor
de Rav-factoren én voor individuele bedrijven met zo’n Rav-factor. Zie hiervoor tabel 4 en Figuur 4
van de bijlage.
In de afzonderlijke meetpunten van de dataset kan geen onderscheid gemaakt worden tussen de
bijdrage van de meetnauwkeurigheid en de spreiding door specifieke bedrijfsfactoren
(voersamenstelling, management) op het moment van meten. In het meetprotocol wordt wel vereist
dat gasconcentraties in duplo worden gemeten. Dat geeft informatie over de toevalscomponent van
de meetfout. De aangeleverde data van de meetrapporten zijn echter gebaseerd op duplogemiddelden. Verwacht mag worden dat het kwaliteitsniveau van de meetmethoden in het
meetprotocol bij adequate kalibratie-procedures overigens zodanig is dat meetnauwkeurigheid een
ondergeschikte rol speelt t.o.v. de spreiding door bedrijfsfactoren.

b) Kan met een wetenschappelijk onderbouwde generieke onzekerheidsfactor met de benodigde
zekerheid (geen redelijke wetenschappelijke twijfel) worden aangetoond dat er bij bedrijfswijzigingen
of -uitbreidingen met emissie-arme stalsystemen geen depositietoename plaatsvindt?
- Tot welke mate van (on)zekerheid van de emissie (en daarmee depositie) leidt toepassing van de
factor in specifieke gevallen?
- Welk percentage biedt, vanuit wetenschappelijk oogpunt, voldoende zekerheid?
In het eerste deel van deze notitie is aangegeven hoe vanuit beschikbare gegevens een 95%betrouwbaarheidsbovengrens voor de ratio tussen staltypes bij bedrijfsomschakeling kan worden
berekend (NB . Deze bovengrens verdisconteert zowel variatie tussen bedrijven als variatie tussen
meettijdstippen op hetzelfde bedrijf. In landbouwkundig onderzoek wordt bij statistische toetsing
vaak met een 95% betrouwbaarheid gewerkt, maar dat is een conventie. Even zo goed kan dit
berekend worden voor 99% betrouwbaarheid. Voor het vaststellen van het betrouwbaarheidsniveau
in de geschetste situatie van vergunningverlening is in feite aan de orde welk risico op toename van
emissie voor het betreffende bedrijf men acceptabel vindt. Dat is een maatschappelijke
belangenafweging die beleidsmatig moet worden gemaakt want 100% zekerheid bestaat niet. Wij
kunnen daarom niet aangeven welk betrouwbaarheidspercentage ‘voldoende zekerheid’ biedt.
c) Hoe verhoudt het onderzoek dat op dit moment wordt uitgevoerd door de WUR, over de spreiding
die optreedt door gebruik van emissie-arme stalsystemen in de praktijk, zich tot de hier onderzochte
inherente onzekerheid? (kwalitatief)
Recente meetdata van enkele emissiearme melkveestallen gemeten uit het onderzoek van Schep et
al. ( 2022, WLR-rapport 1299) levert een beeld met een spreiding van gelijke omvang als die
berekend in de analyse van Rav-metingen in de bijlage. Ook het proportionele karakter van spreiding
is hierin zichtbaar.
d) Hoe verhouden de correctiefactoren uit de NEMA-studie zich tot de onzekerheidsfactor ten
behoeven van toestemmingverlening voor de Wnb? (kwalitatief).
In de NEMA-systematiek worden de Rav-factoren gecorrigeerd voor de ontwikkelingen in het
landelijk gemiddelde TAN-excretieniveau van de verschillende diercategorieën t.o.v. het TAN-excretie
niveau tijdens het jaar van de meting van een Rav-factor. Het gaat om geringe procentuele correcties
per jaar die gelijk zijn voor alle stalsystemen binnen een diercategorie. Eventuele correctie voor
trends in TAN-excretie hebben naar verwachting een gering effect op de berekening van de
onzekerheidsfactor volgens de voorgestelde methode omdat de metingen in de dataset betrekking
hebben op een beperkt tijdvak.

Figuur 1: Dag-nummers voor elk bedrijf waarop de NH3 emissie is gemeten. Het dag-nummer van de
eerste meting is daarbij op 0 gezet.

Figuur 2: NH3-emissies op elk bedrijf waarbij alternerende kleuren zijn gebruikt om de bedrijven
makkelijker te kunnen onderscheiden. De horizontale strepen geven het gemiddelde per bedrijf.

2. Statistische analyse
Een modified Park test voor heteroscedasticiteit (1) geeft aan dat de variantie van de NH3-emissies
ongeveer evenredig is met het gemiddelde in het kwadraat. Om de variantie te stabiliseren is daarom
een statistische analyse op de log-getransformeerde emissies uitgevoerd.
In de statistische analyse moet een onderscheid gemaakt worden tussen (1) variatie tussen bedrijven
met hetzelfde staltype en (2) variatie tussen herhaalde metingen op hetzelfde bedrijf. De eerste
variatiebron hangt samen met verschillen in bedrijfsspecifieke kenmerken. De tweede variatiebron hangt
samen met weersomstandigheden (zoals temperatuur, windsnelheid, luchtvochtigheid) en bijvoorbeeld
de precieze stalbezetting op het meetmoment. Het gebruikte statistische model om de NH3 emissies te
analyseren is daarom als volgt
Log(NH3-emissie) = staltype + bedrijf + meting

(1)

Hierin representeert “staltype” een te schatten parameter 𝜇𝑅𝑎𝑣 voor elk staltype en zijn bedrijf en meting
random effecten met gemiddelde 0 en variantie respectievelijk 𝜎𝐵2 en 𝜎𝑀2 . De interpretatie van de variantie
componenten op de oorspronkelijke niet-getransformeerde schaal is in termen van variatiecoëfficiënten.
De variatiecoëfficiënt tussen bedrijfsgemiddelden met dezelfde Rav code is gelijk aan √exp(𝜎𝐵2 ) − 1. De
variatiecoëfficiënt tussen metingen binnen een bedrijf is gelijk aan √exp(𝜎𝑀2 ) − 1. De parameters van dit
model worden geschat met REML.

a.

Betrouwbaarheidsintervallen voor populatie medianen

De statistische analyse resulteert in schattingen voor de populatie gemiddelden 𝜇𝑅𝑎𝑣 voor de Rav factoren
op de log schaal en bijbehorende standaardafwijkingen. Schattingen voor populatie medianen op de
oorspronkelijke NH3 emissie schaal worden dan gegeven door exp(𝜇̂ 𝑅𝑎𝑣 ) waarin exp(∙) de exponentiële
functie representeert. Op dezelfde wijze kan een betrouwbaarheidsinterval voor een Rav factor op de logschaal terug-getransformeerd worden naar de oorspronkelijke schaal. Tevens kunnen verschillen tussen
geschatte Rav factoren op de log-schaal getransformeerd worden naar ratio’s van Rav factoren op de
oorspronkelijke emissie schaal. Ook voor deze ratio’s kunnen betrouwbaarheidsintervallen worden
berekend. Daarbij zijn vooral ratio’s ten opzichte van het referentie staltype A1.100 interessant.

b.

Predictie-intervallen voor een individuele bedrijf

Voor een individuele bedrijf met een bepaald staltype kan een predictie interval voor de NH3 emissie
worden berekend. Daarbij moet rekening gehouden worden met (1) de onzekerheid waarmee het
populatie gemiddelde 𝜇𝑅𝑎𝑣 op de log-schaal wordt geschat en (2) de inherente variantie 𝜎𝐵2 tussen
bedrijven. De onzekerheid in de schatting van het populatie gemiddelde 𝜇𝑅𝑎𝑣 kan gereduceerd worden
door gegevens voor meer bedrijven te verzamelen. De inherentie variatie tussen bedrijven 𝜎𝐵2 kan niet
gereduceerd worden; deze wordt hooguit nauwkeuriger geschat als er meer gegevens beschikbaar zijn.
2
Als 𝜎𝑅𝑎𝑣
de variantie is van de schatting voor 𝜇𝑅𝑎𝑣 , beiden op de log-schaal, dan is de variantie van de
2
schatting van de Rav factor voor een individueel bedrijf gelijk aan 𝜎𝑅𝑎𝑣
+ 𝜎𝐵2 . Op basis hiervan kan, op de

gebruikelijke wijze, een predictie interval worden berekend voor de Rav factor van een individueel
bedrijf. Voor een verschil tussen twee Rav factoren op een individueel bedrijf geldt een variantie van
2
2(𝜎𝑅𝑎𝑣
+ 𝜎𝐵2 ). De ondergrens voor de variantie op individueel bedrijfsniveau wordt bepaalt door 𝜎𝐵2 en deze
wordt niet gereduceerd indien meer gegevens beschikbaar zijn.

c.

Vergunningsverlening

De bovengrens van het predictie-interval voor de ratio tussen een nieuw staltype en het referentie
staltype A1.100 voor een individueel bedrijf kan gebruikt worden in de vergunningsverlening voor een
individueel bedrijf. Indien deze bovengrens kleiner is dan 100% dan bestaat er “voldoende” zekerheid
dat de NH3-emissie bij overgang op het nieuwe staltype niet toeneemt (bij een gelijke bezetting). De
mate van zekerheid wordt gedefinieerd door de betrouwbaarheid, bijvoorbeeld 95% of 99%, van het

predictie interval. Dit is een beleidskeuze. Hier wordt gerekend met een (eenzijdige) betrouwbaarheid
van 95% voor de bovengrens en dus een betrouwbaarheid van 90% voor de tweezijdige intervallen.

d.

Technische uitleg voor het berekenen van predictie-intervallen

Voor de berekening van een predictie interval voor een individueel bedrijf wordt gebruik gemaakt van het
model Log(NH3-emissie) = Rav-Factor + Bedrijf. Hierin representeert Rav-Factor de schatter voor de
2
gemiddelde Rav factor 𝜇𝑅𝑎𝑣 met variantie 𝜎𝑅𝑎𝑣
en Bedrijf een random bedrijfseffect met gemiddelde 0 en
2
variantie 𝜎𝐵2. Het gemiddelde voor een individueel bedrijf is dan 𝜇𝑅𝑎𝑣 met variantie 𝜎𝑅𝑎𝑣
+ 𝜎𝐵2. Een 90%
2
predictie interval op de logschaal wordt dan gegeven door 𝜇𝑅𝑎𝑣 ± 𝑡 √𝜎𝑅𝑎𝑣
+ 𝜎𝐵2 waarin 𝑡 het tweezijdige

90% percentage punt is van de t-verdeling met aantal vrijheidsgraden behorende bij de schatter voor 𝜎𝐵2.
Op de oorspronkelijke schaal wordt het interval dan gegeven door de e-macht van de eindpunten. Met
data voor 81 bedrijven verdeeld over 16 verschillende staltypen is uitgegaan van 65 (=81-16) graden
van vrijheid voor de t-verdeling. Er wordt hier gesproken over een predictie-interval in plaats van een
betrouwbaarheidsinterval omdat het een voorspelling betreft voor een individueel bedrijf in plaats van
een interval voor een parameter zoals 𝜇𝑅𝑎𝑣 .
Een verschil tussen twee Rav factoren voor een individueel bedrijf wordt gegeven, op log-schaal, door
(Rav-Factor + Bedrijf)1 – (Rav-Factor + Bedrijf)2 waarin het subscript de twee verschillende Rav factoren
2
aangeven. Het gemiddelde van die expressie is 𝜇𝑅𝑎𝑣,1 − 𝜇𝑅𝑎𝑣,2 en de variantie is gelijk aan 2 (𝜎𝑅𝑎𝑣
+ 𝜎𝐵2 ).
Hiermee kan een predictie interval voor de ratio voor een individueel bedrijf op de oorspronkelijke schaal
worden berekend. Het random effect Bedrijf in het verschil valt niet weg omdat een relatief hoge (of
lage) emissie van het oude staltype op een individuele bedrijf niet noodzakelijkerwijs gevolgd wordt door
een relatief hoge (of lage) emissie van het nieuwe staltype op datzelfde bedrijf. Het is immers
grotendeels onbekend wat de oorzaak is van de variatie tussen metingen. Wellicht wordt een deel van de
tussen bedrijven variatie verklaard door bijvoorbeeld het rantsoen. Bij een gelijkblijvend rantsoen zou
dan gerekend kunnen worden met een kleinere 𝜎𝐵2 en dus kortere intervallen. Welk deel van 𝜎𝐵2
aangerekend kan worden aan het rantsoen kan echter niet gekwantificeerd worden. De afgeleide
intervallen voor de ratio voor een individueel bedrijf zijn daardoor mogelijk conservatief.

3. Resultaten
De waargenomen NH3-emissies zijn conform bovenstaande model (1) geanalyseerd. Figuur 3 geeft aan
dat de spreiding van de residuen van het aangepaste model min of meer onafhankelijk is van de gefitte
waarde waarmee aan de belangrijke voorwaarde van homogeniteit van variatie is voldaan.

Figuur 3: Residuen tegen gefitte waarden voor het statistische model op log-getransformeerde
NH3-emissies. De rode punten zijn voor het referentie staltype A1.100.

1.

De modified Park test

Een algemene klasse van variantiefuncties voor waarnemingen wordt gegeven door 𝑉𝑎𝑟 = 𝜎 2 𝜇𝑃 waarin
𝑉𝑎𝑟 de variantie van de metingen per bedrijf, 𝜇 het gemiddelde van een bedrijf, 𝑃 een positieve macht en
𝜎 2 een proportionaliteit constante. 𝑃 = 0 correspondeert met de normale verdeling, 𝑃 = 1 met de quasiPoisson verdeling en 𝑃 = 2 met de gamma verdeling die goed vergelijkbaar is met de lognormale
verdeling. De variantiefunctie kan geschreven worden als het log-lineaire regressiemodel log(𝑉𝑎𝑟) =
log(𝜎 2 ) + 𝑃 log(𝜇). Een gegeneraliseerd lineair model (GLM) met de gamma verdeling en een log-link kan
dan gebruikt worden om de macht 𝑃 te schatten uit de waargenomen varianties en de log-gemiddelden.
Daarbij wordt de gebruikelijke veronderstelling gemaakt dat de varianties een gamma verdeling volgen.
Dat resulteert in een schatting 2.39 voor 𝑃 met standaardafwijking 0.23. Het is dus redelijk om een
waarde 𝑃 = 2 aan te nemen. Dat impliceert een gamma verdeling voor de emissie waarnemingen, of min
of meer equivalent een lognormale verdeling.
Manning, W. G., & Mullahy, J. (2001). Estimating log models: To transform or not to transform? Journal
of Health Economics, 20, 461–494. https://doi.org/10.1016/S0167-6296(01)00086-8
Park, R. E. (1966). Estimation with heteroscedastic error terms. Econometrica, 34, 888.
https://doi.org/10.2307/1910108