9 Bilag:

Luftforureningen og luftvejseffekter hos fynske børn

9.1 Statistisk model

Antag, at man til hvert tidspunkt t h<r målt en række covariater x_1,t,…,x_k,t samt en respons y_t. I det foreliggende tilfælde, er y_t antallet af indlæggelser og covariaterne x_i,t er f.eks NO_X koncentration, SO₂ koncentration, temperatur, luftfugtighed, antal græs- eller birkepollen samt tid. Den anvendte model specificerer, at

y_t | x_i,s,i = 1,…,ks ≤ t ~ po(λ_t)

er uafhængige med:

log(λ_t) = ƒ₁₀(x_1,t) + … + ƒ_ls(x_1,t-s) + ƒ₂₀(x_2,t) + … + ƒ_ks(x_k,t-s)

hvor ƒ_it'erne er funktioner.

Modellen ovenfor er ikke fuldstændig specificeret, f.eks kan man højst bestemme funktionerne ƒ_it op til addition af en konstant. Typisk reparerer man dette, ved at indføre et konstantled samt kræve, at funktionerne summer til 0 over de observerede værdier af covariaterne, dvs. modellen bliver

log(λ_t) = α + ƒ₁₀(x_1,t) + … + ƒ_ls(x_1,t-s) + ƒ₂₀(x_2,t) + … + ƒ_ks(x_k,t-s)

hvor nu ∑ƒ_it = 0 for alle i.

Heller ikke denne restriktion er tilstrækkelig. I konkrete situationer vil man ofte være i stand til at fitte data perfekt, simpelthen fordi klassen af alle funktioner er for stor. Dette kan man så søge at reparere på forskellige måder. En mulighed er at reducere de mulige funktioner ƒ_it, ved at kræve at de ligger i et passende endelig dimensionalt rum. Traditionelle valg inkluderer trinvis konstante funktioner, polynomier og regressionssplines. En sådan model vil man sædvanligvis fitte ved at minimalisere den negative log-likelihood, hvor så ƒ_ij'erne ligger i det valgte rum.

En mere fleksibel model opnås, hvis man istedet blot kræver at ƒ_it'erne skal være glatte funktioner, men dette leder igen tilbage til problemet med en for stor funktionklasse. En smart mulighed er her at begrænse de mest ekstreme funktioner i denne klasse, f.eks kan man straffe de mindre glatte funktioner. En måde at gøre dette er at indføre et strafled J(ƒ_ij) til likelihood'en. Her skal J(ƒ_ij) være valgt, så mere glatte funktioner giver mindre værdier af J(ƒ_ij) end mindre glatte funktioner. Til et sådant valgt par, en familie af funktioner og et strafled, skal man så ikke længere minimalisere likelihooden, men istedet

l(ƒ_ij , y_t , x_lt , … , x_kt) + J(ƒ_ij)

Pointen er, at man nu kan tillade en uendelig dimensional klasse af funktioner, hvilket giver mulighed for en mere nøjagtigt modellering.

I det forhåndenværende tilfælde, modelleres counts'ne y_t som uafhængige poisson observationer, så udtrykket bliver:

y_tlog(λ_t(x_1,t , … , x_k,t-s)) - λ_t(x_1,t , … , x_k,t-s) - log(y_t!) + J(ƒ_1,0 , … , ƒ_k,s)

9.1.1 Udglatningsparametre

Der er mange forslag til nøjagtigt hvilke funktionsklasser der skal anvendes, f.eks polynomier, polynomielle splines der er stykvis polynomielle funktioner begrænset af kontinuitets- og glathedsbetingelser i delepunkterne samt tensor splines til flerdimensionale kovariater. Mere komplicerede konstruktioner som thin-plate splines og tensor thin-plate splines er også anvendt.

I de foreliggende analyser er der anvendt såkaldt cykliske splines til

modellering af sæsoneffekten og kubiske splines til øvrige ikke

parametriske effekter.

Selv efter fastlæggelse af funktionklassen resterer det vigtigste valg, nemlig bestemmelsen af strafleddet. De selv samme funktioner, kan tilknyttes forskellige rest led, hvilket igen har direkte indflydelse på inferencen. Medens man selvfølgelig kan forestille sig vilkårligt indviklede straf led, er opbygningen dog oftest mere simpel:

- Til hvert led ƒ_ij i=1,…k j=0,…s vælges en straf J_ij

- En samlet straf konstrueres som J(ƒ_ij) = ∑θ_ijJ_ij

Typisk er de enkelte straf led valgt, så de dobbelt afledte af funktionerne ikke bliver for store, fx kunne man have J_ij = ∫(ƒ”(x))²dx

De nye parametre θ_ij tjener til at afveje de enkelte led mod hinanden, men selv denne simplificering efterlader problemer i modelspecifikationen, for hvordan vælges strafparametrene θ_ij. Med mere end bare et par funktioner bliver ad hoc metoder svære at håndtere. Der er behov for en datadrevet procedure. En mulighed er, at forsøge en minimering af en størrelse der på passende måde måler afstanden fra den fittede model til den rigtige, fx den symmetriske Kullback-Leibler afstand

