Kisterületi becslések: rövid áttekintés a korszerű módszerekről

Mihályffy László,

a Központi Statisztikai Hivatal ny. statisztikai főtanácsadója E-mail: laszlo.mihalyffy@ksh.hu

Kisterületi becslések:

rövid áttekintés a korszerû módszerekrôl

Ennek a dolgozatnak az a célja, hogy használható információt nyújtson azoknak a statisztikusoknak, szakstatisztikusoknak és módszertanosoknak egyaránt, akik mun- kájuk során első ízben találkoznak kisterületi becslési feladatokkal. Merész vállalko- zás – mondhatja az Olvasó –, ha figyelembe vesszük, hogy Rao ([2003], [2015]) munkája két kötetben mintegy 700 oldalon szinte száz százalékban lefedi a tudo- mányág eddigi eredményeit.

1. A kisterületi becslések fogalma

A mintavételi eljárások területén a kisterületi becslések iránti igényt a kezdetek- től, azaz a múlt század hetvenes éveitől kezdve mind a mai napig a következő kö- rülmény határozza meg. Adva van egy véges sokaság és ennek egy, a reprezentatív megfigyelés követelményeinek megfelelő mintája, ide értve azt is, hogy a minta a teljes véges sokaságra nézve kellő pontosságú becsléseket eredményez. Gyakran előfordul az a jelenség, hogy a felhasználókat a teljes véges sokaság egyes részeire vonatkozó becsült adatok is érdeklik, azonban a tekintett részsokaságokhoz tartozó részminták közül több vagy kevesebb nem eredményez kellő pontosságú becslést, mivel a mintanagyság nem elég nagy. Példák erre a hazai társadalomstatisztikai min- ták, amelyek országos szinten szükségszerűen teljesítik a pontossági követelménye- ket, megyei szinten a nagyobb egységek (Budapest és Pest megye, esetleg Borsod- Abaúj-Zemplén) kivételével azonban nem. Ebben a kontextusban tehát a megyék többsége kisterület.

Az említett példákhoz hasonlóan a kisterületek vizsgálatára is kiterjedő felvételek számottevő részében a terület nagysága a földrajzi terület nagyságával kapcsolatos,

ez azonban nincs minden esetben így. Ha a vizsgált véges sokaság a teljes népesség vagy annak 18 éveseket és azoknál idősebbeket magában foglaló része, akkor kiste- rületet képezhetnek például a mélyszegénységben élők, vagy valamilyen krónikus betegségben szenvedők, tehát a társadalom olyan kisebb létszámú csoportjai, ame- lyekből egy szokásos méretű¹ reprezentatív mintába csak nagyon kevés egység (sze- mély vagy háztartás) kerül.

2. A kisterületi becslések módszerei, feladattípusok

Ez a fejezet néhány, a kisterületi becslések szempontjából fontos fogalom (külső információ, implicit és explicit modellek) tisztázása mellett rövid utalást tartalmaz a bayesi statisztika alapjaira, mivel a kisterületi becslési módszerek egy része ezekre épül. A bayesi módszerek alkalmazásának szempontjából nélkülözhetetlen numeri- kus integrálási eljárások (Markov-lánc Monte-Carlo-módszer stb.) ismertetése terje- delmi okok nem volt lehetséges.

2.1. Hogyan javíthatjuk a becsült kisterületi adatok pontosságát?

Minden kisterületi becslési feladat kiinduló pontja és egyúttal nélkülözhetetlen kelléke egy véges sokaságból kiválasztott reprezentatív minta. A tervezett felvétel szempontjából a teljes sokaságnak azok a részei lesznek kisterületek, amelyekből az adott minta nem tartalmaz kellő pontosságú adat becsléséhez elegendő számú megfi- gyelést. Sem a mintavételi tervet, sem a mintát nem változtathatjuk meg, ezért

– kisterületi adataink közvetlen (azaz, mintából származó) becslé- sének javításához külső információra van szükségünk, és mivel

– adataink javítás, módosítás eredményei, kisterületi becsléseink szerencsés kivételektől eltekintve torzítottak lesznek.

A külső információt, pontosabban ennek forrásait a következőképpen értelmez- zük. Tekintsünk egy ismérvet, például 18 éves vagy idősebb, keringési rendszer be- tegségében szenvedő személy, és legyen S az országos minta, s pedig egy (ebben az

1 Mint ismeretes, a minta nagysága nem a (véges) sokaság nagyságától, elemszámától függ, hanem a be- csült adatok pontosságára vonatkozó követelményektől. A munkanélküliségi ráta adott pontosságú becsléséhez egy kis országban (például Ausztriában) körülbelül ugyanakkora mintára van szükség, mint egy nagy országban (például Németországban).

esetben földrajzi értelemben is) kisterület részmintája. A kisterület szempontjából a külső információ alábbi forrásai jöhetnek szóba.

– az S teljes minta, illetve az abból becsült adat;

– egy az S mintától különböző S’ minta, amely nem tartalmazza az említett ismérvre vonatkozó megfigyeléseket, megfigyelt ismérveiből azonban olyan segédváltozókat lehet meghatározni, hogy azok és a ke- resett adat között (18 éves vagy idősebb személy keringési rendszer betegségével) szoros korreláció áll fenn. Ebben az esetben a segédvál- tozókat az s minta alapján is elő kell tudni állítani;

– a keresett kisterületi adatnak egy korábbi felmérésből származó, nagy pontosságú értéke.

A teljes minta, amely részhalmazként tartalmazza a kisterületek mintáit, a külső információ egyik legegyszerűbb forrása. A külső információ forrásának másik cso- portját alkotják a népszámlálások és az egyéb cenzusok adatállományai, továbbá a regiszterek.

