Eötvös Loránd Tudományegyetem Természettudományi Kar
Statisztikai kísérlettervezés és elemzés BSc Szakdolgozat
Készítette:
Ungár Anikó Matematika BSc, Matematikai elemz® szakirány
Témavezet®:
Pröhle Tamás Valószín¶ségelméleti és Statisztika Tanszék
Budapest 2016
Tartalomjegyzék
1. Bevezetés
2
2. Véletlenítés
4
2.1. Zajfaktorok hatásának minimalizálása . . . . . . . . . . . . . . . . . . . . . 2.2. Véletlenen alapuló próbák . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 5
3. Teljesen véletlenített tervek
11
4. Többfaktoros eljárások
20
5. Tobit-modellek
28
6. Összegzés
38
3.1. A teljesen véletlenített tervek modelljei . . . . . . . . . . . . . . . . . . . . . 11 3.2. Elemzés az ANOVA módszerrel . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.3. Kontrasztok . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.1. A kétfaktoros eset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2. Faktoriális elemzés . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.1. Tobin eredeti modellje . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 5.2. Általánosítások . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 5.3. A sampleSelection csomag . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
1
1. fejezet
Bevezetés A tudományos kutatásokban és az ipar különböz® területein kitüntetett fontosságúak a kísérletek. Kísérleteket azért végzünk, mert bizonyos kérdésekre szeretnénk választ kapni. Tipikus kérdések lehetnek például a következ®k: • Biztonságos és hatásos-e egy újonnan kifejlesztett gyógyszer? • Hogyan készítsünk tápanyagoldatot, ha a termesztett gabona min®ségét szeretnénk
maximalizálni?
• Mennyi egy háztartási gép várható élettartama?
A kísérletezés során kiválasztunk kísérleti egyedeket, amelyekre alkalmazunk valamilyen el®re meghatározott eljárásokat, majd rajtuk az eljárások hatását mérjük. A mért eredményeket válasz változóknak nevezzük, és ezek elemzése után vonunk le következtetéseket. Lényeges különbség a kísérletek és a meg gyelések eredményein alapuló tanulmányok között, hogy a kísérleteket mi magunk irányítjuk, beállításait szándékosan változtatjuk annak érdekében, hogy a változás hatását vizsgálni tudjuk, a meg gyelések eredményeit viszont nem tudjuk befolyásolni. Meg gyelések adatai például az id®járási tényez®k, a pénzügyi vagy demográ ai adatok. Olyan kísérleteket szeretnénk végezni, amelyek eredményeib®l értelmes és helyes következtetéseket vonhatunk le a kutatásunkkal kapcsolatban. Ezen a ponton azonban számos kérdés felmerül. Hogyan lehet az eljárásokat és azok hatását jól elkülöníteni? Milyen módszerekkel lehet a válaszokat hatékonyan elemezni? Mekkora a hiba a kísérletben? Ehhez hasonló kérdésekre adnak választ a statisztikai kísérlettervezés (angolul Design of Experiment ) módszerei. Irányított kísérleteket nem minden kutatásnál lehet végezni, ugyanakkor a modellalkotáshoz és a meg gyelt adatok elemzéséhez ilyen esetben is jól használhatóak a matematikai statisztika eszközei. Szakdolgozatom célja, hogy bemutasson néhányat a legfontosabb kísérlettervezési és elemzési módszerek közül. A 2. fejezet a véletlen jelent®ségér®l szól. Amellett, hogy véletlenítéssel csökkenteni lehet a kísérlet hibáját, a véletlen hipotézisek önmagukban is megfelel® alapot adhatnak a következtetésekhez. Ezt konkrét példán fogom bemutatni. A 3. fejezetben már általánosan írok a véletlenített kísérlettervekr®l. Ismertetem a két leggyakrabban vizsgált modellt és becslést adok a modellek paramétereire, ezután pedig részletezem, hogy miért jól m¶köd® elemzési módszer az ANOVA. Az elemzés kiegészítéseként írok röviden a kontrasztokról is. A 4. fejezetben azzal b®vülnek a kísérletek, hogy már nemcsak egy tényez® befolyásolja az eredményeket, hanem több tényez® különböz® beállításainak hatását vizsgálom. Az 5. fejezet eltér az ®t megel®z®kt®l, mert nem irányított kísérletekr®l szól, 2
hanem egy meg gyelések eredményein alapuló gazdasági modellt és annak általánosításait mutatja be. A dolgozatban szerepl® ábrák és programok mindegyike az R statisztikai programmal készült, illetve a vizsgált adathalmazok is megtalálhatók az R-ben ([12]).
3
2. fejezet
Véletlenítés Egy kísérletet meghatároznak a vizsgált eljárások, a kísérleti egyedek, a mért válasz változók, valamint az egyedeket az eljárásokhoz rendel® módszer. Egy kísérlet véletlenített, ha az egyedeket és az eljárásokat véletlenszer¶en rendeljük egymáshoz. A kísérlet eredményét befolyásoló tényez®ket faktoroknak nevezzük, és a kísérlet során ezek beállításait változtatjuk. A beállítások az adott faktor szintjei. A tervezés kezdetén meghatározzuk azokat a faktorokat, amelyek várhatóan legjobban befolyásolják majd a kísérletek eredményét. A többi faktort zajfaktornak tekintjük, és ezek hatását minimalizálni szeretnénk. Ebben a fejezetben két olyan a problémáról írok, melyeket - ha felmerülnének a tervezés során - a véletlen segítségével oldhatunk meg. Az irányított kísérletekr®l szóló fejezetek Gary W. Oehlert 2010-es könyvén ([1]) alapulnak, nagy segítséget jelentettek továbbá Pröhle Tamás és Zempléni András tanár urak jegyzetének ([2]) kapcsolódó fejezetei. 2.1. Zajfaktorok hatásának minimalizálása
A hatások keveredésének azt nevezzük, amikor egy faktornak vagy eljárásnak a hatása nem különböztethet® meg egy másik faktornak vagy eljárásnak a hatásától. A véletlenítésnek könnyen látható, ugyanakkor igen fontos szerepe van abban, hogy ezt a problémát megel®zzük. Például egy kísérletben kétfajta írógépet hasonlítunk össze, és azt szeretnénk megállapítani, hogy melyikkel lehet gyorsabban dolgozni. Tíz titkárn®nek kell ugyanazt a szöveget legépelnie. A mért válasz a gépeléshez szükséges id®. Hogyan végezzük el a kísérletet? Dönthetünk úgy, hogy öt titkárn® az els®, öt pedig a második fajta gépen dolgozik, és a két csoportban összehasonlítjuk a válaszok átlagát. Ekkor azonban a titkárn®ket gyorsaságuk szerint véletlenszer¶en kell az egyes csoportokba sorolni, hiszen ha az öt leggyorsabb kerülne egy csoportba az öt leglassabbal szemben, akkor nem lehetne eldönteni, hogy például az els® csoportban azért kaptunk jobb válasz értékeket, mert az els® típusú gép jobb, vagy azért, mert a gyorsabb titkárn®k bármelyik gépen gyorsabban dolgoznak. Dönthetünk úgy is, hogy mindkét gépen legépelik a szöveget egymás után. Ekkor az a kérdés, hogy melyikkel kezdjék. Ha a második gépen jobb eredményeket kapunk, akkor gondolhatunk arra, hogy már ismert volt a szöveg és ezért gyorsabban tudtak haladni. Ha rosszabb eredményeket kapunk, annak az is lehet az oka, hogy fáradtabbak, amikor másodjára kell gépelniük. Mivel nem tudjuk, hogy melyik tényez® befolyásolja az eredményeket, ebben az esetben is véletlenszer¶en kell kiválasztani öt titkárn®t, aki az els® gépen kezd, és a többi pedig a második gépen. 4
Ebben a példában a számunkra érdekes faktorok a gépek, a zajfaktorok pedig a titkárn®k tulajdonságai (gyorsaság, szöveg memorizálása, fáradtság). A véletlenítés abban segít, hogy olyan faktorok torzító hatását is kisz¶ri, amelyekr®l esetleg nem is tudtuk, hogy fontosak, és emiatt nem vettük volna gyelembe ®ket. Éppen ezért a véletlenítés alkalmas módszer a zajfaktorok hatásának minimalizálására. 2.2. Véletlenen alapuló próbák
A kísérletek eredményét elemz® legtöbb módszer t -próbákat, illetve F -próbákat használ, melyek azon a feltételezésen alapulnak, hogy a hibák, és így a teljes adathalmaz független, és normális eloszlást követ 0 várható értékkel és azonos szórással. A függetlenséget például χ2 -próbával ellen®rizhetjük, a normalitás vizsgálatára pedig alkalmas többek között a Kolmogorov-Szmirnov-teszt, az Anderson-Darling-próba vagy az adott szinten az el®z®ekhez viszonyítva er®sebb Shapiro-Wilk teszt. Azonban nem biztos, hogy a feltételek minden valós kísérletnél teljesülnek. Ebben az esetben is a segít a véletlen megközelítés: a véletlen nullhipotézisen alapuló próbákat egy konkrét példán keresztül szeretném bemutatni. El®ször de niáljunk néhány fontos alapvet® fogalmat a statisztikai próbákkal kapcsolatban. Legyen adott egy X = (X1 , ..., Xn ) n elem¶ minta, melynek valódi eloszlását Pϑ jelöli, valódi paramétere, ϑ pedig a Θ paramétertér egy ismeretlen eleme. A mintából kiszámolunk valamilyen T (X1 , ..., Xn ) próbastatisztika értéket.
2.2.1. De níció. A hipotézisek: H0 : H1 :
ϑ ∈ Θ0 ϑ ∈ Θ1
ahol Θ = Θ0 ∪ Θ1 a paramétertér diszjunkt felbontása.
2.2.2. De níció. Egy próba elfogadási tartománya χe , kritikus tartománya χk , és χe ∪ χk az összes lehetséges próbastatisztika értékek diszjunkt felbontása.
Ha T (X) ∈ χe , akkor H0 -t elfogadjuk, ha T (X) ∈ χk , akkor H0 -t elutasítjuk.
2.2.3. De níció. Els®fajú hibát követünk el, ha a nullhipotézis igaz, mégis elutasítjuk. Ennek valószín¶sége: Pϑ (T (X) ∈ χk )
ϑ ∈ Θ0 .
Másodfajú hibát követünk el, ha a nullhipotézis hamis, mégis elfogadjuk. Ennek valószín¶sége: Pϑ (T (X) ∈ χe )
ϑ ∈ Θ1 .
2.2.4. De níció. Ha az els®fajú hiba kisebb, mint α, akkor az α szintet a próba terjedel-
mének nevezzük.
2.2.5. De níció. A próba er®függvénye annak a valószín¶sége, hogy nem követünk el má-
sodfajú hibát:
β(ϑ) = 1 − Pϑ (T (X) ∈ χe ) = Pϑ (T (X) ∈ χk )
ϑ ∈ Θ1 .
2.2.6. De níció. A p-érték annak a valószín¶sége, hogy a nullhipotézis igaz volta mellett
egy másik kiszámolt próbastatisztika érték legalább annyira extrém, mint a meg gyelt érték. 5
2.2.7. De níció. Legyen
X1 , ..., Xn minta Fϑ eloszlásból, ahol ϑ a keresett valós paraméter. A (T1 (X), T2 (X)) intervallumot 1 − α megbízhatóságú kon dencia-intervallumnak
nevezzük, ha
∀ϑ : Pϑ (T1 (X) < ϑ < T2 (X)) ≥ 1 − α.
2.2.8. De níció. Legyen a tatlan becslés ϑˆ-ra, ha
ϑ paraméternek ϑˆ egy becsült értéke. Egy T0 (X) érték torzíEϑ (T0 (X)) = ϑˆ ∀ϑ ∈ Θ.
A következ® az R programban is megtalálható Student's Sleep Data adathalmaz, mely Student 1908-as cikkében ([4]) szerepelt példaként. A táblázat kétféle gyógyszer altató hatását mutatja (az alvással töltött órák számának növekedése a kontroll adatokkal összehasonlítva) tíz páciensen. ID 1 2 3 4 5 6 7 8 9 10
extra(group1) 0.7 -1.6 -0.2 -1.2 -0.1 3.4 3.7 0.8 0.0 2.0
extra(group2) 1.9 0.8 1.1 0.1 -0.1 4.4 5.5 1.6 4.6 3.4
Az adatokat célszer¶ mindig vizuálisan is megjeleníteni, mert el®zetes információkat ez alapján is tudunk szerezni. A 2.1 ábra például sugallja, hogy a 2-essel jelölt altatószer hatása nagyobb, mint az 1-essel jelölté. Kétféle próbával fogom tesztelni, hogy megegyezike a két csoportban a válaszok átlaga.
2.1. ábra. Box-plot a sleep adathalmazra 6
Megoldás páros t -próbával: Akkor használhatjuk a páros t -próbát, ha a függetlenségre és normalitásra vonatkozó feltételek teljesülnek. Ekkor kivonjuk a 2-es esetben kapott válaszból az 1-es esetben kapott választ minden páciensre, és az így nyert 10 új értékre egymintás t -próbát végzünk. Jelölje az n különbséget d1 , d2 , ..., dn , és feltételezzük, hogy ezek független, normális eloszlású valószín¶ségi változók µ várható értékkel és azonos szórásnégyzettel, melyet σ 2 jelöl. A nullhipotézisünk az, hogy a várható érték egy el®re meghatározott értékkel egyenl®, most µ0 = 0-val, és ezzel szemben kétoldali alternatív hipotézis áll: H0 : µ = µ 0 = 0 H1 : µ 6= 0.
A statisztikai próbáknál mindig a nullhipotézist tartjuk fontosabbnak, ezért az els®fajú hibát tekintjük súlyosabbnak, ennek valószín¶ségét írjuk el®. Tetsz®leges mintaelemszámra és terjedelemre megtalálhatjuk a t -eloszlás megfelel® értékét ([5], 8. táblázat). A nullhipotézist elutasítjuk, ha a kapott t -érték nagyobb a táblázatbelinél. Az egymintás t -próba képlete t=
d¯ − µ0 √ s/ n
ahol d¯ az adatok (itt a meg gyelt különbségek) átlaga, az s pedig a korrigált tapasztalati szórás: v u u s=t
n
1 X ¯ 2. (di − d) n−1 i=1
A különbségek: d = (1.2, 2.4, 1.3, 1.3, 0.0, 1.0, 1.8, 0.8, 4.6, 1.4), az átlag d¯ = 1.58, µ0 = 0, n = 10 és s ≈ 1.23, ahonnan t ≈ 4.0621. 10 elem¶ mintára 95%-os biztonsági szinten t = 2.262, 98%-os biztonsági szinten t = 2.821, és még 99%-os biztonsági szint esetén is t = 3.250. Ez azt jelenti, hogy az els®fajú hiba valószín¶ségét maximum 0.01-nek választva, ezáltal a másodfajú hiba valószín¶ségét jelent®sen növelve is elutasítjuk a nullhipotézist, ami szigni káns eredmény. A páros t -próba eredménye az R-ben is megkapható, x jelöli az 1-es, y pedig a 2-es esetben mért válaszokat: Paired t-test data: y and x t = 4.0621, df = 9, p-value = 0.002833 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: 0.7001142 2.4598858 sample estimates: mean of the differences 1.58
A kapott p ≈ 0.0028 érték er®sen szigni káns eredmény, jelent®s bizonyíték H0 ellen, a kapott 95%-os kon dencia-intervallum szintén arra utal, hogy nem igaz a nullhipotézis. 7
Megoldás véletlen próbával: Miben különbözik a fentiekt®l egy véletlenítésen alapuló vizsgálat? A véletlen nullhipotézis az, hogy a kétféle gyógyszer teljesen ekvivalens, és a mért válaszokból nem következtethetünk arra, hogy melyik gyógyszert használta a páciens. Például az els® páciens válasz értékei 0.7 és 1.9 voltak, mely értékeket mi megcímkéztünk úgy, mint az 1-es és a 2-es gyógyszer hatása. A véletlen nullhipotézis alapján az els® páciens válaszai 0.7 és 1.9 lettek volna, függetlenül attól, hogy melyik gyógyszert nevezzük 1-esnek és melyiket 2-esnek. A válaszok különbségét tekintve ez azt jelenti, hogy az el®bbi páros t -próbánál az 1.9−0.7 = 1.2 értékkel számoltunk, de ugyanekkora valószín¶séggel kaphattunk volna 0.7 − 1.9 = −1.2 értéket is, és ugyanez igaz a többi különbségre is. Ezért a páros t -próbához hasonló, de véletlen nullhipotézisen alapuló vizsgálatnál a különbségek abszolút értékét vesszük, és mindegyikhez véletlenszer¶en rendelünk hozzá egy el®jelet. Az el®jelek függetlenek egymástól, és ugyanakkora a pozitív és a negatív el®jel valószín¶sége. A véletlen próba elvégzéséhez választunk valamilyen statisztikát: például a véletlen el®jellel vett értékek összegét vagy átlagát tekintjük, és meghatározzuk ezeknek az értékeknek az eloszlását. Ha ezen eloszlás mellett a meg gyelt érték extrémnek mondható, akkor elutasítjuk a nullhipotézist. Bármilyen statisztikát is választunk, a véletlen el®jellel vett értékek eloszlásához meg kellene határoznunk az adott statisztika értéket minden lehetséges el®jel-kombinációra. Kétféle lehet az el®jel minden különbség esetén, így bármilyen x elemszámú mintára 2x értéket kellene kiszámolnunk, ez pedig még viszonylag kis mintaelemszámokra, számítógéppel dolgozva is hosszú id® lehet. Ha túl sok számítást kellene elvégeznünk, akkor nem az összes értékb®l határozzuk meg, hanem véletlen mintavételezéssel becsüljük az eloszlást. Ezen okok miatt az ilyen típusú vizsgálatok általában nem a legkönnyebben megvalósíthatóak, ugyanakkor megvan az az el®nyük, hogy semmilyen plusz feltételnek nem kell teljesülnie az alkalmazhatóságukhoz. Általában tehát célszer¶ a hagyományos módszereket alkalmazni, akkor érdemes a véletlen nullhipotézisen alapuló próbákra gondolni, amikor a függetlenség vagy normalitás nem teljesül. A konkrét feladat megoldásában a t -statisztika értékek eloszlását vizsgálom. A megoldás R-beli kódja a következ®: X <- matrix( c(0.7, 1.9, -1.6, 0.8, -0.2, 1.1, -1.2, 0.1, -0.1, -0.1, 3.4, 4.4, 3.7, 5.5, 0.8, 1.6, 0.0, 4.6, 2.0, 3.4),byrow=TRUE,ncol=2) xy <- as.numeric(t(X)) xi <- seq(1,by=2,length=10) yi <- (1:20)[-xi] x <- xy[xi]
8
y <- xy[yi] t.orig <- t.test(y,x,alt=c("two.sided"),paired=TRUE)$stat set.seed(0) n <- 1000 t.rnd <- rep(NA,n) for(k in 1:n) { xi <- seq(1,by=2,length=10)+rbinom(10,1,.5) yi <- (1:20)[-xi] t.rnd[k] <- t.test(xy[yi],xy[xi],alt=c("two.sided"), paired=TRUE)$stat } hist(t.rnd,breaks=100) abline(v=t.orig,col="red",lwd=2)
Az adatokat az X mátrix tárolja, els® oszlopában az 1-es, a másodikban a 2-es altatószer esetében mért válaszok. A meg gyelt t -értéket rögzítjük a t.orig nev¶ változóban. Ezután n = 1000-szer kiszámolunk egy véletlen t -értéket (t.rnd) úgy, hogy most a kiválasztott xi indexhez véletlenszer¶en vagy 0-t vagy 1-et adunk. Ezáltal minden páciens mért válaszait véletlenszer¶en címkézzük 1-esnek vagy 2-esnek, így a különbségük el®jele véletlen lesz. Az 1000 darab véletlen t -érték eloszlását a 2.2 ábra mutatja.
2.2. ábra. A véletlen t -értékek eloszlása a sleep adathalmazra 9
A függ®leges piros egyenes a t.orig érték, mely a nullhipotézis igaz volta mellett extrém érték. A véletlen p -értéket úgy határozhatjuk meg, hogy a kapott, t.orig-nál nagyobb értékek számát elosztjuk az összes kiszámolt érték számával. Ez a konkrét példánkban 0, aminek az az oka, hogy a meg gyelt különbségek mindegyike pozitív vagy 0, tehát bármely más el®jel-kombinációra csak kisebb értéket kaphatunk. A nullhipotézist tehát a véletlenen alapuló próba után is elutasítjuk. Általában igaz, hogy a véletlen p -értékek közel vannak a függetlenséget és normális eloszlást feltételez® próbák p -értékeihez. Emiatt érdemes a véletlen nullhipotézisen alapuló próbákat használni, amikor ezek a feltételek nem teljesülnek.
10
3. fejezet
Teljesen véletlenített tervek Az el®z® fejezetben láttuk, hogy a zajfaktorok hatásának minimalizálása érdekében lényeges, hogy az eljárások és a kísérleti egyedek egymáshoz rendelése véletlenszer¶en történjen. Azokat a kísérleteket is érdemes véletleníteni, amelyeknél ez els®re nem t¶nik fontosnak, hiszen minimális munkával sok kés®bbi problémát megel®zhetünk. A legegyszer¶bb véletlenített kísérletterv több eljárás összehasonlítására a teljesen véletlenített terv (Completely Randomized Design ), ennek struktúrájáról és elemzésér®l fogok írni ebben a fejezetben. 3.1. A teljesen véletlenített tervek modelljei
A teljesen véletlenített tervek struktúrája a következ®. Adott g összehasonlítandó eljárás és N kísérleti egyed, melyek úgy osztunk el, hogy 1. megadjuk a csoportok méretét, jelölje ezeket n1 , n2 , ..., ng úgy, hogy n1 +n2 +...+ng = N teljesül; 2. véletlenszer¶en kiválasztunk n1 egyedet, amelyekre az 1. eljárást alkalmazzuk, a maradék N − n1 közül véletlenszer¶en kiválasztunk n2 -t, amelyekre a 2. eljárást alkalmazzuk, és így tovább. Ebben az elosztásban az N egyed minden hozzárendelése minden n1 , n2 , ..., ng méret¶ csoporthoz egyforma valószín¶ség¶. Fontos, hogy a teljes véletlenítés csak azt határozza meg, hogy az egyedeket hogyan rendeljük az eljárásokhoz, de természetesen egy kísérletnek más véletlenített tulajdonságai is lehetnek. A teljesen véletlenített kísérlettervek standard elemzésénél azt szeretnénk megállapítani, hogy a különböz® kísérleti csoportokban mért válaszok várható értékei megegyeznek-e, és ha nem, akkor melyek térnek el egymástól. A továbbiakban feltételezzük, hogy a hibák független, normális eloszlású valószín¶ségi változók 0 várható értékkel és azonos szórásnégyzettel. Jelölje yij a j -edik válasz változót az i-edik kísérleti csoportban (i = 1, ..., g, j = 1, ..., ni ). Az várható értékekre kétféle modellt fogok vizsgálni. Az els® modellben azt feltételezzük, hogy minden kísérleti csoportban azonos a válaszok várható értéke, melyeket µ-vel jelölünk. Ezt nevezzük redukált modellnek. A hibára tett feltételezésekb®l az is következik, hogy a teljes adathalmaz független, normális eloszlású, azonos szórásnégyzettel melyet σ 2 fog jelölni. A modell tehát: yij ∼ N (µ, σ 2 ).
Más formában felírva: yij = µ + εij
11
(3.1)
ahol εij valamilyen független, normális eloszlású zaj 0 várható értékkel és σ 2 szórásnégyzettel. A második modellben azt feltételezzük, hogy a válaszok várható értéke különböz® minden eljárás esetén, jelölje ezeket az értéket µi . Ezt teljes modellnek nevezzük. (A redukált modell valójában ennek egy speciális esete.) A feltételezéseink alapján yij ∼ N (µi , σ 2 ),
másképpen: yij = µi + εij .
(3.2)
A csoportonként különböz® átlagot a következ®képpen is kifejezhetjük: µ i = µ ∗ + αi
(3.3)
ahol a konstans µ∗ egyfajta általános hatás, αi pedig az i. kezelés hatása. Ez a felírás megszorítások nélkül túlparaméterezés, az elemzés során viszont az αi hatások vizsgálata is fontos lesz, err®l a 3.3. szakaszban írok. A µ∗ -ot többféleképpen is megválaszthatjuk. Egy klasszikus választás a csoportonkénti várható értékek átlaga, ezzel azt mondjuk, hogy a csoportok hatása átlagosan 0: g X µi
g
(3.4)
αi = 0.
(3.5)
∗
µ =
i=1
ebb®l a választásból azonban következik, hogy g X i=1
Ennek egy alternatívája az, hogy a µ∗ a csoportonkénti várható értékek súlyozott átlaga, a súlyok pedig az egyes csoportok elemszámai. A csoportok hatása átlagosan 0, gyelembe véve a csoportokban lév® egyedek számát. ∗
µ =
g X n i µi i=1
N
(3.6)
Az αi értékekre ez is megszorítást jelent: g X
ni αi = 0.
(3.7)
i=1
Az ismeretlen paramétereket becsülni szeretnénk. El®ször vezessünk be néhány jelölést a szokásos módon. A • mindig a megfelel® alsó indexen való összegzést jelöli. A meg gyelések összege az i. kísérleti csoportban: yi• =
ni X
yij .
(3.8)
j=1
A meg gyelések átlaga az i. kísérleti csoportban: ni 1 X yi• y¯i• = yij = . ni ni j=1
12
(3.9)
Az összes meg gyelés összege: y•• =
g X ni X
g X
yi• ,
(3.10)
g ni 1 XX y•• yij = . N N
(3.11)
yij =
i=1 j=1
i=1
és az összes meg gyelés átlaga: y¯•• =
i=1 j=1
A becsléshez a legkisebb négyzetek módszerét alkalmazzuk. A 3.2 egyenlet alapján felírhatjuk a teljes modellt mátrixos formában: (3.12)
Y = XM + ε,
ahol Y ∈ RN ×1 a meg gyelt yij értékeket tartalmazó mátrix, M ∈ Rg×1 a csoportonkénti várható értékek vektora, X ∈ RN ×g együtthatók mátrixa, amit tervmátrixnak is nevezünk, ε ∈ RN ×1 pedig a hibák mátrixa.
y11 1 : : y1n1 1 y21 0 : : y2n = 0 2 . . .. .. yg1 0 : : ygng
0 : 0 1 : 1
··· : ··· ··· : ···
0 : 0
··· : ···
.. .
0
.. .
0 ε11 : : 0 ε1n1 0 µ ε 1 21 : µ2 : . + 0 .. ε2n2 . .. .. . µg εg1 1 : : 1
(3.13)
εgng
Az M mátrix becslése a 3.12 egyenlet alapján (3.14)
c = (X T X)−1 X T Y M
Mivel
1 0 XT X = . .. 0
··· ···
1 0 ··· 0 1 ···
0 ··· 1 ···
0 ··· 0 ···
··· ···
0 0 ···
0 ···
1 ···
.. .
.. . ···
.. . ···
.. . ···
1 : 1 0 0 0 : 0 .. . . . 1 . 0 : 0
0 : 0 1 : 1
··· : ··· ··· : ···
0 : 0
··· : ···
.. .
.. .
0 : 0 0 n1 0 · · · 0 n2 · · · : = . .. . . 0 .. . . .. . 0 0 ··· 1 : 1
0 0 .. . ng
így (X T X)−1 diagonális mátrix, amelynek f®átlójában a csoportok elemszámának reciprokai szerepelnek. A további számolásokat elvégezve kapjuk, hogy c= M
y1• n1
y2• n2
13
...
yg• ng
(3.15)
Minden csoport várható értéke tehát az adott csoportban meg gyelt értékek átlaga. Így a becslések: µ ˆi = y¯i• (3.16) (3.17)
µ ˆ = y¯••
A 3.6 egyenletb®l µˆ∗ =
g X ni µ ˆi i=1
N
=
g X ni y¯i• i=1
N
=
y•• = y¯•• . N
(3.18)
Ez ugyanaz a becslés, amivel a µ értéket becsültük, ezért elhagyom a ∗ megkülönböztetést, és a továbbiakban a µ jelölést használom mindkét értelemben. A kezelések hatásai, αi -k felírhatók αi = µi − µ alakban, tehát a becslés: (3.19)
α ˆi = µ ˆi − µ = y¯i• − y¯•• .
Az utolsó becsülend® paraméter a σ 2 , amit a következ®képpen írhatunk fel: Pg 2
σ ˆ =
i=1
Pni
j=1 (yij
N −g
− y¯i• )2
.
(3.20)
Azért N − g szerepel a nevez®ben, mert g lineáris összefüggés van az N érték között. A fenti becslések mind torzítatlan becslések a megfelel® paraméterre.
Pontbecslésnek nevezzük azt, amikor egy ismeretlen paramétert egyetlen számértékkel becsülünk, ahogy a fentiekben tettük. Egyetlen konkrét szám azonban nem ad információt arról, hogy mennyire jó a becslésünk, mennyire ingadozhat a becsült érték a valódi érték körül. Ezért fontos intervallumbecsléseket is készíteni. A vizsgálatunkban a µ paraméter értéke a legfontosabb, ezért ehhez készítek most példaként ε megbízhatóságú kon dencia-intervallumot. Mivel a teljes adathalmaz független és normális eloszlású, ezért a µ becslésére használt y¯•• is normális eloszlású µ várható értékkel σ és √ szórással, az N
y¯•• − µ √ σ/ N
pedig már standard normális eloszlású. Mivel a standard normális eloszlás eloszlásfüggvénye az y tengelyre szimmetrikus, ezért a kon dencia-intervallumot a 0-ra szimmetrikusan keressük −uε , +uε kon dencia-határokkal. A 2.2.7 de níció alapján y¯•• − µ √ ≤ uε = 2φ(uε ) − 1 = 1 − ε P − uε ≤ σ/ N
(3.21)
ahol φ a standard normális eloszlás eloszlásfüggvénye. Az egyenl®tlenség-rendszer rendezése után kapjuk, hogy σ σ P y¯•• − uε √ ≤ µ ≤ y¯•• + uε √ = 1 − ε. N N
(3.22)
Fontos, hogy a kon dencia-intervallumnak ez a becslése csak akkor jó, ha a σ elméleti szórás ismert. Amikor csak a mintából kiszámolt s∗ tapasztalati szórás ismert, akkor az el®bbieknek megfelel®en y¯•• − µ √ s∗ / N
14
változót írhatjuk fel. Err®l megmutatható, hogy N − 1 szabadsági fokú Student-eloszlású, tehát y¯•• − µ √ ≤ tε = 1 − ε, (3.23) P − tε ≤ σ/ N
s∗
P y¯•• − tε √
s∗ ≤ µ ≤ y¯•• − tε √ = 1 − ε. N N
(3.24)
A szükséges uε és tε értékeket a megfelel® függvények táblázatában találhatjuk meg, például Lukács Ottó könyvében ([5]). 3.2. Elemzés az ANOVA módszerrel
Az elemzés során a csoportonkénti várható értékek egyez®ségér®l szeretnénk információt nyerni. A legismertebb ilyen elemzési eljárás az ANOVA (Analysis of Variance ), amely a szórások vizsgálatával keresi a választ a várható értékek egyez®ségére. A nullhipotézisünk az, hogy a válaszok várható értéke minden csoportban megegyezik. A 3.2 egyenlet alapján ez ekvivalens azzal, hogy minden kezelés hatása megegyezik. Ezért a nullhipotézis felírható H0 :
α1 = ... = αg = 0
formában, az alternatív hipotézis pedig: ∃i : αi 6= 0.
H1 :
A döntéshez kétféle távolságnégyzet becslést készítünk, amelyeket majd F -próbával hasonlítunk össze. Távolság alatt egyel®re számértékek különbségét értjük. Az els® a csoporton belüli távolságnégyzet összeg, ez azt jelenti, hogy minden egyes adatpontnak az ® csoportjának átlagától vett négyzetes eltérését összegezzük. Ez valójában a hiba négyzete a megfelel® csoportban, ezért SSE -vel jelölöm (Error Sum of Squares ): g X ni X SSE = (yij − y¯i• )2 .
(3.25)
i=1 j=1
A másodikat csoportok közötti négyzetösszegnek nevezzük, itt az egyes csoportok átlagának az összes meg gyelés átlagától vett négyzetes eltérését összegezzük. Ez tulajdonképpen a kezelések hatásának négyzetösszege az egyes csoportok elemszámával súlyozva ezért SST rt vel fogom jelölni (Treatment Sum of Squares ): SST rt =
g X ni X
2
(¯ yi• − y¯•• ) =
i=1 j=1
g X
2
ni (¯ yi• − y¯•• ) =
i=1
g X
ni α ˆ i2 .
(3.26)
i=1
A két négyzetösszeg összege éppen a teljes távolságnégyzet összeg (az egyes válaszok eltérése az összes válasz átlagától), ezt SST -vel fogom jelölni (Total Sum of Squares ): SST =
g X ni X
(yij − y¯•• )2 =
i=1 j=1
=
g X ni X
(yij − y¯i• + y¯i• − y¯•• )2 =
i=1 j=1 g X ni X i=1 j=1
2
(yij − y¯i• ) + 2
g X ni X i=1 j=1
15
g X ni X (yij − y¯i• )(¯ yi• − y¯•• ) + (¯ yi• − y¯•• )2 i=1 j=1
Mivel minden csoportban a válaszok csoportátlagtól való eltérésének összege 0, ezért g X g ni ni X X X (¯ yi• − y¯•• ) (yij − y¯i• )(¯ yi• − y¯•• ) = (yij − y¯i• ) = 0 i=1 j=1
i=1
j=1
tehát a kétszeres tag kiesik, és (3.27)
SST = SSE + SST rt
teljesül. Vizsgáljuk most a fenti négyzetösszegeket más szemszögb®l. Az N darab yij meg gyelés értéke egy N -dimenziós pontnak felel meg, ezt jelölte Y . A teljes modell szerint a válaszok várható értéke csoportonként különböz®, tehát egy Yˆ szintén N -dimenziós pontot kapunk, melynek az els® n1 koordinátája az els® csoportban mért válaszok átlaga, a következ® n2 koordinátája a második csoportban mért válaszok átlaga, és így tovább. A redukált modell szerint nincs különbség az egyes csoportok várható értékei között, tehát egy olyan N -dimenziós pontot kapunk, melynek minden koordinátája az összes meg gyelés átlaga. Jelölje ezt a pontot Y¯ .
y¯11 : y¯1n 1 Y = ... y¯g1 : y¯gng
y¯1• : y¯1• ˆ Y = ... y¯g• : y¯g•
y¯•• : y¯•• ¯ Y = ... y¯•• : y¯••
Így az Y Yˆ szakasz hossza SSE , az Yˆ Y¯ szakasz hossza SST rt , Y Y¯ hossza pedig SST . Az X tervmátrix i. oszlopa az RN egy ni dimenziós alterének térátlója. A g altér térátlói által kifeszített alteret G-vel jelölöm, és a teljes modell szerinti altérnek nevezem. A redukált modell szerinti altér, melyet T -vel jelölök, a teljes N -dimenziós tér térátlója által kifeszített altér. Ekkor teljesül, hogy T ⊂ G ⊂ RN
továbbá igaz, hogy Yˆ az Y vetülete a G-n, Y¯ pedig egyrészt Y vetülete a T -n, másrészt Yˆ vetülete T -n. Ezekb®l következik, hogy Y Yˆ ⊥Yˆ Y¯
tehát Y Yˆ Y¯ derékszög¶ háromszög, és a Pitagorasz-tétel alapján Y Yˆ 2 + Yˆ Y¯ 2 = Y¯ Y 2 ,
ami a fentebb kiszámolt SST = SSE + SST rt egyenl®séggel ekvivalens. Mivel minden kísérleti csoportban a csoportátlagtól való eltérések összege 0, ezért bármely ni − 1 közülük meghatározza az utolsót. Másképpen: a hiba ni − 1 szabadsági fokú minden P csoportban, és i (ni − 1) = N − g szabadsági fokú a teljes adathalmazban. Feltettük 16
korábban, hogy a kezelések átlagos hatása 0, tehát a g eljárás közül bármelyik g − 1 meghatározza az utolsót, azaz az eljárások szabadsági foka g − 1. Vezessük be az M ST rt =
SST rt , g−1
M SE =
SSE N −g
(3.28)
jelöléseket. M ST rt és M SE függetlenek egymástól, és mindkett® χ2 eloszlású rendre g−1 és N − g szabadsági fokkal. Ezért a hányadosuk a nullhipotézis igaz volta mellett F-eloszlású g − 1, N − g szabadsági fokokkal: F =
M ST rt ∼ Fg−1,N −g . M SE
(3.29)
A nullhipotézis tehát F -próbával ellen®rizhet®: F < Fα esetén elfogadjuk, F > Fα esetén elutasítjuk, ahol Fα a fenti F-eloszlás 1 − α kvantilise. Amikor az alternatív hipotézis igaz, akkor az F -statisztika nem centrális F-eloszlást követ. A nem centrális F -eloszlás szabadsági fokai megegyeznek az eddig vizsgált (centrális) F eloszlás szabadsági fokaival, azonban van még egy nemcentralitási paramétere is, melyet ζ jelöl és a következ®képpen adhatjuk meg: 2 i ni αi , σ2
P ζ=
(3.30)
ahol az αi értékek az egyes csoportok tényleges hatásai. A nemcentralitási paraméter azt méri, hogy milyen messze vannak a kezelési átlagok az egyenl®ségt®l, az y¯i• szórásához viszonyítva. A centrális F -eloszlásra ζ = 0, és minél nagyobb a ζ értéke, annál valószín¶bb, hogy elutasítjuk H0 -t. Az alábbi ANOVA táblázat az R programban megtalálható Chicken Weights by Feed Type adathalmazra vonatkozik, a 3.1 ábrán pedig a hozzá tartozó boxplot látható. Analysis of Variance Table Response: weight Df Sum Sq Mean Sq F value Pr(>F) feed 5 231129 46226 15.365 5.936e-10 *** Residuals 65 195556 3009 --Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
A kísérlet során azt vizsgálták, befolyásolja-e a csirkék növekedését az, hogy milyen táplálékkal etetik ®ket. Az adathalmaz N = 71 meg gyelést tartalmaz, ennyi csibét osztottak véletlenszer¶en g = 6 kísérleti csoportba, és minden csoportban különböz® táplálékot kaptak az egyedek. A válaszok a hat hét elteltével mért tömegek.
17
3.1. ábra. Box-plot a chickwts adathalmazra
Az ábra is sugallja, hogy nem igaz a nullhipotézis, a táblázatból pedig kiolvashatóak a korábban meghatározott szabadsági fokok: a faktoroké g − 1 = 5, a hibák szabadsági foka (reziduális) N − g = 65. Ilyen szabadsági fokokra az F-próba kritikus értéke 2.37 körüli, 95%-os biztonsági szinten ([5], 9. táblázat). A kapott F -érték 15.365 >> 2.37, ami alapján elutasítjuk a nullhipotézist és azt a következtetést vonjuk le, hogy a csirkék növekedésének mértéke jelent®sen függ attól, hogy milyen tápanyagot kapnak. 3.3. Kontrasztok
Az ANOVA módszerrel megállapíthatjuk, ha nem ugyanaz a válaszok várható értéke minden kísérleti csoportban. Azonban egy ilyen elemzésb®l azt nem tudjuk meg, hogy mely eljárások különböznek a többit®l, és milyen mértékben. Szükséges tehát vizsgálni a kezelések hatását. Erre a vizsgálatra alkalmasak a kontrasztok ([6]). Ha a csoportok várható értékére valamilyen speciális tulajdonságot feltételezünk, ami számértékekkel kifejezhet®, akkor a nullhipotézis kontrasztként megadható. Például egy kísérletben három eljárás hatását szeretnénk összehasonlítani. Azt feltételezzük, hogy a második csoportban mért válaszok átlaga kétszerese lesz az els® csoportban mért válaszok átlagának, a harmadik csoportban pedig az átlag háromszorosa lesz az els® csoport átlagának. Ez kontrasztokkal kifejezve:
18
C1 1
C2 2
C3 3
A nullhipotézist kifejez® számértékeket kontraszt együtthatóknak nevezzük akkor, ha az összegük 0. A fenti példában az együtthatók összege 6, átlaguk 2. Ha mindegyik értékb®l kivonjuk az átlagot, akkor megkapjuk a nullhipotézisünkhöz tartozó kontraszt együtthatókat: C1 −1
C2 0
C3 1
3.3.1. De níció. A kontraszt kísérleti csoportok átlagának a lineáris kombinációja. A következ®képpen jelöljük: g X
wi µi = w({µi })
i=1
ahol a wi együtthatók kielégítik a
Pg
i=1 wi
= 0 feltételt.
De niálhatjuk a kontrasztokat úgy is, mint a csoportok hatásainak lineáris kombinációja. Jelölje ezt a lineáris kombinációt w({αi }). Mivel a kontraszt együtthatók összege 0, így minden rögzített x konstansra teljesül: w({αi }) =
g X
wi α i = x
i=1
g X
wi +
g X
i=1
wi αi =
i=1
g X
wi (x + αi )
i=1
és x = µ helyettesítéssel adódik, hogy w({αi }) = w({µi }).
(3.31)
Tehát valóban az eljárások hatásait hasonlítjuk össze. A nullhipotézis akkor igaz, ha a fenti lineáris kombinációk értéke 0 körüli. Az R-ben négyféle beépített kontraszt található: 1. Helmert-kontraszt (contr.helmert): a vizsgált faktor egy szintjét mindig az el®z® szintek átlagával hasonlítja össze. 2. Polinomiális kontraszt (contr.poly): azt keresi, hogy milyen rendben függ a válasz a vizsgált faktor hatásától. 3. Összeg típusú kontraszt (contr.sum): az utolsó kontraszt együttható az el®z®ek összegének −1-szerese. 4. Eljárás típusú kontraszt (contr.treatment): a faktor minden szintjét egy el®re megadott alapszinttel hasonlítja össze. A fentieken kívül mátrix formájában tetsz®leges kontraszt megadható.
19
4. fejezet
Többfaktoros eljárások Az el®z® fejezetben bemutatott teljesen véletlenített terveknél g eljárást rendeltünk véletlenszer¶en az N egyedhez, azonban az eljárások szerkezetér®l eddig nem mondtunk semmit. Faktoriális struktúráról beszélünk akkor, amikor a g eljárás mindegyike kett® vagy több faktor szintkombinációiból tev®dik össze. Ezeket a kombinált eljárásokat faktor-szint kombinációknak vagy faktoriális kombinációknak nevezzük. A véletlenítést továbbra is használjuk, csak most már az eljárásokra úgy gondolunk, hogy faktoriális szerkezet¶ek. Ha többfaktoros a kísérlet, akkor vizsgálhatjuk a faktorok hatását külön, egymás után, azonban sokkal hasznosabb faktoriális elemzést végezni. Ennek okáról és az elemzés módszereir®l írok b®vebben ebben a fejezetben. 4.1. A kétfaktoros eset
Kezdésként a kétfaktoros eljárásokat vizsgálom. A két faktort A és B jelöli, a faktoroknak rendre a és b szintje van. A válaszokat yijk jelöli, ahol yijk a k-adik válasz abban az eljárásban, amit az A faktor i-edik és a B faktor j -edik szintje határoz meg (i = 1, ..., a és j = 1, ..., b). A válaszokat 4.1. táblázatnak megfelel®en lehet általánosan megjeleníteni. 1.
2.
y111 : y11n y211 : y21n
y121 : y12n y221 : y22n
:
:
:
...
a.
ya11 : ya1n
ya21 : ya2n
...
1. 2.
...
... ...
b. y1b1 : y1bn y2b1 : y2bn : yab1 : yabn
4.1. táblázat. Válasz változók kiegyensúlyozott kétfaktoros eljárás esetén
4.1.1. De níció. Egy többfaktoros kísérlet kiegyensúlyozott, ha minden faktor-szint kom-
binációra ugyanannyi választ mérünk.
20
Mostantól feltesszük, hogy az adatok kiegyensúlyozottak, és n választ mérünk minden faktor-szint kombinációra. A kétfaktoros esetben összesen g = ab eljárásunk van, ab − 1 szabadsági fokkal az eljárások között. Az egyes kísérleti csoportokban a várható értéket µij -vel fogom jelölni. A korábbi egyfaktoros kísérlettervekhez hasonlóan y¯ij• a meg gyelt átlag az ij kísérleti csoportban, y¯i•• és y¯•j• a meg gyelt átlaga minden olyan válasznak, amit az A faktor i-edik, illetve a B faktor j -edik szintjén mérünk. El®ször vizsgáljuk külön az A és B faktorokat. Hagyjuk gyelmen kívül a B faktort, és tekintsük a kísérletet, ami így egy teljesen véletlenített terv az A faktorra. Ebben a tervben a eljárás szerepel. Az egyfaktoros eljárások vizsgálatánál kapott eredmények alapján a válaszok átlagának várható értéke az i. sorban µ + αi
(4.1)
amib®l következik, hogy az αi kezelési hatásokra a X
αi = 0
(4.2)
i=1
teljesül. Kiszámolhatjuk ezekre az eljárásokra a csoportok közötti négyzetösszeget, amit SSA -val jelölök, és a − 1 szabadsági fokú lesz. SSA =
a X
bn(ˆ αi )2
(4.3)
i=1
Ugyanígy, ha az A faktor hatását hagyjuk gyelmen kívül, és egy teljesen véletlenített tervet tekintünk a B faktorra, és b eljárásunk van. A válaszok átlagának várható értéke a j. oszlopban µ + βj (4.4) ahol feltesszük, hogy a βj kezelési hatásokra b X
βj = 0.
(4.5)
j=1
A csoportok közötti négyzetösszeget most SSB jelöli, és b − 1 szabadsági fokú. SSB =
b X
an(βˆi )2
(4.6)
j=1
4.1.2. De níció. Az αi és βj hatásokat a megfelel® faktorok f®hatásainak nevezzük, SSA
és SSB pedig a f®hatás négyzetösszegek.
A f®hatás azt a szórást magyarázza, amit egyetlen adott faktor szintjeinek különbsége okoz, vagyis amikor a válasz táblázatban sorról sorra, vagy oszlopról oszlopra lépünk. Azonban a válaszok átlaga függhet az A faktor szintjét®l, amikor a B faktor szintjei között váltunk, és fordítva is ugyanígy.
4.1.3. De níció. Azt mondjuk, hogy a faktorok között kölcsönhatás van, ha az egyik faktor szintjeinek hatása függ a másik faktor szintjeit®l. A kölcsönhatást αβij -vel jelöljük. 21
A kölcsönhatás négyzetösszeget SSAB -vel jelölöm, és a következ®képpen számolhatjuk: SSAB =
a X b X
c )2 . n(αβ ij
(4.7)
i=1 j=1
Példaként most az R program ToothGrowth adathalmazát használom. A kísérletben tengerimalacok fogának növekedését vizsgálták két faktor hatása mellett. A dose faktor a bevitt C-vitamin szintjei (0.5, 1, és 2 mg), a supp faktor szintjei pedig a bevitel módja, vagyis az, hogy narancslevet (OJ) vagy aszkorbinsavat (VC) kaptak a tengerimalacok. A mért válasz a fogak hossza (len). A 4.1 ábrán az látható, hogy a supp faktor egyes szintjein milyen értékeket mértünk a dose faktor szintjeinek függvényében.
4.1. ábra. Coplot a ToothGrowth adathalmazra
Ennél szemléletesebben jeleníthetjük meg a hatásokat, ha úgynevezett interaction plotot, azaz kölcsönhatást mutató ábrát készítünk (4.2 ábra). Ez tulajdonképpen a fenti két gra kon megjelenítése egy ábrán, ebb®l azonban könnyebben tudunk következtetéseket levonni. A konkrét példában az x tengelyen a dose faktor szintjei vannak, és minden szint felett két pont, mely a mért hosszat ábrázolja a supp faktor két szintjén. Ha a megfelel® pontok összekötésével körülbelül párhuzamos szakaszokat kapunk (például a dose faktor 0.5 és 1 értékei között), akkor feltehet®en nincs kölcsönhatás a faktorok között. A 4.2 ábráról 22
viszont azt olvashatjuk le, hogy 1 mg-nál nagyobb dózis esetén az aszkorbinsav jobban befolyásolja a fogak hosszát, mint a narancslé.
4.2. ábra. Interaction plot a ToothGrowth adathalmazra A faktoriális eljárási struktúrák legnagyobb el®nye éppen az, hogy ha van a faktorok között kölcsönhatás, akkor faktoriális elemzéssel ezt a kölcsönhatást becsülni lehet. Az egyszerre csak egy faktor hatását vizsgáló, egymás utáni kísérletek ugyanakkor nem alkalmasak ilyen becslésre, ez pedig komoly hibákhoz vezethet, ha a faktorok közötti kölcsönhatás jelent®s. Annak érdekében, hogy az egyes faktorok f®hatásait megkapjuk, érdemes faktoriális elemzést végezni akkor is, amikor nincs kölcsönhatás. 4.2. Faktoriális elemzés
Az el®z® szakaszban vizsgált hatások alapján a modellünk: yijk = µ + αi + βj + αβij + εijk ,
(4.8)
ahol µ az egyfaktoros terveknél is vizsgált általános átlag, αi az A faktor f®hatása az i. szinten, βj a B faktor f®hatása a j. szinten, αβij a két faktor kölcsönhatása az A faktor i. és a B faktor j. szintjén, εijk pedig független, normális eloszlású, 0 várható érték¶ véletlen zaj. A válasz várható értéke az ij eljárásra: (4.9)
Eyijk = µ + αi + βj + αβij .
A paraméterekre teljesülnek a következ® megszorítások: 0=
a X i=1
αi =
b X
βj =
j=1
a X i=1
23
αβij =
b X j=1
αβij .
(4.10)
A paraméterek becslését az egyfaktoros eset általánosításával kapjuk. Az eredmények: µ ˆ = y¯•••
(4.11)
α ˆ i = y¯i•• − µ ˆ = y¯i•• − y¯•••
(4.12)
βˆj = y¯•j• − µ ˆ = y¯•j• − y¯•••
(4.13)
ˆ = y¯ij• − µ αβ ˆ−α ˆ i − βˆj = y¯ij• − y¯i•• − y¯•j• + y¯••• ij
(4.14)
A számunkra fontos kérdés továbbra is az, hogy megegyezik-e a válaszok várható értéke az egyes kísérleti csoportokban. A kérdés megválaszolásához tehát a továbbiakban is az ANOVA módszert használjuk, most azonban több nullhipotézist kell tesztelnünk. A tesztelés F-próbákkal történik, melyekhez az (4.15)
SSA + SSB + SSAB = SST rt
egyenl®séget használjuk fel, mivel az eljárások hatását bontottuk fel három részre. A 4.2 táblázat összefoglalja, hogy mely nullhipotéziseket kell tesztelnünk, és milyen számolással.
H0
SS
Szabadsági fok
MS
F -érték
∀i : αi = 0 ∀j : βj = 0 ∀i, j : αβij = 0
SSA SSB SSAB
(a − 1) (b − 1) (a − 1)(b − 1)
SSA /(a − 1) SSB (b − 1) SSAB /(a − 1)(b − 1)
M SA /M SE M SB /M SE M SAB /M SE
4.2. táblázat. Nullhipotézisek és ellen®rzésük kétfaktoros eljárások esetén
El®ször a kölcsönhatást érdemes tesztelni. Ha elfogadjuk azt a nullhipotézist, hogy nincs a kölcsönhatás a faktorok között, akkor vizsgálhatjuk a f®hatásokat. A ToothGrowth adathalmaz struktúrájából látjuk, hogy csak a supp szerepel faktorként, a dose pedig numerikus típusú változó. Ez alapján a válaszokra háromféle modellt készíthetünk. Az els® esetben (M1 modell) a dose nem faktor típusú változó. A második esetben (M2 modell) a dose változót faktornak választjuk, és összeg típusú kontrasztot állítunk be hozzá. A harmadik esetben (M3 modell) azt is gyelembe vesszük, hogy a dose valójában rendezett faktor, ilyen esetben pedig polinomiális kontraszt az alapértelmezés. M1 modell: supp faktor, dose numerikus változó:
Ez tehát egy egyfaktoros kísérlet, a faktornak két szintje van, illetve egy lineáris hatás befolyásolja még az eredményt. Ezért az ANOVA táblázat a következ®képpen néz ki:
24
Analysis of Variance Table Response: len Df Sum Sq Mean Sq F value Pr(>F) supp 1 205.35 205.35 12.3170 0.0008936 *** dose 1 2224.30 2224.30 133.4151 < 2.2e-16 *** supp:dose 1 88.92 88.92 5.3335 0.0246314 * Residuals 56 933.63 16.67 --Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
Ebben az esetben a következ® modell írja le a válasz értékeket: M1 <- lm(len~ supp*dose,ToothGrowth) m a b c
<<<<-
coef(M1)[1] c(0,coef(M1)[2]) coef(M1)[3] c(0,coef(M1)[4])
yhat <- fitted(M1) supp <- as.numeric(ToothGrowth$supp) dose <- ToothGrowth$dose ycalc <- m+a[supp]+b*dose+c[supp]*dose round(yhat-ycalc,5) summary(M1)
El®ször is készítünk egy lineáris modellt, majd ennek együtthatóit használjuk fel. A válasz függ egy m konstanstól, melyet a lineáris modellb®l megkaptunk (itt m=11.55), ehhez a supp faktor szintjét®l függ®en ehhez vagy 0-t, vagy szintén a lineáris modellb®l megkapott −8.225 értéket adunk, hozzá kell adnunk továbbá a dose faktor egy lineáris függvényét, és a dose faktor egy supp-tól is függ® lineáris függvényét. Ez jó modell, hiszen az M1 lineáris modell illesztett értékeib®l az ycalc (a számolt modell) értékeit kivonva, öt tizedesjegyre kerekítve 0 értékeket kapunk. M2 modell: supp és dose is faktor: Analysis of Variance Table Response: len Df Sum Sq Mean Sq F value Pr(>F) supp 1 205.35 205.35 15.572 0.0002312 *** dose 2 2426.43 1213.22 92.000 < 2.2e-16 *** supp:dose 2 108.32 54.16 4.107 0.0218603 * Residuals 54 712.11 13.19 --Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1
25
A dose numerikus típusú változót faktor típusúvá alakítva kapjuk a fenti ANOVA táblázatot. A szabadsági fokokra most már a 4.2 táblázatban szerepl® értékeket kaptuk. A modellben a dose faktorhoz összeg típusú kontrasztot állítunk be. ToothGrowth$dose <- factor(ToothGrowth$dose) contrasts(ToothGrowth$dose) <- "contr.sum" M2 <- lm(len~ supp*dose,ToothGrowth) m a b c
<<<<-
coef(M2)[1] c(0,coef(M2)[2]) c(coef(M2)[3],coef(M2)[4],-(coef(M2)[3]+coef(M2)[4])) matrix(c(0,coef(M2)[5],0,coef(M2)[6], 0,-(coef(M2)[5]+coef(M2)[6])),2,3) rownames(c)<-c("OJ","VC") colnames(c)<-c(".5","1","2") supp<-as.numeric(ToothGrowth$supp) dose<-as.numeric(ToothGrowth$dose) c12<-NULL for(j in 1:60) c12<-c(c12,c[supp[j],dose[j]]) yhat <- fitted(M2) ycalc <- m+a[supp]+b[dose]+c12 round(yhat-ycalc,5) summary(M2)
M3 modell: supp faktor, dose rendezett faktor:
Ennél a modellnél nem állítunk be kontrasztot a dose faktorhoz, hanem az alapértelmezett polinomiális kontrasztot vizsgáljuk. ToothGrowth$dose <- ordered(ToothGrowth$dose) M3 <- lm(len~ supp*dose,ToothGrowth) m <- coef(M3)[1] a <- c(0,coef(M3)[2]) b1 <- contr.poly(3)[,1]*coef(M3)[3] b2 <- contr.poly(3)[,2]*coef(M3)[4] c1 <- contr.poly(3)[,1]*coef(M3)[5] c2 <- contr.poly(3)[,2]*coef(M3)[6] supp <- as.numeric(ToothGrowth$supp) dose <- as.numeric(ToothGrowth$dose) yhat <- fitted(M3) ycalc <- m+a[supp]+b1[dose]+b2[dose]+(supp==2)*(c1[dose]+c2[dose]) round(yhat-ycalc,5) summary(M3)
26
Az eredményekb®l az látható, hogy a dose faktornak lényegében csak lineáris hatása van. Fontos megjegyezni, hogy a fentiekben csak a kétfaktoros esetet vizsgáltam, de a módszerek tetsz®leges számú faktorra általánosíthatók. A faktoroknak a f®hatásuk mellett magasabb rend¶ kölcsönhatásai is lehetnek. Ha például A, B, C és D faktorok tartoznak a kísérlethez, melyek f®hatásait rendre α, β , γ , δ jelöli, akkor a kétfaktoros esetben használt jelöléssel a αβ , αγ , αδ , βγ , βδ és γδ két megfelel® faktor kölcsönhatása, ezek a másodrend¶ hatások. A harmadrend¶ hatásokat αβγ , αβδ , αγδ és βγδ jelöli, és azt a hatást fejezik ki, amit a megfelel® három faktor adott szintbeállítása együtt okoz. Ugyanígy a negyedrend¶ hatás, αβγδ azt a hatást fejezi ki, amit a négy faktor beállításai együttesen okoznak. Ha minden lehetséges faktor-szint kombináción elvégezzük a kísérleteket, akkor teljes faktoriális tervekr®l beszélünk. A gyakorlatban viszont, ha sok faktor hatását kell vizsgálnunk, nem mindig lehet minden faktor minden szintjén mérést végezni (túl sok id® kellene hozzá vagy túl költséges lenne). Ilyenkor részleges faktoriális tervekkel dolgozunk, és csak bizonyos el®re meghatározott faktor-szint kombináción végzünk méréseket. Természetesen a kihagyott kísérletek hibákhoz vezethetnek, például ha két hatás minden kísérletben ugyanazon a szinten szerepel (ezt nevezik alias struktúrának ), akkor sehogyan nem lehet ®ket megkülönböztetni. A tervezés szakaszában tehát ezt is gyelembe kell vennünk, és úgy kell beállítanunk a kísérleteket, hogy a hatások lehet®leg ne keveredjenek.
27
5. fejezet
Tobit-modellek A bevezet® fejezetben megkülönböztettem a kísérleteket a meg gyelések eredményein alapuló tanulmányoktól, az azt követ® fejezetekben pedig arra láttunk példákat, hogy az ipari és tudományos kutatásokban hogyan használhatók jól a tervezett kísérletek. Sok esetben azonban nem lehet olyan kísérleteket tervezni, amiben minden információ rendelkezésre áll a kérdéseink megválaszolásához. Például mérnöki kutatásokban fontos kérdés lehet egy gép várható m¶ködési ideje, azaz várható élettartama, de mi csak azt tudjuk meg gyelni, hogy az adott pillanatban m¶ködik-e még a gép. Ugyanígy orvosi kutatásokban is rendkívül fontos kérdés a különböz® kezelések utáni várható élettartam. Ezek a konkrét példák túlélési modellekként is ismertek, szerepük pedig nemcsak az iparban, hanem közgazdasági vagy szociológiai kutatásokban is fontos, amikor olyan jelenségek várható elhúzódását vizsgálják, mint például a munkanélküliség, munkahelyváltások között eltelt id®k, vagy házasságkötés és gyermekvállalás között eltelt id®k. Gazdasági kutatások témája lehet továbbá egy család élelmiszer fogyasztása, amely még az alacsony jövedelemszintek esetén sem csökkenhet egy bizonyos konstans érték alá. Az ilyen kérdéseket vizsgáló modelleknek két f® csoportja van, ezeket csonkított, illetve cenzorált modelleknek nevezzük. Csonkított modellr®l beszélünk akkor, ha a változó értékét valamilyen okoknál fogva egyáltalán nem tudjuk meg gyelni. Az okok elvileg jól de niáltak, ugyanakkor lehetséges, hogy az okok sem meg gyelhet®k. Cenzorált modellekr®l pedig akkor beszélünk, ha a változó valamely okoknál fogva egy konstans értékre változik. Tipikusan ilyen modellek körébe tartoznak az úgynevezett készlet-modellek. A készlet mértékét számtalan körülmény befolyásolhatja (fedezet, igény, stb), ugyanakkor sohasem válhat negatívvá. Az ilyen típusú problémák kezelésére James Tobin 1958-as cikkében ([7]) bevezetett egy jó modellt, melynek azóta sok általánosítása terjedt el. Ebben a fejezetben szeretném bemutatni Tobin eredeti modelljét és annak általánosításait is. Az eddigi fejezetekkel ellentétben most már nem irányított kísérleteket vizsgálunk, hanem meg gyelések eredményeit elemezzük.
5.1. Tobin eredeti modellje
A modell bevezetése el®tt szeretnék de niálni két fogalmat, amelyeket ebben a fejezetben használni fogok.
28
5.1.1. De níció. Legyen X1 , ..., Xn független, azonos eloszlású minta. A likelihood-függvény diszkrét minta esetén
L(ξ, ϑ) = Pϑ (X = ξ) =
n Y
Pϑ (Xi = ξi )
i=1
folytonos eloszlású minta esetén pedig L(ξ, ϑ) = fϑ (ξ) =
n Y
fϑ (ξi )
i=1
ahol fϑ a mintához tartozó s¶r¶ségfüggvény.
5.1.2. De níció. A ϑ paraméter maximum likelihood becslése ϑˆ = T (X) ∈ Θ, ha ˆ = max L(X, ϑ). L(X, ϑ) ϑ∈Θ
Tobin saját modelljét a korlátos függ® változók modelljének nevezte ([7]), és a mikrogazdaságtanban fontos jelenségek leírására alkalmazta. Azt vette észre, hogy a háztartásokra vonatkozó közgazdasági kutatásokban számos változóra igaz, hogy korlátos (lehet alsó vagy fels® korlátja), és a meg gyelt válasz értékek jelent®s része a korlátra esik, a megmaradt válaszok pedig széles intervallumban mozognak a korlát megfelel® oldalán. Ezeket nevezzük korlátos változóknak. Ez a jelenség látható például akkor, amikor azt vizsgálják, hogy a háztartások a jövedelmükt®l függ®en mennyit költenek fogyasztási cikkek bizonyos kategóriáira. Élelmiszert például minden jövedelemszinten kell, hogy vásároljanak, ugyanakkor ha az a kérdés, hogy mennyit költenek luxuscikkekre, akkor az alacsonyabb jövedelemszinteknél számos 0 válaszértéket fogunk kapni. Természetesen 0-nál kevesebbet nem lehet költeni, így a 0 itt alsó korlát. Olyan példát is találunk, ahol az alsó korlát nem szükségszer¶en 0, s®t nem is azonos minden háztartásra: vizsgáljuk azt, hogy az adott háztartásban egy év alatt mennyit változott a megtakarítás összege. A változás lehet negatív is, de abszolút értéke nem lehet nagyobb, mint amennyi az év elején rendelkezésre állt. Az ilyen jelleg¶ kutatásokban az a kérdés, hogy milyen kapcsolat van a korlátos változó és a többi magyarázó változó között. Figyelembe kell vennünk a meg gyelt válaszok összetételét, hiszen egy magyarázó változó befolyásolhatja mind a korlát válaszok valószín¶ségét, mind a nem-korlát válaszok értékeit. Kétféle megközelítés lehetséges: • Ha csak a korlát és nem-korlát válaszok valószín¶ségét vizsgálnánk függetlenül a nem-
korlát válaszok értékeit®l, akkor a probit regresszióval megfelel® statisztikai modellt kapnánk. Ez a mi esetünkben nem hatékony megközelítés, mert így rendelkezésre álló információt hagynánk gyelmen kívül.
• Ha csak a válasz változó értékét vizsgálnánk, és nem vennénk gyelembe, hogy a kor-
lát körül sok válasz tömörül, akkor a többszörös regresszió lenne hatékony statisztikai eljárás. A többszörös regressziós modell feltételei azonban nem teljesülnek, mert a válaszok a korlát értékt®l nem térhetnek el bármilyen irányban.
Tobin tehát a fenti két módszer helyett dolgozott ki egy harmadikat, melyet neve és a hasonló probit modellek után Tobit modellnek neveznek. 29
Legyen W a korlátos függ® változó L alsó korláttal, ami lehet különböz® érték is az egyes háztartásokra. Legyen Y a független (X1 , X2 , ..., Xm ) változók lineáris kombinációja: Y = β0 + β1 X1 + β2 X2 + ... + βm Xm .
(5.1)
A hipotézisünk az, hogy W függ Y -tól. A háztartások eltérhetnek egymástól olyan okokból is, amelyet az Xi független változók és az L alsó korlát nem magyaráz. Ezeket az eltéréseket ε jelöli, melyr®l feltesszük, hogy független, normális eloszlású, 0 várható érték¶, σ 2 szórásnégyzet¶ valószín¶ségi változó. A háztartások viselkedésére feltesszük: ( L ha Y − ε < L, W = Y − ε ha Y − ε ≥ L.
(5.2)
Ez a modell a következ®képpen értelmezhet®: a fogyasztás az L alsó korlátnál kevesebb nem lehet, a háztartás rendelkezésre álló lehet®ségeit az Y írja le, de a család ezzel nincs egészen pontosan tisztában, ezt csak ε hibával értékeli. Emiatt a család fogyasztása akkor is csak L, ha a család rendelkezésre álló forrása L+ε, és ha a család fogyasztása meghaladja az L küszöbértéket, akkor is csak Y − ε-t fogyaszt el a lehet®ségeib®l. Máshogy értelmezve a modell a következ®képpen is magyarázható: a család lehet®sége a státusza alapján Y , ezzel szemben Y − ε-t érzékel, emiatt a fogyasztásának mértékér®l az Y − ε alapján dönt. Ha Y − ε kevesebb, mint a küszöb, akkor a fogyasztás a küszöbbel egyenl®, ha Y − ε nagyobb mint a küszöb, akkor minden lehet®ségüket felhasználják: Y − ε a fogyasztás. Az ε eloszlása alapján a W eloszlása meghatározható a következ®k szerint: P r(W < x) =
( 0 Φ
Y −x σ
ha x ∈ (−∞, L) ha x ∈ [L, ∞)
(5.3)
ahol Φ(x) a standard normális eloszlás eloszlásfüggvényét, Φ(x) = 1 − Φ(x) a standard normális eloszlás túlélésfüggvényét jelöli. A fenti képlet azt is jelenti, hogy P r(W = L) = Φ
Y − L , σ
(5.4)
továbbá a W s¶r¶ségfüggvénye az x > L értékekre: 1 Y − x ϕ . σ σ W várható értéke adott Y és L mellett: Y − L Y − L Y − L E(W ) = LΦ +YΦ + σϕ σ σ σ
ahol felhasználtuk, hogy:
ϕ0 (x) = −xϕ(x).
(5.5)
(5.6) (5.7)
A W és egy vagy több Xi független változó kapcsolatára vonatkozó hipotézisek likelihoodhányados próbával tesztelhet®k. 5.2. Általánosítások
A Tobit-modellek egy részletes összefoglalója Takeshi Amemiya 1984-es cikke ([8]), aki a likelihood függvény alakja alapján sorolja öt osztályba az ismert modelleket. 30
Vezessük be a szükséges jelöléseket a következ® módon: y jelöli a háztartás kiadását valamilyen tartós fogyasztási cikkre, y0 az elérhet® legolcsóbb cikk ára, z minden más kiadás, x pedig a jövedelem. A háztartásokról feltételezzük, hogy maximalizálni szeretnék a hasznosságukat, melyet U (y, z) jelöl. A hasznosságra két megszorítás teljesül: y+z ≤x
az úgynevezett költségvetési megszorítás, míg vagy
y ≥ y0 y=0
a határ megszorítás. Tegyük fel, hogy a költségvetési megszorítással, de a határ megszorítást gyelmen kívül hagyva a hasznossági függvény y ∗ mellett maximális, feltételezve, hogy y ∗ = β1 + β2 x + u, (5.8) ahol u a nem meg gyelhet® változók összessége, melyek befolyásolják a hasznossági függvényt. Ekkor az eredeti probléma megoldása, melyet y jelöl, felírható ( y∗, y= 0, vagy
ha y ∗ > y0 y0 , ha y ∗ ≤ y0 .
(5.9)
Az y0 lehet különböz® érték a háztartásokra. A likelihood függvényt n független meg gyelésre a következ® módon írhatjuk fel: L=
Y
Fi (y0i )
0
Y
fi (yi ),
(5.10)
1
ahol Fi és fi rendre eloszlás- ésQs¶r¶ségfüggvénye, és 0 jelenti a szorzást azokra az i ∗ indexekre, amelyekre yi ≤ y0i és 1 jelenti a szorzást azokra az i indexekre, amelyekre yi∗ > y0i . yi∗
Q
A Tobin által bevezetett eredeti modell lényegében ugyanez volt, ezt fogjuk I. típusú, vagy standard Tobit-modellnek nevezni. I. típusú Tobit-modell:
yi∗ = x0i β + ui , i = 1, 2, ..., n ( yi∗ , ha yi∗ > 0 yi = 0, ha yi∗ ≤ 0
(5.11) (5.12)
ahol ui -k független, normális eloszlású valószín¶ségi változók 0 várható értékkel és azonos szórásnégyzettel. A likelihood függvény az I. modellre: L=
Y Y [1 − φ(x0i β/σ)] σ −1 ϕ[(yi − x0i β)/σ], 0
(5.13)
1
ahol φ és ϕ rendre a standard normális változó eloszlás- és s¶r¶ségfüggvénye. Ez a modell a cenzorált regressziós modellek közé tartozik. Csonkított regressziós modellr®l beszélünk, ha yi∗ ≤ 0 esetén sem xi , sem yi nem meg gyelhet®. Ennek a modellnek a likelihood függvénye: Y L= φ(x0i β/σ)−1 σ −1 ϕ[(yi − x0i β)/σ]. (5.14) 1
31
Az likelihood függvények típusait az alábbi táblázat szerint lehet általánosan összefoglalni. Típus
Likelihood függvény
I. II. III. IV. V.
P (y1 < 0)P (y1 ) P (y1 < 0)P (y1 > 0, y2 ) P (y1 < 0)P (y1 , y2 ) P (y1 < 0, y3 )P (y1 , y2 ) P (y1 < 0, y3 )P (y1 > 0, y2 )
5.1. táblázat. Likelihood függvények Amemiya-féle csoportosítása
A táblázatban minden yj (j = 1, 2, 3) érték N (x0j βj , σj2 ) eloszlású, és P az adott eloszlásból származó valószín¶ség, vagy s¶r¶ségfüggvény-érték. Az y1 értéke alapján határozunk meg kategóriákat, és az azonos kategóriába es® meg gyelések P -értékeit szorozzuk össze. Így például az I. típusú Tobit-modellben P (y1 < 0)P (y1) a rövidített jelölése a Q Q ∗ < 0) P (y f (y1i ) kifejezésnek, ahol f1 az N (x01i β1 , σ12 ) eloszlású valószín¶ségi vál1i 1 0 1 tozók s¶r¶ségfüggvénye.
II. típusú Tobit-modell:
Ez a modell a James J. Heckman nevéhez f¶z®dik. 1979-es cikkében ([9]) a korlátos függ® változók és a csonkított modellek mellett úgynevezett szelekciós modellekkel is foglalkozott. Ezeknek a lényege az, hogy a választ valamilyen meg gyelt változóval (vagy változókkal) szeretnénk magyarázni, ugyanakkor van a háttérben egy olyan változó, ami a modellünket jelent®sen eltorzíthatja. Ezt nevezzük szelekciós változónak. Erre egy egyszer¶ példa, ha azt vizsgáljuk, hogy házaspárok mennyit költenek kulturális programokra (például színház, koncert, kiállítás). A kiadásukról, mint meg gyelt válasz értékr®l feltételezzük, hogy a jövedelmükkel és az iskolázottságukkal magyarázható. Felírhatunk tehát egy modellt, ami ezekt®l a meg gyelt értékekt®l függ. Ugyanakkor a családi helyzetükt®l is függ az, hogy mely házaspárokra lesz jó a modell, hiszen hiába van meg az anyagi lehet®ség és az érdekl®dés, ha például több kisgyermekr®l, esetleg id®s családtagokról kell gondoskodni, akkor lehet, hogy egyáltalán nem költenek kulturális programokra egyszer¶en azért, mert nem szívesen mennek el otthonról. Ebben a példában tehát a a meg gyelt változók a jövedelem és az iskolázottság, a szelekciós változó pedig a családi helyzet, amely kiválasztja, hogy mely házaspárokra lesz jó a modell. A többiek esetében pedig 0 értéket kapunk a magyarázó változók értékét®l függetlenül. A modellt a következ®képpen írhatjuk fel: (5.15) (5.16)
∗ y1i = x01i β1 + u1i , ∗ y2i
=
x02i β2
( ∗, y2i y2i = 0,
+ u2i ∗ >0 ha y1i ∗ ≤ 0, ha y1i
32
i = 1, 2, ..., n
(5.17)
ahol u1i -k és u2i -k kétdimenziós normális eloszlásból származó független változók, melyek várható értéke 0, szórásnégyzetük rendre σ12 és σ22 , a kovarianciájuk pedig σ12 . Feltéte∗ értékeket meg gyelhetjük (ezek a szelekciós változók), y ∗ -t (a tényleges lezzük, hogy y1i 2i ∗ > 0. Feltételezzük továbbá, hogy magyarázó változót) pedig csak akkor ismerjük, ha y1i ∗ ≤ 0. x1i -k meg gyelhet®k minden i-re, de x2i -t nem kell ismernünk, ha y1i Az I. típus ennek tulajdonképpen egy speciális esete volt, amikor a szelekciós változó egyszerre a tényleges magyarázó változó is volt (csak a háztartások jövedelmét®l függött, hogy ∗ lehet negatív is. A mennyit költenek). Az I. típussal ellentétben a magyarázó változó, y2i likelihood függvény L=
Y
∗ P (y1i ≤ 0)
Y
0
ahol
∗ f (y2i |y1i>0 )P (y1i > 0),
(5.18)
1
azon tényez®k szorzata amely i-k esetén y2i = 0 és 1 azoké, amely i-k esetén ∗ > 0) pedig az y változó y ∗ > 0 feltétel melletti feltételes eloszlását jelöli. y2i 6= 0, f (y2i |y1i 2i 1i Q
Q
0
III. típusú Tobit-modell: ∗ -ot is csak akkor gyelhetjük meg, ha Ez a modell abban különbözik a II. típustól, hogy y1i pozitív. Ez azt jelenti, hogy a magyarázó változó továbbra is függ a szelekciós változótól is, de lehet, hogy a szelekciós változó értékét sem ismerjük.
(5.19) (5.20)
∗ y1i = x01i β1 + u1i , ∗ y2i
=
x02i β2
+ u2i
( ∗, y1i y1i = 0, ( ∗, y2i y2i = 0,
∗ >0 ha y1i ∗ ≤ 0, ha y1i
i = 1, 2, ..., n
∗ >0 ha y1i ∗ ≤ 0, ha y1i
i = 1, 2, ..., n
(5.21) (5.22)
ahol u1i -k és u2i -k kétdimenziós normális eloszlásból származó független változók, melyek várható értéke 0, szórásnégyzete pedig rendre σ12 és σ22 , a kovarianciájuk pedig σ12 . A modell likelihood függvénye: L=
Y
∗ P (y1i ≤ 0)
0
Y
f (y1i , y2i ),
(5.23)
1
∗ és y ∗ együttes eloszlása. ahol f (y1i , y2i ) az y1i 2i
IV. típusú Tobit-modell:
Ez a modell a III. típustól abban különbözik, hogy gyelembe veszünk egy újabb magyarázó ∗ -t, amit csak akkor gyelhetünk meg, ha y ∗ ≤ 0. Ez azt jelenti, hogy a változót, y3i 1i szelekciós változó értékét®l függ®en gyelhetjük meg az egyik, illetve a másik magyarázó változót.
33
A modell a következ®képpen írható fel: (5.24) (5.25) (5.26)
∗ y1i = x01i β1 + u1i , ∗ y2i ∗ y3i
= =
x02i β2 x03i β3
+ u2i , + u3i
( ∗, y1i y1i = 0, ( ∗, y2i y2i = 0, ( ∗, y3i y3i = 0,
∗ >0 ha y1i ∗ ≤ 0, ha y1i
i = 1, 2, ..., n
∗ >0 ha y1i ∗ ≤ 0, ha y1i
i = 1, 2, ..., n
∗ ≤0 ha y1i ∗ > 0, ha y1i
i = 1, 2, ..., n
(5.27) (5.28) (5.29)
ahol u1i -k u2i -k és u3i -k független változók egy háromdimenziós normális eloszlásból. A modell likelihood függvénye a következ®képpen írható fel: L=
YZ 0
0
−∞
∗ ∗ f3 (y1i , y3i )dy1i
Y
f2 (y1i , y2i ),
(5.30)
1
∗ , y ) az y ∗ és y ∗ együttes eloszlása, f (y , y ) pedig y ∗ és y ∗ együttes eloszlása. ahol f3 (y1i 3i 2 1i 2i 1i 3i 1i 2i
V. típusú Tobit-modell:
Az V. típus a IV. típus átalakításával kapható: elhagyjuk az y1i -re vonatkozó egyenletet. A szelekciós változó tehát ismert, és az ® értékét®l függ, hogy melyik magyarázó változót használjuk. A modell ebben az esetben: (5.31) (5.32) (5.33)
∗ y1i = x01i β1 + u1i , ∗ y2i ∗ y3i
= =
x02i β2 x03i β3
( ∗, y2i y2i = 0, ( ∗, y3i y3i = 0,
+ u2i , + u3i ∗ >0 ha y1i ∗ ≤ 0, ha y1i
i = 1, 2, ..., n
∗ ≤0 ha y1i ∗ < 0, ha y1i
i = 1, 2, ..., n
(5.34) (5.35)
ahol u1i -k u2i -k és u3i -k független változók egy háromdimenziós normális eloszlásból. A modellhez tartozó likelihood függvény: L=
YZ 0
0
−∞
∗ ∗ f3 (y1i , y3i )dy1i
YZ 1
∞
∗ f2 (y1i , y2i )dy1i ,
0
ahol f3 és f2 de níciója ugyanaz, mint a IV. modellnél. 34
(5.36)
5.3. A sampleSelection csomag
Az R sampleSelection csomagja, melyet például Toomet és Henningsen cikkéb®l ([10]) ismerhetünk meg, éppen a fentiekben bemutatott modellezésre alkalmas. Megtalálható benne többféle Tobit-modell, de a legnépszer¶bbek talán a Heckman-féle (kétlépcs®s) eljárások, melyeket az el®z® szakaszban is említettem. A fent említett [10] cikkb®l szeretnék két példát hozni a következ®kben. El®ször a II. típusú Tobit-modell szerint készítünk becslést, melyhez mi generáljuk az adatokat. Az R-beli kód a következ®: rm(list=ls()) set.seed(0) eps <- rmvnorm(500, c(0,0), matrix(c(1,-0.7,-0.7,1), 2, 2)) xs <- runif(500) xo <- runif(500) ysx <- xs + eps[,1] > 0 yox <- xo + eps[,2] ys <- ysx yo <- yox*(ysx > 0) M <- selection(ys ~ xs, yo ~ xo) summary(M)
El®ször is kétdimenziós normális eloszlású zajokat generálunk −0.7 korrelációval. Ezután egyenletes eloszlású magyarázó változókat generálunk, xs a szelekciós változót magyarázó változó, xo pedig a meg gyelt választ magyarázó változó. ysx logikai változó, melynek értéke az eps zajjal terhelt szelekciós változó pozitív el®jelét®l függ. A szelekciós választ ys a meg gyelt választ yo jelöli. Ismert érték a szelekciós változó (ys==ysx), a szelekciós változót (ysx-t) magyarázó változó (xs), a célváltozónak (yox-nek) a szelektált értéke (yo), és a célváltozó magyarázó változója (xo). A modellt a selection függvény készíti el. Az eredmény: -------------------------------------------Tobit 2 model (sample selection model) Maximum Likelihood estimation Newton-Raphson maximisation, 5 iterations Return code 1: gradient close to zero Log-Likelihood: -712.3163 500 observations (172 censored and 328 observed) 6 free parameters (df = 494) Probit selection equation: Estimate Std. error t value Pr(> t) (Intercept) -0.2228 0.1081 -2.061 0.0393 * xs 1.3377 0.2014 6.642 3.09e-11 ***
35
Outcome equation: Estimate Std. error t value Pr(> t) (Intercept) -0.0002265 0.1294178 -0.002 0.999 xo 0.7299070 0.1635925 4.462 8.13e-06 *** Error terms: Estimate Std. error t value Pr(> t) sigma 0.9190 0.0574 16.009 < 2e-16 *** rho -0.5392 0.1521 -3.544 0.000394 *** --Signif. codes: 0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 --------------------------------------------
Második példaként az R-ben a sampleSelection csomag használatával elérhet® Mroz87 adathalmazt vizsgálom ([11]). Az adathalmaz 753 házasságban él® n®r®l tartalmaz adatokat, például a gyermekeik száma (külön 5 év alatt, illetve 6 és 18 év között), képzettségi szintjük, életkoruk, korábbi munkatapasztalat években kifejezve, a munkanélküliség aránya a lakóhelyükön, stb. A vizsgálat arra kereste a választ, hogy hogyan függ ezekt®l az adatoktól a n®k jövedelme. Az adathalmazt egy Mroz87.Rdata nev¶ fájlba mentettem, és ezzel dolgoztam. Az lfp változót (labour force participation ) szeretnénk modellezni, vagyis a munkaer®-piaci jelenlétet. Az R-beli modell: load("Mroz87.Rdata") Mroz87$kids <- ( Mroz87$kids5 + Mroz87$kids618 > 0 ) M <- heckit( lfp ~ age + I( age^2 ) + faminc + kids + educ, wage ~ exper + I( exper^2 ) + educ + city, data = Mroz87, method = "2step" ) summary(M)
A Mroz87$kids logikai változó, amely azt mondja meg, hogy van-e gyermek a családban (összeadtuk a külön 5 év alatti, illetve 6-18 év közötti gyermekek számát minden családban). Az lfp változót az életkor másodfokú függvényével (age), a család jövedelméve (faminc), a gyermekek jelenlétével, (kids), és a képzettségi szinttel (educ, években) modellezzük. Modellezzük továbbá a keresetet is (wage) a tapasztalat (exper), a képzettség (educ), illetve a nagyvárosi jelenlét (city) függvényében. A kapott eredmény: -------------------------------------------Tobit 2 model (sample selection model) 2-step Heckman / heckit estimation 753 observations (325 censored and 428 observed) 14 free parameters (df = 740)
36
Probit selection equation: Estimate Std. Error (Intercept) -4.157e+00 1.402e+00 age 1.854e-01 6.597e-02 I(age^2) -2.426e-03 7.735e-04 faminc 4.580e-06 4.206e-06 kidsTRUE -4.490e-01 1.309e-01 educ 9.818e-02 2.298e-02
t value -2.965 2.810 -3.136 1.089 -3.430 4.272
Pr(>|t|) 0.003127 0.005078 0.001780 0.276544 0.000638 2.19e-05
** ** ** *** ***
Outcome equation: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.9712003 2.0593505 -0.472 0.637 exper 0.0210610 0.0624646 0.337 0.736 I(exper^2) 0.0001371 0.0018782 0.073 0.942 educ 0.4170174 0.1002497 4.160 3.56e-05 *** city 0.4438379 0.3158984 1.405 0.160 Multiple R-Squared:0.1264, Adjusted R-Squared:0.116 Error terms:
Estimate Std. Error t value Pr(>|t|) invMillsRatio -1.098 1.266 -0.867 0.386 sigma 3.200 NA NA NA rho -0.343 NA NA NA --------------------------------------------
A 325 cenzorált változó itt azt jelenti, hogy ennyi n® nem dolgozott a meg gyeltek közül, a többi 428 pedig pozitív számú órát dolgozott, amit meg gyeltek. Az eredményekb®l az látható, hogy a munkaer®piaci jelenlét csak a család teljes jövedelmét®l (faminc) nem függ szigni kánsan, a többi tényez®t®l igen. Ugyanakkor a dolgozó n®k jövedelme a vizsgált tényez®k közül egyedül a képzettségi szintjükt®l (educ) függ szigni kánsan.
37
6. fejezet
Összegzés Szakdolgozatomban az ipari és tudományos kísérletek tervezésével, gazdasági modellek vizsgálatával és az elemzés statisztikai módszereivel foglalkoztam. A 2. fejezetben azt mutattam be, hogy a kísérletek tervezésénél mennyire fontos a véletlen szerepe. Példákat hoztam arra, hogy milyen hibákhoz vezethet, ha nem véletlenítjük a kísérletünket, és azt is láttuk, hogy hogyan lehet egy véletlen nullhipotézisen alapuló próbát elvégezni akkor, amikor a szokásos elemzési módszerek feltételei nem teljesülnek. A 3. fejezetben bemutattam a két leggyakrabban vizsgált egyfaktoros modellt, a redukált és a teljes modellt, és becsléseket adtam a modellek paramétereire. Az elemzéshez az ANOVA módszert használtam. Láttuk azt is, hogy az ANOVA önmagában nem elegend®, hiszen általában nem csak annyit szeretnénk megállapítani, hogy van-e különbség a kísérleti csoportok között, hanem azt is, hogy pontosan mely csoportok között és mekkora különbségek vannak. Ennek a kérdésnek a megválaszolásához de niáltam a kontrasztokat, és egyszer¶ példán keresztül bemutattam azt is, hogy hogyan lehet kontraszt együtthatók segítségével megadni a nullhipotézist. A 4. fejezetben azt vizsgáltam, hogy hogyan lehet az egyfaktoros esetben felírt modellt több faktorra általánosítani. Részletesebben a kétfaktoros esetr®l volt szó. Láttuk, hogy amikor több faktor befolyásolja a kísérlet eredményét, akkor nem elég csak a faktorok hatását külön vizsgálni, hiszen kölcsönhatás is gyakran fellép a különböz® faktor-szint kombinációk között. A kétfaktoros esethez képest további általánosításokról és a részleges faktoriális tervekr®l csak említésszer¶en írtam, mert ezeknek a megfelel® bemutatása már meghaladná ennek a szakdolgozatnak a kereteit. Az 5. fejezetben azt láttuk, hogy nem minden kutatási területen lehet kísérleteket tervezni. A közgazdaságtanban és a szociológiában például legtöbbször meg gyelések adataival kell dolgozni, ami tovább nehezedik, ha olyan változók befolyásolják a válasz változót, amelyeket nem is tudunk meg gyelni. A dolgozatomban bemutattam modellek egy olyan családját, amelyek korlátos függ® változókat, csonkított és cenzorált változókat, illetve szelekciós változókat használnak a válasz változók magyarázásához. Minden fejezetben láttunk R-beli példákat is, az 5. fejezetben pedig külön szakaszban írtam a sampleSelection csomagról, amely kifejezetten a bemutatott modellek számítógépes megvalósításához készült. Mindent összevetve elmondható, hogy bármilyen kutatáshoz elengedhetetlen a kísérletezés vagy adatgy¶jtés, valamint a kapott válasz értékek elemzése. A matematikai statisztika módszereinek ismeretével ezeket a lépéseket hatékonyan tudjuk elvégezni, melyhez számos statisztikai program, például az R is nagy segítséget nyújt. 38
Irodalomjegyzék [1] Gary W. Oehlert: A First Course in Design and Analysis of Experiments (2010) http://users.stat.umn.edu/∼gary/Book.html
[2] Pröhle Tamás - Zempléni András: Többdimenziós statisztika számítógépes módszerei. Elektronikus jegyzet. (2013.06.28.) http://www.cs.elte.hu/∼zempleni/tobbdim_stat.pdf
[3] Csiszár Vill®: Statisztika. Elektronikus jegyzet. (2009. május 6.) http://www.cs.elte.hu/∼villo/esti/stat.pdf
[4] Student: The Probable Error of a Mean, Biometrika, Vol. 6, No. 1 (March, 1908), pp. 1-25 [5] Lukács Ottó: Matematikai statisztika, M¶szaki Kiadó (2006) [6] Hervé Abdi - Lynne J. Williams: Contrast analysis, Encyclopedia of Research Design (Szerk.: Neil J. Salkind), SAGE Publications (2010), pp. 243-251 https://www.utdallas.edu/∼herve/abdi-contrasts2010-pretty.pdf
[7] James Tobin: Estimation of Relationships for Limited Dependent Variables, Econometrica, Vol. 26., No. 1 (Jan., 1958) pp. 24-36 [8] Takeshi Amemiya: TOBIT MODELS: A SURVEY, Journal of Econometrics 24 (1984) pp. 3-61. [9] James J. Heckman: The Common Structure of Statistical Models of Truncation, Sample Selection and Limited Dependent Variables and a Simple Estimator for Such Models, Annals of Economic and Social Measurement, Vol. 5, Number 4, Oct. 1976, pp. 475492. [10] Ott Toomet - Arne Henningsen: Sample Selection Models in R: Package sampleSelection, Journal of Statistical Software, Jul. 2008, Vol. 27, Issue 7 https://www.jstatsoft.org/article/view/v027i07
[11] Thomas A. Mroz: The Sensitivity of an Empirical Model of Married Women's Hours of Work to Economic and Statistical Assupmtions, Econometrica, Vol. 55, Issue 4 (Jul. 1987), pp. 765-799. [12] https://cran.r-project.org/
39
Köszönetnyilvánítás Ezúton is szeretném hálás köszönetemet kifejezni témavezet®mnek, Pröhle Tamás tanár úrnak, aki amellett, hogy gyelmembe ajánlotta ezt a témát számos konzultációs lehet®séggel és hasznos tanáccsal is hozzájárult ahhoz, hogy a dolgozatom elnyerje végleges formáját.
40