1/n ∑(λ_t(θ) - λ_t) (log(λ_t(θ)) - log(λ_t))

hvor λ_t(θ) er det til θ_ij'erne hørende estimat til tid t og λ_t er den sande værdi til tid t. Eftersom de sande λ_t'er jo er ukendte, må man forlade sig på at minimere et estimat af ovenstående størrelse. I [Wood, 2000] anvendes følgende algoritme:

- Den 'straffede' likelihood approksimeres med en første ordens Taylorudvikling omkring λ₀, hvilket erstatter leddet l(ƒ_ij , y_t , x_lt , … , x_kt) + J(ƒ_ij) med et kvadratisk led.

- For faste valg af strafparametrene θ_ij er minimeringen af denne approksimation et kvadratisk problem med kendt løsning, men spørgsmålet om hvilke strafparametre der skal anvendes resterer stadigvæk. Dette kan fx løses ved at minimalisere størrelsen

y˜^T(I - A)y˜ + 1/2tr(A)

hvor y˜_t = log(λ₀) + y/λ₀ - 1 og A er den af θ_ij'erne afhængende matrice der minimerer den kvadratiske sum ovenfor. Begrundelsen for dette valg er bla. at i normalfordelingstilfældet er dette led en estimator for den symmetriske Kullback Leibler afstand. I denne situation kan man også intuitivt fortolke de to led som estimatorer af henholdsvis varians og bias.

- Grundet den oprindelige approksimation af indsættes de nu opnåede estimater af λ_t'erne og θ_ij'erne i den kvadratiske tilnærmelse og det første skridt gentages.

Disse trin gentages nu indtil konvergens opnås. Der er indicier, men endnu intet bevis, for at de hermed fundne strafparametre approksimativt vil minimalisere den symmetriske Kullback Leibler afstand. Denne metode er implementeret i R-pakken mgcv og anvendt her.

Baysiansk tilgangsvinkel:

I den ovenstående tilgangsvinkel tillader man funktionerne ƒ_ij at variere i et meget stort rum. Fitningen bliver så muliggjort af, at man vælger at staffe meget volatile funktioner i dette rum. Man kan sige, at hele modellen arbejder ud fra en a priori antagelse om, at funktioner er mere sandsynlige jo glattere de er.

Denne tankegang minder til forveksling om den der ligger bag Baysiansk statistik: Man starter med sine nuværende antagelser, der kan være mere eller mindre godt funderet, og ser hvordan de forhåndenværende data ændrer disse antagelser. Mere præcist betyder dette, at man opstiller en model for responsen y_t og postulerer en (prior) fordeling på de ukendte parametre i denne. Fx kunne man vælge at opstille modellen

log(λ_t) = αNOX_t-1 + ƒ(influenza_t)

i hvilken parameteren α og funktionen ƒ( ) er ukendte og derfor skal tillægges fordelinger. Et muligt valg af fordelinger kunne være, at α var normaltfordelt med en middelværdi på 0.001 og en varians der selv følger en gammafordeling med stor varians. Disse valg afspejler i rimeligt omfang vores nuværende antagelser om parameteren α, nemlig at den forventede effekt er dårligt bestemt, men forventes at være omkring en promille. Funktionen ƒ( ) vil vi heller ikke tillade os at mene ret meget mere om end at den er nogenlunde glat.

Dette kan man formalisere ved at kræve at anden ordens differencerne

ƒ(influenza_t) - 2ƒ(influenza_t-1) + ƒ(influenza_t-2)

ikke bliver alt for store. En mulighed kunne her igen være, at disse differencer er normaltfordelte med middelværdi 0 og en varians der igen følger en gammafordeling, der selv har stor varians som ovenfor. Her skal man så undersøge hvor meget det præcise valg af denne sidste varians i gammafordelingen påvirker estimationen i den oprindelige model.

Denne model er fittet med sammenlignelige resultater for enkelte forureningskomponenter.

9.2 Omregningsfaktorer

9.2.1 TSP → PM10

Hvis man ser på det spring der skete ved overgangen mellem TSP og PM10 så faldt niveauerne med ca 1/3 - så en omregningsfaktor på 0.67 er det bedste estimat for niveauerne af PM10 før disse målinger startede. Hvis man ser på hvad der faktisk opsamles med TSP metoden så er der tale om et cut-off som varierer noget med vindhastigheden fra PM25 op til PM80 i ekstreme tilfælde. www2.dmu.dk/1_Viden/2_miljoe-tilstand/3_luft/4_maalinger/5_niveauer/6_Partikler/partikler_generelt.asp

9.2.2 Mg/m³ → ppm

Ozon: 1,98 mikrogram O₃/m³ per ppb

NO₂: 1,9 mikrogram NO₂/m³ per ppb

NO: 1,24 mikrogram NO/m³ per ppb

SO₂: 1,33 mikrogram SO₂/m³ per ppb

Benzen 3,5 mikrogram C₆H₆ per ppb