2.2. A kisterületi becslési eljárások modelljei és módszerei

Minden kisterületi becslési eljárás egy modellen alapul, még akkor is, ha a feladat egyszerűsége miatt a modell jelenléte nem nyilvánvaló – ilyenkor implicit modellről beszélünk. Időrendben főként a korábbi eljárások tartoznak ehhez a csoporthoz. A későbbi és egyben korszerűbb eljárások az explicit modellekre épülnek, ezek szinte kivétel nélkül regressziós modellek. Az explicit modellek további két nagy csoportja a lineáris és a nem lineáris, az utóbbiaknak több fajtája van: bináris-, béta-bináris-, Poisson-gamma-, logit-normál- stb. modellek. Ezek lehetnek egy- vagy több lépcső- sek²; a 3. fejezetben bemutatott explicit modellek valamennyien két lépcsősek, az első lépcső a mintavételi hiba eloszlásával kapcsolatos.

Becslési eljárás, tehát a számítási algoritmus szerint az implicit modellre épülő el- járások többségénél – de nem mindegyiknél – kizárólag a mintavételi terv(ek) által meghatározott becslési módszerekre van szükség, vagyis az értékösszegek, átlagok, és arányok becslésére. Az explicit modellek alkalmazására épülő módszereket a célul kitűzött becslés jellemzőire utaló angol nyelvű megnevezés rövidítésével jelöljük:

– BLUE (best linear unbiased estimator – leghatásosabb lineáris torzítatlan becslő függvény),

2 Ebben a vonatkozásban a „lépcső” kifejezés mellett a „szint” is használatos.

– BLUP (best linear unbiased prediction – leghatásosabb lineáris torzítatlan előrejelzés)³,

– EBLUP (empirical best linear unbiased prediction – empirikus leghatásosabb lineáris torzítatlan előrejelzés),

– HB (hierarchikus Bayesi becslés), – EB (empirikus Bayesi becslés).

Mivel gyakorlatilag minden kisterületi becslés torzított, nyilvánvaló az a köve- telmény, hogy a keresett becslésnek mind a szórásnégyzete, mind pedig a torzítása lehetőség szerint kismértékű legyen. Ezt a célt az MSE (mean squared error – átlagos négyzetes hiba) minimalizálására való törekvés fejezi ki, erre utalnak a BLUE, BLUP, illetve EBLUP elnevezések is. (Gyakorlati megfontolásból alanyesetben az MSE rövidítést használjuk, ragozott formában azonban célszerűbbnek látszik a ma- gyar nyelvű kifejezés alkalmazása.) A lehetőség szerinti minél kisebb MSE elérése természetesen nemcsak az explicit modelleknél cél, hanem az implicit modellek esetén is.

Fontos felhívni a figyelmet a következő körülményre. Az implicit modelleken alapuló kisterületi eljárások többsége az adott mintavételi terv szerint véges sokasá- gokon működik, ezért a szórásnégyzetre, valamint a torzításra a mintavételi tervnek megfelelő becslést szolgáltat. Az explicit modellekre épülő eljárások ezzel szemben a modellt meghatározó eloszlás szórásnégyzetének megfelelő pontossági paramétert rendelik a becsült mutatóhoz, ami általában eltér a mintavételi terv szerinti szórás- négyzettől. Vannak azonban olyan modellek is, amelyek kezelik ezt a problémát.

1. táblázat A modelltípusok és a becslési eljárások kapcsolata

A modell típusa Használható becslési módszerek

Implicit modell Becslés a mintavételi terv alapján Bizonyos esetekben kalibrálás Általánosított lineáris modell Standard regressziós technika

Hierarchikus Bayes-módszer Empirikus Bayes-módszer Nem lineáris modellek Hierarchikus Bayes-módszer

Empirikus Bayes-módszer

3 A BLUP- és az EBLUP-becsléseknél a (lineáris) regresszió paramétereinek a becslése „BLUE”, a para- méterek felhasználásával meghatározott becslések, előrejelzések tulajdonsága „BLUP”. Az EBLUP csak köze- lítően torzítatlan, lásd a 3.4. alfejezetet.

A kisterületi becslési módszerek többségét eredetileg általánosabb célokra dol- gozták ki, a kisterületi feladatokra való alkalmazásukhoz esetenként kisebb-nagyobb módosításra volt szükség. Kivétel ez alól az implicit modellekhez tartozó módszerek közül néhány, amelyeket kizárólag csak kisterületi becslési feladatokban alkalmaz- nak; ezek közül kettőt ismertetünk a 3. fejezetben, nevezetesen a szintetikus és a kombinált becslést. Feltételezhető, hogy a táblázatban felsorolt becslési módszerek a kétféle Bayes-módszer kivételével széles körben ismertek, ezért a következő alfeje- zetben csupán erről a kettőről lesz szó.

2.3. A hierarchikus Bayes-módszer

Kiinduló pontunk a következő feladat: adott egy véges sokaságból származó min- tán megfigyelt n számú megfigyelés x  (x₁, , ..., ) ,x₂ x_n ^T becsüljük ezek alapján a sokaság valamely θ paraméterét.⁴ A mintavételi eljáráson alapuló hagyományos eljárás szerint a megfelelő képlet és az x₁, , ..., x₂ x_n megfigyelések segítségével kiszámítjuk a paraméter ˆθ becslését, és lehetőség szerint a mintavételi hibáját is valamilyen valószínűségi szinten.

A bayesi szemlélet szerint θ nem konstans érték, hanem valószínűségi változó, amelynek eloszlását kell becsülnünk. A következő lépéseket kell tennünk.

1. Rögzítünk egy f θ₁( ) sűrűségfüggvényű a priori eloszlást, ami lehet korábbi tapasztalatokon alapuló informatív adat, de lehet például a minden információt nélkülöző „lapos”, „nem igazi” megoszlás

1( )

f θ  1 sűrűségfüggvénnyel.

2. Meghatározzuk az f₂( | )x θ feltételes sűrűségfüggvénnyel jel- lemzett likelihood eloszlást, amely az összes olyan θ eloszlást képvi- seli, amelyek mellett az x megfigyelés létrejöhet. Az explicit modelle- ken alapuló kisterületi becslési eljárásoknál a likelihood függvény adottság.

3. Bayes tételének értelmében előállítjuk az

₃ ^{2 1}

0

( | ) ( ) ( | )

( ) f θ f θ f θ  xf

x x /1/

4 A dolgozatban a matematikai képletek a MathType egyenletszerkesztő jelölési konvenciójának felelnek meg, a félkövér betűk a vektort vagy mátrixot, a skalárokat dőlt latin vagy görög betűk jelölik. A felső T index transzponálást jelent. A paraméter lehet skalár vagy vektor, az utóbbi esetben x komponensei is vektorok.

sűrűségfüggvényű posterior eloszlást, amely a hagyományos eljárással meghatározott ˆθ paramétert helyettesíti.

A posterior eloszlás ˆθ^B várható értéke a hagyományos eljárással adódó ˆθ érték Bayes becsléssel meghatározott megfelelője, az /1/ képlet figyelembevételével a következő képlettel számíthatjuk ki:

θˆ^B



θ f θ₃( | )x dθ  c θ f



₂( | ) ( )x θ f θ dθ₁ , /2/

ahol c 1 f₀( )x és f x₀( )a mintabeli megfigyelések vektorának feltétel nélküli sűrűségfüggvénye az x helyen. Ha θ nem skalár mennyiség, akkor a /2/ jobb oldalán szereplő kifejezés sokdimenziós integrál lehet, és mivel gyakorlati alkalmazásokban az integrandust csak ritkán lehet zárt alakban előállítani, ˆθ^Bvalamint a θ feltételes szórásnégyzetének kiszámítása rendkívül számításigényes feladat. A Markov-lánc Monte-Carlo-módszerek megjelenése előtt a numerikus integrálás akkori metodiká- jával nem mindig lehetett a /2/ jobb oldalán szereplő kifejezést számszerűsíteni, és ez a körülmény akkortájt gátat szabott a bayesi módszerek elterjedésének. Napjainkban már nem ez a helyzet, bár a szóban forgó integrálok kiszámítása még a modern eljá- rások mellett is számításigényes.

2.4. Az empirikus Bayes-módszer

Ez az eljárás abban különbözik a hierarchikus Bayes-módszertől, hogy a prior de- finiálásánál felhasználják a x megfigyeléseket vagy azok egy részét, miáltal a szóban forgó eloszlás elveszíti „prior” jellegét. Ez kétségtelenül azzal az előnnyel jár, hogy a /2/ relációban szereplő integrálokkal kapcsolatos számítástechnikai problémák jelen- tős mértékben kisebbek lesznek, ennek azonban az az ára, hogy az eredeti Bayes- módszer filozófiája sérül, lemondunk a prior eloszlás által képviselt véletlenszerű ingadozások vizsgálatáról. Mindenesetre az empirikus Bayes-módszert szokás az

„igazi” (tehát hierarchikus) Bayes-módszer közelítésének tekinteni.

Az „igazi” prior használatától való eltérésnek nagyjából két útja van: egyrészt a prior által hordozott információt az /1/ összefüggésben megfogalmazott szabálytól eltérő módon építjük be a likelihood függvénybe, másrészt ún. konjugált priort hasz- nálunk. Ez azt jelenti, hogy a választható prior eloszlások körét úgy szűkítjük, hogy a becsülni kívánt paramétereknek ugyanolyan típusú matematikai kifejezésben kell szerepelniük a prior sűrűségfüggvényében, mint a likelihood eloszlás sűrűségfügg- vényében.

3. Példák kisterületi modellekre, becslésekre

Az alfejezetek tartalma a következő: 3.1.–3.3. szintetikus becslés, kombinált becslés, struktúrát megőrző becslés; 3.4. az EBLUP-eljárás, kisterületi adatok becslé- se általános vegyes lineáris regressziós modellel; 3.5.–3.8. példák az empirikus bayesi módszer alkalmazására. Az összes példa konjugált prior alkalmazásán alapul.

Más típusú empirikus bayesi módszer, illetve a hierarchikus Bayes-módszer bemuta- tására elsősorban terjedelmi szempontok miatt nem került sor.

3.1. Szintetikus becslés

Feltesszük, hogy adott egy mintavételi terv, továbbá az U véges sokaságból a tervnek megfelelően kiválasztott s₁ és s₂ minta, melyek elemszáma n₁ illetven₂, n1<< n₂.

s

₁ egy kisterület mintája, s₂ egy nagyobb tartományé, és egybe is eshet a teljes mintával. A sokaság i eleméhez az x y z_i, , ,_{i i} stb. értékek tartoznak. X és X az xi ismérvhez tartozó sokasági átlag és értékösszeg, x, ˆX, illetve x, ˆX pedig X és X becslése az s₁, illetve az s₂ minta alapján. Az y_i, z_i stb. ismérvekhez tar- tozó sokasági értékeket és a mintákból származó becsléseiket analóg módon értel- mezzük. Az U sokaság N elemszámának a mintákból származó becsléseit ˆ -vel,N

illetve ˆ -velN jelöljük.

A szintetikus becslés azon a feltevésen alapul, hogy egy értékváltozónak az átlaga egy kisterületen és egy nagyobb területen megegyezik. Ha ebben az esetben a tekin- tett átlag becslése a kisterületen a minta csekély elemszáma miatt nem kellő pontos- ságú, akkor a kisterülethez tartozó becslést a nagyobb területhez tartozó megfelelő becsléssel helyettesítjük, feltéve természetesen, hogy az utóbbi pontossága megfele- lő. Ha értékváltozó helyett kategorikus változóval van dolgunk, akkor az átlag helyé- be az arány lép. (Lásd a 2. táblázatot.)

2. táblázat A szintetikus becslés alapesetei

Mutató Közvetlen becslés Szintetikus becslés Megjegyzés

Átlag y _{ysz X}^ˆ_/N^ˆ_ Adatátvétel nagyobb területről

Értékösszeg _Y^ˆ Y^ˆ_sz X_iY^ˆ/X^ˆ Xi a kisterület sokasági adata

Az átlagok és az értékösszegek szintetikus becsléseit értelemszerűen használhat- juk utólagosan rétegzett mintáknál az egyes rétegekben. A szintetikus becslések szó- rásnégyzete rendszerint kis értékű, a torzításuk becslésére viszont nincs minden eset- ben jól működő szabály, az átlagos négyzetes hibájukat ezért csak hozzávetőlegesen tudjuk becsülni. Emiatt napjainkban már csak ritkán használnak szintetikus becslés- sel előállított adatokat.

3.2. Kombinált becslések

Az előző szakasz feltételei mellett legyen ˆ₁

Yi egy kisterületi értékösszeg becslése és ˆ₂

Yi egy ehhez tartozó szintetikus becslés. Tegyük fel, hogy az MSE jól becsülhető.

A két becslés

K 1 2

ˆ_i ˆ_i (1 )ˆ_i

Y  λY   λ Y , 0 λ 1

súlyozott számtani átlagát kombinált becslésnek nevezzük. λ értékét úgy kell meg- határozni, hogy a kombinált becslés átlagos négyzetes hibája minimális legyen. Ele- mi szélsőérték-számítási eljárással a következő eredmény adódik:

 

ˆK

MSE Yi közel minimális, ha

2

1 2

(ˆ ) 1

ˆ ˆ 1 ,

( ) ( )

i

i i i

MSE Y

λ  MSE Y MSE Y  F

 

ahol F_i MSE Y(ˆ_i1) MSE Y(ˆ_i2). A kombinált becslés nyilván stabilabb, mint a szinte- tikus, hiszen amennyiben

ˆ

₁

Y

i szórása nagy, akkor

ˆ

_K

Y

i az

ˆ

₂

Y

i komponenshez lesz közel, ha viszont

ˆ

₂

Y

i torzító hatása nagy, akkor

ˆ

_K

Y

i

ˆ

₁

Y

i -hez lesz közel. Az MSE szempontjából

ˆ

_K

Y

i mindkét komponensénél jobb. Gyakorlati alkalmazása szem- pontjából jelenleg nem tartozik a gyakran használt módszerek közé.

3.3. Struktúrát megőrző becslések

A SPREE (structure preserving estimator – struktúrát megőrző becslés) módszer a Deming–Stephan [1940] szerzőpáros nevével fémjelzett, iteratív arányos közelíté-

sek módszerének az általánosítása. Tekintsünk egy három kategóriaváltozóval defi- niált háromdimenziós kereszttáblát, amelynek celláiban nemek, korcsoportok és földrajzi egységek által meghatározott x_ijk létszámú személyek vannak. Tekintsük továbbá az

1 Nk

ij ijk

k

a x







, i 1, 2, j 1, 2, ..., N_j, /3/

1 Nj

ik ijk

j

b x







, i 1, 2, k1, 2, ..., N_k, /4/

2

1

jk ijk

i

c x







, j 1, 2, .., N_j, k 1, 2, ..., N_k /5/

egyenleteket. a_ij az -ediki nemhez ( i 1 vagy 2) és a -edikj korcsoporthoz tartozó személyek száma a tekintett tartományban (ország vagy régió stb.), b_ik az

-edik

i nemhez tartozó és a -adikk településen lakó személyek száma, és végül c_jk a -edik

j korcsoporthoz tartozó és a -adikk településen lakó személyek száma.

Ha az x_ijk cellagyakoriságok népszámlálási adatok, akkor a /3/–/5/ egyenletekkel definiált a_ij, b_ik és c_jk teljesítik a rájuk vonatkozó konzisztencia-feltételeket. Te- gyük fel, hogy a népszámlálás után például egy mikrocenzusban megállapították a_ij és b_ik továbbszámított értékét, az a_ij és a b_ik létszámokat, a c_jk továbbszámított értékét azonban nem, és most ez a feladat. A következőképpen célszerű eljárni.

Meg kell oldani az összesen 2  N_j 2  N_k számú

1 k ˆ

N

ij ijk

k

a x



 



és

1 j ˆ

N

ik ijk

j

b x



 



/3’/–/4’/

egyenletet kalibrálás segítségével, induló értékként a népszámlálási x_ijk gyakorisá- gokat választva. Minden egyes iterációs lépés két részből áll:

1. lépés: a 2  N_j 2  N_k egyenlet helyett még csak ugyanennyi

„<, >” vagy éppen „=”jelet viselő relációnk van, ezek jobb oldalát megfelelő r_ij vagy r_ik tényezővel szorozzuk, hogy egyenlőséget kap- junk;

2. lépés: minden (i, j, k) indexhármasra össze kell gyűjteni az aktuális ˆx_ijk r_ij vagy r_ik szorzóját, ezek egyszerű számtani átlaga lesz ˆx_ijk frissített egységes korrekciós tényezője. Következik az 1.

lépés.

Az eljárás konvergál, ˆx_ijk végleges értékét az /5/ egyenlőségbe helyettesítve megkapjuk c_jk becslését:

2

1

ˆ_jk ˆ_ijk

i

c x







, j 1, 2, .., N_j, k 1, 2, ..., N_k.

A /3/–/5/ egyenletekből a megoldotthoz hasonló további feladatokat lehet megfo- galmazni. A SPREE-módszer nem ment ki a divatból, a „Stakeholders of Statistics”

című konferencia (Budapest, 2016. október) egyik előadásában is hivatkoztak rá.

Mivel speciális kalibrálási eljárásnak tekinthető, a becsült adatok szórásnégyzetének becslését a Deville–Särndal [1992] által javasolt eljárással határozhatjuk meg; a SPREE ennek következtében olyan kisterületi módszer, amely közelítően torzítatlan becslést eredményez.

3.4. EBLUP-eljárás: kisterületi adatok becslése lineáris regressziós modellel

Tekintsük a következő modellt:

ˆ θ_i  x β^T_i  z v_{i i}e_i i 1, 2, …, ,m /6/

ahol

– x_i a magyarázóváltozók, β pedig a regressziós együtthatók p dimenziós vektora,

– v_ia kisterületi hatások azonos eloszlású független változói nulla várható értékkel és σ²szórásnégyzettel,

– e_i nulla várható értékű mintavételi hibák adott ψ_i szórásnégyzet- tel, egymástól és a v_i hatásoktól függetlenek,

– θ_i= ( )g Y_i az y változó átlagának transzformált értéke az i-edik kisterületen,

– ˆ a θ_i θ_i mintából származó közvetlen becslése, – z_i adott pozitív konstans.

Tekintsük θ_i következő becslését:

θ_i^H γ θ_{i i}ˆ (1   γ_i)x β^T_i , /7/

ebben γ_i z σ_i^{2 2} (z σ_i^{2 2}ψ_i) és a regressziós együtthatók (β β₁, ₂, ..., β_p)^T β becs- lését a /6/ egyenlet

ˆθ Xβ +Zv  e, /6’/

mátrix-vektor alakjából számítjuk az általánosított legkisebb négyzetek módszerével (Aitken algoritmusa) a következőképpen. A /6’/ egyenletben X m  p méretű mát- rix, sorai az x₁^T,x₂^T, ..., x^T_m és Z  diag( , , ..., z₁ z₂ z_m), V  ( , , ..., v₁ v₂ v_m)^T, …,

1 2

( ,e e , ..., e_m)^T



e .

Ezekkel a jelölésekkel

1 1 1ˆ

( ^TΩ^ )^{ T}Ω^



β X X X θ,

ahol Ω  diag((z σ₁^{2 2}ψ₁), (z σ₂^{2 2}ψ₂), ..., (z σ_m^{2 2}ψ_m)). Az általánosított legki- sebb négyzetek módszerének tulajdonságaiból következően β a β legjobb lineáris torzítatlan becslése, az általánosan használt angol rövidítéssel BLUE.

A /6/ modell hierarchikus modell összevont alakban. Az első lépcsője

ˆ θ_i  θ_i e_i, i 1, 2, …, ,m /6a/

ez azt mutatja, hogy ˆθ_i mintavételes eljárásból származik, a második,

θ_i  x β^T_i  z v_{i i}, i 1, 2, …, ,m /6b/

a „tulajdonképpeni” modell. A /6b/ modell mellett a /7/ képlettel adott θ_i^H a θ leg- jobb lineáris torzítatlan becslése, átlagos négyzetes hibája pedig

MSE θ( _i^H) E θ( _i^Hθ_i)² γ ψ_{i i} (1 γ_i)²x^T_i (X^TΩ^1X)^1x_i, /8/

γi definíciója a /7/ képlet mellett szerepelt. Ez az állítás a /7/ és a /6b/ képletekből kiindulva bizonyítható, az ehhez szükséges számítás azonban nem triviális.

A /6/ modell az általános lineáris vegyes modellcsalád speciális esete, az idetar- tozó modellek valamennyien az

y  Xβ  Zv e /9/

modell speciális esetei. Itt

– y a mintából származó megfigyelések n dimenziós vektora;

– X adott n  p dimenziós, maximális rangú mátrix;

– Z adott n h dimenziós, maximális rangú mátrix;

– v független eloszlású, nulla várható értékű kisterületi hatások vektora, G kovariancia mátrixszal;

– e független eloszlású, nulla várható értékű mintavételi hibák vek- tora, R kovariancia mátrixszal.

G és R bizonyos δ δ₁, , ..., ₂ δ_q paraméterek függvényei. y variancia-kovariancia mátrixa R  ZGZ^T. A θ_i^H becslés BLUP tulajdonságára vonatkozó eredmény a /9/

modell megfelelő eredményének speciális esete.

A /6/ modell kisterületi szintű, a kisterületekhez tartozó, közvetlenül a mintából származó ˆθ_i becsléseket dolgozza fel. A /9/ modell speciális esetei között vannak mintavételi egység szintűek is, ezen belül olyanok is vannak, amelyek többlépcsős mintavételből származó megfigyeléseket is tudnak kezelni. Találkozhatunk olyan modellváltozatokkal is, amelyeknél az MSE =E θ( _i^Hθ_i)² számításánál csak a min- tavételi terv szerinti várható értéket veszik figyelembe, így az eredmények jobban megfelelnek a mintavételes statisztikusok ízlésének. A kisterületi átlagok becslése mellett az értékösszegek becslésére alkalmas modellek sem hiányoznak, és végül figyelmet fordítanak az idősorok kisterületi problémáira is.

A szórásnégyzetek, kovarianciák δ paraméterektől való függése a részletesen tárgyalt /6/ modellnél, illetve a megoldását jelentő θ_i^H becslésnél sem hiányzik, igaz, hogy itt csak egyetlen ilyen paraméter van, éspedig δ  σ (helyesebb lett volna a σ_i jelölést használni). Az /7/ és /8/ explicit képletek σ ismerete nélkül nem használha- tók, becslésére több módszer is létezik. Az egyik legegyszerűbb a Fay–Herriot [1979] által javasolt iterációs eljárás.

Legyen ^{2 2 2}

1

( ) ( ˆ ) ( )

m

T

i i i i

i

h σ θ ψ z σ







x β  , β β (σ²), σ² értékét a ( 2)

h σ  m  p egyenletből határozzuk meg. σ^2,10, t 1, 2, 3, … esetén legyen

2, 1 2, 2, 2,

( ( )) ( )

t t t t

σ ^ σ  m  p h σ h σ ,

ahol h σ( ²) a h σ( ²) függvény deriváltjának közelítése. Ha a határérték nem pozitív, legyen σˆ² 0. A kapott határértéket a /7/ egyenlőségbe helyettesítve – β kifejezésébe is – megkapjuk az elméleti BLUP-becslés, θ_i^H gyakorlatban használható értékét:

ˆθ_i^H γ θˆ_{i i}ˆ (1  γˆ_i)x β^T_i ˆ. /10/

A ˆθ_i^H becslést empirikus BLUP-, röviden EBLUP-becslésnek nevezzük. Reg- ressziós módszerek alkalmazásánál a /9/ modellcsalád egyedeinél hasonló módon járunk el, mint a /6/ modellnél, ideértve a BLUP–EBLUP terminológia használatát is.

A δ δ₁, , ..., ₂ δ_q paraméterek meghatározására a maximum likelihood és a korlátozott maximum likelihood módszereket lehet alkalmazni.

3.5. Kisterületi adatok becslése empirikus Bayes-módszerrel lineáris regressziós modell alapján

A becslési eljárás alapja ugyanaz a hierarchikus modell, mint az előző alfejezet- ben, azzal a különbséggel, hogy az egyes lépcsőkben a hibatagok viselkedésével kapcsolatban a normális eloszlást feltételezzük:

ˆ |θ θ_{i i} N θ ψ( ,_{i i}), i 1, 2, ..., ,m θ_i N(x β^T_i , z^{2 2}_iσ ), i 1, 2, ..., m. /11/–/12/

Válasszuk a /11/–/12/ modell mellett a prior megoszlást úgy, hogy a posterior várható értéke ^ˆ_i^B a ^ˆ_i és az x β^T_i kifejezések súlyozott átlaga legyen:

E θ θ( _i| , ,ˆ_i βσ²)  γ θ_{i i}ˆ (1  γ_i)x β^T_i . /13/

Mivel a ^ˆ_i^B várható érték a θ_i változó Bayesi becslése, a _i paraméter értékét az a követelmény határozza meg, hogy ^ˆ_i^B átlagos négyzetes hibája, E θ( _i^Bθ_i)² mini-

mális legyen. Vonjuk le ^ˆ_i^B-ből – /13/ jobb oldalából – a θ_i γ θ_{i i}  (1  γ θ_i) _i kifejezést. A különbség négyzete

B 2 2 ˆ 2 ˆ 2 2

(θ_i θ_i)  γ θ_i( _i θ_i)  2 (1 γ_i  γ θ_i)( _i θ_i)(x β^T_i θ_i) (1 γ_i) (x β^T_i θ_i) =

2 2 2 2

2 (1 ) ( ) (1 ) ( )

i i i i i i i i i i

γ e γ γ e z v γ z v

       .

Az utóbbi kifejezés várható értéke E θ( _i^Bθ_i)² γ ψ_i² _i  0 (1  γ_i)^{2 2 2}z σ_i , ami- nek minimuma a γ_i  z σ_i^{2 2} (z σ_i^{2 2}ψ_i) helyen van, a minimum értéke pedig γ ψ_{i i}. Ezzel tehát kimutattuk ^ˆi^B minimumtulajdonságát, átlagos négyzetes hibája, γ ψ_{i i} kisebb, mint az előző alfejezetben taglalt θ_i^H BLUP-tulajdonságú becslésé (lásd a /8/

képletet). Ennek az az oka, hogy ^ˆi^B meghatározásánál a /11/–/12/ feltételeket az i indexnek csak egy értéke mellett használtuk, így nincs bizonyíték arra, hogy a reg- ressziós paraméterek konkrétan meg nem határozott β vektora a többi egyenlet ese- tén is megfelelő lett volna. Ezt a hiányosságot úgy korrigálhatjuk, hogy a számításo- kat i 1, 2, , m értéke mellett elvégezve, az i számú, különböző σ és β para- métereket alkalmas iterációs eljárással „összefésüljük”. A /10/ összefüggéshez ha- sonlót kapunk:

ˆ_iEB ˆ_{i i}ˆ (1 ˆ_i) ^T_i ˆ θ  γ θ   γ x β.

Ezt a becslést az empirikus BLUP = EBLUP mintájára empirikus bayesi becslés- nek mondjuk és az „EB” rövidítéssel jelöljük; mind az értéke, mind pedig az átlagos négyzetes hibája megegyezik a θ_i^Hbecslés EBLUP változatának, a ˆθ_i^H becslésnek a megfelelő adataival (/10/ képlet). Az MSE becsléséhez bizonyos esetekben numeri- kus integrálásra vagy Monte-Carlo-módszerekre van szükség.

3.6. Konjugált prior alkalmazása lineáris regressziós modell esetén

A címben szereplő eljárás a hierarchikus Bayes-módszer változata (Wikipedia [2017]), amelynél a prior eloszlás választásának lehetőségét leszűkítjük oly módon, hogy a sűrűségfüggvény típusa hasonlítson a mintaelemeknek az alkalmazott model- lel meghatározott eloszlásához. A bayesi szigorral szembeni engedménynek pozitív hozadéka, hogy ezáltal a posterior eloszlás függvényére explicit képletet kapunk. A prior kiválasztásától eltekintve az eljárás egyes lépései megfelelnek a hierarchikus Bayes-módszer szabályainak.

Tekintsük az

 

y Xβ e

modellt, ahol y (y₁, , ..., y₂ y_n) ,^T X n k méretű, k  n, maximális rangú mátrix, e  ( , , ..., ) ,e e_{1 2} e_n ^T e_i N(0, σ²), i 1, 2, , n. Szükségünk lesz β legkisebb négyzetek szerinti becslésére: βˆ ( X X^T )^1X y^T .

A modellhez a következő likelihood függvény tartozik:

2 2 /2

2

( | , , ) ( ) exp 1 ( ) ( )

2

n T

f σ σ

σ

  

    



y X β  y Xβ y Xβ .

Annak érdekében, hogy a β – βˆ különbség normális eloszlású legyen, a likelihood függvényt a következőképpen alakítjuk át:

( | , , 2) f y X βσ 

2 /2 2 2 2 ( )/2

2

1 ˆ ˆ

( ) exp( 2 )( ) exp ( ) ( )( ) ,

2

v n v T T

σ vs σ σ

σ

    

    



 β β X X β β

ahol vs² ( y Xβˆ)^T( y Xβˆ) és v  n k.

Ennek alapján a következő prior eloszlást célszerű választani:

2 2 2

1( , ) 2( ) 3( | ) f βσ  f σ f β σ , ahol

f σ₂( ²)(σ²)^{v02 1} exp(v s_{0 0}² 2σ²) egy inverz gamma eloszlás sűrűségfüggvénye,

 

2 2 2 2

3( | ) ( ) ^k exp ( 0)^T 0( 0) 2

f β σ  σ ^  β μ Λ β  μ σ

pedig egy normális eloszlásé, ahol Λ₀ ˆβ kovariancia-mátrixának inverze, az ún.

pontossági mátrix, ami itt a prior hatékonyságának a jelzőszáma. A μ₀ megválasztá- sa a felhasználó feladata.

Az előzők birtokában a posterior eloszlását a következőképpen adhatjuk meg:

2 4( , | , )

f βσ y X  f( | , , y X β σ²) f σ₂( ²)f₃( |β σ²)

2 /2 2 2 2 2

0 0 0

(σ )^n exp( (  y  Xβ) ( ^T y Xβ) 2σ )(σ )^k exp( (  β  μ )^TΛ( β  μ ) 2σ )

(0 1)

2 2

(σ )^{ a } (exp( b σ0 ))

  ;

az a₀ v₀ 2 és b₀ v s_{0 0}² jelölések az eloszlás inverz gamma típusára hívják fel a figyelmet.

Az f₄( ,βσ²,| , )y X sűrűségfüggvény formailag tovább egyszerűsíthető, de erre itt nincs szükség. Ennek az alfejezetnek – miképpen a dolgozat egészének is – a célja nem az Olvasó felhasználói ismeretekkel való felruházása, hanem – ahogyan a beve- zetőben is említettük – a kisterületi becslések módszertanának rövid bemutatása.

3.7. Empirikus Bayes-módszer alkalmazása (0, 1) értékű megfigyelések esetén, segédváltozók nélkül

Tekintsük a következő kétlépcsős, mintavételi egység szintű hierarchikus modellt:

| FAE

ij i

y p ⁵ Bernoulli (p_i), y_ij  0 vagy 1, j 1, 2, ..., ,n_i i 1, 2, ..., m /14/

p_i FAE béta ( , )α β , α  0, β  0. /15/

A feladat ebben az esetben a p_i valószínűségek becslése az y_ij megfigyelések alapján. A béta eloszlás sűrűségfüggvénye

1 1

( )

( | , ) (1 ) ,

( ) ( )

β α

i i i

Γ α β

f p α β p p

Γ α Γ β



 

   α 0, β  0,

Γ(.) a gamma függvény. Feladatunkban (f p α β_i| , ) lesz a konjugált prior elosz- lás sűrűségfüggvénye. Az y_i 



_jy_ij jelölés bevezetésével az első lépcső így mó- dosul:

( _i| _i) ⁱ _i^yi(1 _i)^{n yi}

i

f y p n p p

y

  

     , /14’/

5 Független azonos eloszlású.

( _i| _i)

f y p a likelihood sűrűségfüggvénye. (f p α β_i| , )és ( | )f y_i p_i szorzatából a

 ,,

| i i y

p FAE béta(yii,niyii)

összefüggés adódik, tehát az előzőekkel összhangban a posterior is béta eloszlást követ. A pontbecslés a következő:

pˆ ( , ) _i^B α β  E p y α β( _i| _i, , )  (y_i α) ( n_i  α  β) és /16/

2

( )( )

( | , , )

( 1)( )

i i i

i i

i i

y α n y β

Var p y α β

n α β n α β

  

      .

Gyakorlati célokra az α és β paramétereket a momentumok módszerével be- csülhetjük, éspedig a ˆp 



_kn p_kˆ_k



_kn_k mintabeli súlyozott átlag és az

2 ( (ˆ ˆ)2

p k k k k k

s 



n p  p



n súlyozott szórásnégyzet segítségével, a következő egyenletekből:

ˆ ˆ( ˆ) ˆ, α α β  p



^{αˆ  βˆ 1}



^1

   



kn s^{k p2} p^{ˆ 1} p m^{ˆ 1}



p^ˆ

^

¹ p^ˆ

^ 

kn^k ln^l2 kn^k

^

m ¹

^ 

¹

 

             .

Ezeket az értékeket /16/-ba helyettesítve a pˆ_i^EB γ p_iˆ_i  (1 γ p_i)ˆ becslést kap- juk, ami az ˆα és a ˆβ véletlen ingadozását is tükrözi, γ_in_i (n_i α βˆ ˆ). pˆ_i^EB átla- gos négyzetes hibáját jackknife módszerrel lehet meghatározni.

3.8. Empirikus Bayes-módszer alkalmazása a Poisson-gamma-modellre

Tekintsük a következő kétlépcsős modellt:

yi Poisson (e θ_{i i}), azaz (f y e θ_i| _{i i})  exp(e θ e θ_{i i})( _{i i})^yi y_i!

θi FAE gamma ( , ),v α azaz ₁( | , ) ¹ ( )

i

v αθ v

i i

f θ α v α e θ

Γ v

 

 ,

1, 2, ..., ,

i  m m a kisterületek száma, y_i a mintából származó közvetlen becslés, értéke – Poisson-eloszlásról lévén szó – nem negatív egész. A Poisson-eloszlást likelihoodnak, a gamma eloszlást pedig – az empirikus Bayes-módszernek megfele- lően – prior eloszlásnak választva, a posterior eloszlás sűrűségfüggvényét ezek sűrű- ségfüggvényének szorzataként határozhatjuk meg:

2( _i| _i, , )

f θ y v α =f y e θ f θ α v( _i| _{i i}) (_{1 i}| , )exp(  ( α  e θ θ_i) )_{i i}^{vyi 1} .

Mint látható, f θ y v α₂( _i| _i, , ) is egy gamma eloszlás sűrűségfüggvénye y_i  v és ei α paraméterekkel. Ennek következtében θ_i bayesi becslése és a posterior szó- rásnégyzete

ˆ B ( | , , )

i

i

i i

i

y v

θ E θ y α v

e α

  

 , illetve ( | , , ) ₂.

( )

i

i i

i

y v

Var θ y α v

e α

 

 /17/–/18/

Az α és v paraméterek értékére explicit kifejezés nincsen, közelítő eljárás vi- szont többféle is van. A legegyszerűbb közelítő eljárás a következő. Határozzuk meg a mintából származó értékek súlyozott átlagát és az ehhez tartozó szórásnégyzetet a következőképpen:

 

. .

ˆ_e 1 _{i i} ˆ ,_i

θ e e θ

 m s^e2.



i



e e^{i .}

 

θ^ˆi θ^ˆe.



²,

ahol e_. 



_ie_i, és oldjuk meg a v αˆ ˆθˆ_e. és a v αˆ ˆ² s_e²_. θ eˆ_{e. .} egyenletrendszert az ˆα és a ˆv ismeretlenekre. Az így meghatározott értékeket maximum likelihood iterációkkal lehet finomítani. ˆα és a ˆv értékét a /17/ egyenlőségbe helyettesítve a következő, „igazi empirikus” bayesi becslést kapjuk:

EB B

ˆ_i ˆ_i ( , ) ˆ ˆ ˆ_{i i}ˆ (1 _i)ˆ_e., θ θ α v  γ θ  γ θ ahol ˆγ_i  e e_i ( _i αˆ). ˆ^EB

θi egy a mintából származó közvetlen becslés és egy szin- tetikus becslés súlyozott átlaga.

Irodalom

DEMING,W.E.–STEPHAN,F.F. [1940]: On a least squares adjustment of a sampled frequency table when the expected marginal totals are known. The Annals of Mathematical Statistics. Vol. 11.

No. 4. pp. 427–444. https://doi.org/10.1214/aoms/1177731829

DEVILLE,J.-C.–SÄRNDAL,C.-E.[1992]: Calibration estimators in survey sampling. Journal of the American Statistical Association. Vol. 87. No. 418. pp. 376–382.

https://doi.org/10.2307/2290268

FAY,R.E.–HERRIOT,R.A. [1979]: Estimates of income for small places: an application of James- Stein procedures to census data. Journal of the American Statistical Association. Vol. 74. Issue 366a. pp. 269–277. https://doi.org/10.1080/01621459.1979.10482505

RAO, J. N. K. [2003]: Small Area Estimation. John Wiley & Sons, Inc. Hoboken.

http://dx.doi.org/10.1002/0471722189

RAO,J.N.K.–MOLINA,I. [2015]: Small Area Estimation. 2^nd Edition. John Wiley & Sons, Inc.

Hoboken. http://dx.doi.org/10.1002/9781118735855

WIKIPEDIA [2017]: Bayesian Linear Regression. https://en.wikipedia.org/wiki/Bayesian_linear_

regression