Statisztika 1 - Lab 10

A koleszterines adatok, ha kellenek majd.

cholestr<-read.delim("http://www.statsci.org/data/general/cholestr.txt")
cholestg<-read.delim("http://www.statsci.org/data/general/cholestg.txt")

Előjelpróba (Sign test)

A múlt heti anyagban szereplő Wilcoxon (előjeles rang) próbának feltétele, hogy a különbség változó szimmetrikus legyen. A Sign-Rank tesztnek nincs ilyen feltétele.

Azt teszteli, hogy két összefüggő minta (\(X_1,\dots, X_n\) és \(Y_1,\dots, Y_n\)) esetén a különbségek (\(Z_i:=Y_i-X_i,\ i=1,\dots,n\)) háttérváltozójának (\(Z\)) mediánja \(0\)-e. Feltétel, hogy a \(Z_i\) különbségek legyenek függetlenek és folytonos háttérváltozójúak. Az eddigi jelölésekkel a hipotézisek: \[ H_0:\ \mathbb{P}(Z\ge 0)=0.5 \quad\text{vs.}\quad H_1:\ \mathbb{P}(Z\ge 0)\neq 0.5, \] vagy \[ H_0:\ \mathbb{P}(Y\ge X)=0.5 \quad\text{vs.}\quad H_1:\ \mathbb{P}(Y\ge X)\neq 0.5, \] A próbastatisztika és annak eloszlása \(H_0\) esetén: \[ W=\sum_{i=1}^n \mathbb{I}_{\{Y_i>X_i\}} \underset{H_0}{\sim} \text{Binom}{\left(m,0.5\right)}, \] ahol \(m:=\sum_{i=1}^n \mathbb{I}_{\{Y_i=X_i\}}\), tehát azoknak a pároknak a száma, ahol nincs különbség (\(m\le n\)). Vegyük észre, hogy a próbastatisztika csupán a szigorú egyenlőtlenségeket számolja le (valamelyik irányban).

Ha vizsgáljuk az összepárosított mintánkat, akkor azt látjuk, hogy kétfajta kimenetel lehet, mivel az egyenlő párokat nem vesszük figyelembe: \(Y_i<X_i\) vagy \(Y_i>X_i\). A nullhipotézis fennállása (kétoldali esetben) nagyjából azt jelentené, hogy egyik változó sem szignifikánsan nagyobb, mint a másik: \[ \mathbb{P}(Y_i<X_i) = \mathbb{P}(Y_i>X_i) = 0.5, \] vagy bevezetve a \(Q:=\mathbb{I}_{\{Y>X\}}\) változót, a \(\{Q_i\}_{i=1}^m\) minta (csupa 0 és 1-es) alapján akarjuk eldönteni, hogy \(Q\sim \text{Binom}(m,0.5)\) eloszlású-e. Tehát az előjel teszt az alábbi ú.n. binomiális teszt speciális változata.

Binomiális teszt

Egy két kategóriát tartalmazó f.a.e (\(\{Q_i\}_{i=1}^n\)) mintáról (elkódolva pl. \(Q\in\{a,b\}\)) szeretnénk eldönteni, hogy a \(Q\) Bernoulli eloszlású háttérváltozó paramétere a nullhipotézisben adottal egyezik-e meg. Tehát a hipotézisek: \[ H_0:\ \mathbb{P}(Q = a)=p_0 \quad\text{vs.}\quad H_1:\ \mathbb{P}(Q = a)\neq p_0, \] vagy ekvivalensen \[ H_0:\ \mathbb{P}(Q = b)=1-p_0 \quad\text{vs.}\quad H_1:\ \mathbb{P}(Q = b)\neq 1-p_0. \] Legyen \(A\) a sikerek (‘\(a\)’-k) száma \(n\) darab Bernoulli próbából, és \(k=\sum_{i=1}^n \mathbb{I}_{\{Q_i=a\}}\) a sikerek (‘\(a\)’-k) száma a mintában. Tehát \(A\) a binomiális eloszlású háttérváltozó, míg \(k\) a próbastatisztika. Akkor tudjuk, hogy \(H_0\) mellett \[ \mathbb{P}(A = k) \underset{H_0}{=} {n \choose k} p_0^k (1-p_0)^{n-k}, \quad\text{és}\quad \mathbb{E}(A)\underset{H_0}{=}n p_0. \] Ez alapján a \(k\)-ból számolt p-érték kétoldali esetben a szokásos: \[ p\text{-érték} = 2\cdot \min\{\mathbb{P}_{H_0}(A\le k), \mathbb{P}_{H_0}(A\ge k)\}, \] ahol egyébként könnyű látni, hogy a \(\min\) attól függ, hogy \(k<np_0\) vagy \(k>np_0\), vagyis, hogy a \(H_0\) szerinti várható értékhez képest hol vagyunk. Ez a \(p\)-érték a binomiális eloszlásra egzakt kiszámolható, de nehezen, túl számításigényes. Magas elemszám esetén viszont közelíthető különböző folytonos eloszlású statisztikákkal.

A fentiek alapján könnyen látható, hogy az Előjelpróba valójában a különbségváltozó előjelére (pozitív vagy negatív) vett Binomiális próba, ahol \(p_0=0.5\) a nullhipotézisbeli paraméter.

R példa

Sajnos alapértelmezetten nincs beépített fv. az előjelpróbára, de a Binomiális teszt elérhető, így írhatunk köré egy fv.-t, ami az Előjelpróbát végzi el.

sign.test<-function(x, y, data = NULL, ...){
    # x : első minta oszlop/vektorként vagy csak oszlopnévként
    # y : második minta oszlop/vektorként vagy csak oszlopnévként
    # data : (opcionális) data frame, amelynek oszlopa x és y
    # ... : egyéb paraméterek, amik továbbadódnak a binomiális tesztnek
    
    # csak a nemegyenlő és nem hiányzó részhalmazt vizsgáljuk
    # és leszámoljuk hányszor nagyobb y, mint x
    if(is.null(data)){
        noteq<-sum(x != y, na.rm = TRUE)
        bigger<-sum(x < y, na.rm = TRUE)
    } else {
        noteq<-sum(data[[x]] != data[[y]], na.rm = TRUE)
        bigger<-sum(data[[x]] < data[[y]], na.rm = TRUE)
    }
    
    # binomiális tesztet végzünk
    # az első paraméter a  'sikerek' száma
    # a második paraméter a (nemegyenlő) elemek száma
    # a harmadik ... paraméter helyére a főfüggvény hivásakor
    # az ottani ... helyére írt paraméterek fognak tábbadódni
    bobj<-binom.test(bigger, noteq, p=0.5, ...)
    
    # a teszt eredményeit összefoglaló objektumban 
    # a binom teszthez specifikus stringeket átírhatjuk
    bobj$method<-"Sign test"
    names(bobj$statistic)<- "number of times y is greater than x"
    names(bobj$parameter)<- "number of non-equal pairs"
    names(bobj$estimate)<- "probability of y is greater than x"
    names(bobj$null.value)<- "value of Pr(x<y)"
    bobj$data.name<- "paired sample of x and y"
    
    # visszaadjuk a módosított objektumot 
    return(bobj)
}

Megint nézzük meg a betegek második és negyedik napi koleszterinszintjét. Teszteljük, hogy a különbség mediánja 0-e. Alább láthatjuk, hogy mindkétféle megadással működik a függvényünk, és a kimenetünk is az Előjelpróbának megfelelőre íródott át.

sign.test('Day2', 'Day4', cholestr)

## 
##  Sign test
## 
## data:  paired sample of x and y
## number of times y is greater than x = 9, number of non-equal pairs =
## 28, p-value = 0.08716
## alternative hypothesis: true value of Pr(x<y) is not equal to 0.5
## 95 percent confidence interval:
##  0.1587760 0.5235164
## sample estimates:
## probability of y is greater than x 
##                          0.3214286

sign.test(cholestr$Day2, cholestr$Day4)

## 
##  Sign test
## 
## data:  paired sample of x and y
## number of times y is greater than x = 9, number of non-equal pairs =
## 28, p-value = 0.08716
## alternative hypothesis: true value of Pr(x<y) is not equal to 0.5
## 95 percent confidence interval:
##  0.1587760 0.5235164
## sample estimates:
## probability of y is greater than x 
##                          0.3214286

A p-érték alapján 5%-os szignifikancia mellett nem utasíthatjuk el a nullhipotézist, tehát a különbség (y-x) mediánját elfogadhatjuk 0-nak, de pl. 10%-os szignifikancia mellett elutasíthatjuk a nullhipotézist.

Mivel a beépített Binomiális tesztre van beépített egyoldali próba az is működni fog a megírt fv.-ünkkel. Azt szeretnénk tesztelni, hogy ahogy a napok telnek csökkenni kezd a koleszterinszint, tehát hozzáadjuk az alternative="less" paramétert.

sign.test('Day2', 'Day4', cholestr, alternative="less")

## 
##  Sign test
## 
## data:  paired sample of x and y
## number of times y is greater than x = 9, number of non-equal pairs =
## 28, p-value = 0.04358
## alternative hypothesis: true value of Pr(x<y) is less than 0.5
## 95 percent confidence interval:
##  0.000000 0.493789
## sample estimates:
## probability of y is greater than x 
##                          0.3214286

Itt 5%-os szignifikancia mellett elutasíthatjuk a nullhipotézist, hogy \(\mathbb{P}(kol_2<kol_4)\ge 0.5\). Tehát annak a valószínűsége, hogy nőtt az érték, kisebb, mint 0.5; vagyis annak a valószínűsége, hogy csökkent az érték, nagyobb, mint 0.5.

Megj.: az iménti következtetést 1%-os szignifikancia mellett már nem tudnánk megtenni.

Vizsgáljuk most a második napit a 14. napihoz képest. Itt már kisebb lesz az elemszám, sokan nem lettek a 14. napon megmérve. Fv.-ünk szerencsére ezt is kezeli a megfelelő helyen beszúrt na.rm=TRUE paraméterek miatt.

sign.test('Day2', 'Day14', cholestr, alternative="less")

## 
##  Sign test
## 
## data:  paired sample of x and y
## number of times y is greater than x = 4, number of non-equal pairs =
## 19, p-value = 0.009605
## alternative hypothesis: true value of Pr(x<y) is less than 0.5
## 95 percent confidence interval:
##  0.0000000 0.4191202
## sample estimates:
## probability of y is greater than x 
##                          0.2105263

Itt már 1%-os szignifikancia mellett is elutasíthatjuk a nullhipotézist, hogy \(\mathbb{P}(kol_2<kol_{14})\ge 0.5\). Tehát annak a valószínűsége, hogy nőtt az érték kisebb, mint 0.5; vagyis annak a valószínűsége, hogy csökkent az érték nagyobb, mint 0.5.

Megj.: Az Előjelpróbának kevesebb feltétele van, mint a Wilcoxon előjeles rang próbának, de utóbbinak az ereje (\(H_{alt}\) mellett \(H_{alt}\) elfogadásának valószínűsége) nagyobb és kevesebb adatra van szüksége.

Illeszkedésvizsgálat, Normalitás tesztek

Egymintás Kolmogorov–Szmirnov teszt

Az egymintás Kolmogorov–Szmirnov teszttel folytonos háttérváltozók eloszlásának (\(F_X\)) illeszkedését lehet vizsgálni egy a nullhipotézisben megadott eloszláshoz (\(F_0\)). A nullhipotézis tehát az, hogy a megadott eloszlást követi a minta: \[ H_0:\ F_X(x) = F_0(x) \qquad\text{vs.}\qquad H_1:\ F_X(x) \neq F_0(x) \] A próbastatisztikához szükség van a már tanult tapasztalati eloszlásfv.-re. Egy \(X_1, X_2,\dots, X_n\) f.a.e. minta esetén ez: \[ F_X(x)= \mathbb{P}\left( X\le x \right) \approx \frac{\left| \left\{ i: X_i\le x \right\}\right|}{n} =: F_n(x) \qquad \forall x\in \mathbb{R} \] A próbastatisztika pedig: \[ S_n =\sqrt{n}\ D_n = \sqrt{n}\ \sup_x \left|F_n(x)-F_0(x)\right|, \] amelynek az eloszlása \(H_0\) fennállása esetén aszimptotikusan (ahogy \(n\to\infty\)) Kolmogorov-eloszlású.

Több statisztikai programban a Kolmogorov-Szmirnov tesztnél van egy kis csalás. Csak azt kell megadni, hogy milyen eloszláshoz való illeszkedést teszteljen a program. Ezután például normálishoz való illeszkedés esetén a program kiszámolja az empirikus várható értéket és szórást, és az ilyen paraméterű normálishoz hasonlítja az adatokat. Mindezt úgy lehet megfogalmazni, hogy a program egy kicsit módosított Kolmogorov–Szmirnov próbát hajt végre, ami hajlik arra, hogy elfogadjuk a nullhipotézist. Ennek az előnye az, hogy általában lehet tesztelni azt, hogy a változó például normális-e, anélkül, hogy ismernénk a várható értéket és a szórását, hátránya a leírt torzítás. Léteznek az említett torzítást korrigáló statisztikai eljárások.

Végezzünk el egy normalitástesztet (azaz a normális eloszlás kerül a nullhipotézisbe) a 2. napi koleszterinszinten. Először nézzük meg a tapasztalati eloszlásfv. képét. Ez R-ben sokkal egyszerűbb, mint Excelben volt.

Fn<-ecdf(cholestr$Day2)
plot(Fn)

A teszthez elegendő megadni a tesztelendő eloszlás nevét és paramétereit. Ha paramétereket nem adunk meg, akkor azok alapértelmezett értékével tesztel, azaz pl. normális eloszlás esetén (0,1) paraméterekkel.

ks.test(cholestr$Day2, "pnorm", 260, 50)

## Warning in ks.test(cholestr$Day2, "pnorm", 260, 50): ties should not be present
## for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  cholestr$Day2
## D = 0.14488, p-value = 0.5992
## alternative hypothesis: two-sided

Most próbáljuk meg mintából becsülni az eloszlás paramétereit. Ez a korábbi megjegyzés miatt csalás, és a fv. leírása is kiemeli, hogy a K– S teszt becsült paraméteres változatai nincsenek itt implementálva.

ks.test(cholestr$Day2, "pnorm", mean(cholestr$Day2), sd(cholestr$Day2))

## Warning in ks.test(cholestr$Day2, "pnorm", mean(cholestr$Day2), sd(cholestr
## $Day2)): ties should not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  cholestr$Day2
## D = 0.13559, p-value = 0.682
## alternative hypothesis: two-sided

Mindkét esetben elfogadhatjuk, hogy az 2. napi koleszterinszint normális eloszlású.

Ennek a tesztnek rengeteg módosítása van: * Kétmintás: annak tesztelésére, hogy két minta háttérváltozói ugyanazt az ismeretlen eloszlást követik-e * Egyoldali: sztochasztikus dominancia tesztelésére. * Becsült paraméteres változatok: néhány ismertebb eloszlásra létezik a tesztnek olyan változata, amely korrekt maradt becsült paraméterek használata mellett is. * Nem mindenhol folytonos eloszlások tesztelésére készített változat.

Végezetül ábrázoljuk egymáson a tapasztalati eloszlásfv.-t és a normális eloszlás becsült paraméteres változatát.

kol2 <- cholestr$Day2
Fn<-ecdf(kol2)
xx<- seq(min(kol2), max(kol2), length.out = 100)
plot(xx, pnorm(q=xx,
               mean = mean(kol2),
               sd = sd(kol2)),
     type='l', col="blue", xlab = 'x', ylab = "F(x)")
plot(Fn, add=TRUE)

Shapiro–Wilk teszt

A Shapiro–Wilk teszt az előzővel ellentétben kifejezetten normális eloszlás tesztelésére szolgál. A teszt elméleti része túlmutat a kurzus keretein, de a Wikipédián érdemes ránézni még.

Érdemes még megjegyezni, hogy a teszt paraméterek megadása nélkül teszteli a normalitást. Tehát nem egy adott paraméterű normális eloszláshoz való illeszkedést vizsgál, hanem általánosabban nézi a normalitást.

shapiro.test(cholestr$Day2)

## 
##  Shapiro-Wilk normality test
## 
## data:  cholestr$Day2
## W = 0.96894, p-value = 0.5525

Ez alapján is elfogadhatjuk, hogy a 2. napi koleszterinszint normális eloszlású.

Megj.: Szimulációk azt mutatják, hogy népszerűbb normalitástesztek közül ennek az ereje a legnagyobb.

Khí-négyzet (\(\chi^2\)) próba (Chi-squared test)

Többféle khí-négyzet próba létezik, de mikor csak annyit mondunk, hogy khí-négyzet próba, akkor a Pearson-féle khí-négyzet próbára utalunk.

Pearson-féle khí-négyzet próba

A Pearson-féle khí-négyzet próbát kategorikus adathalmazokra alkalmazzuk és azt tudjuk megvizsgálni, hogy a megfigyelt halmazok közötti különbség betudható-e a véletlennek. Pontosabban, a Pearson-féle khí-négyzet próbával háromféle összehasonlítást végezhetünk: illeszkedésvizsgálat (goodness of fit), homogenitás (homogeneity), és függetlenség (independence).

• Az illeszkedésvizsgálat során azt vizsgáljuk, hogy különbözik-e a megfigyelt gyakoriság egy elméleti eloszlástól. Például a megfigyelt dobókocka szabályos vagy cinkelt?
• A homogenitásvizsgálat kettő vagy több csoport/populáció ugyanannak (kategorikus) tulajdonságának/változójának eloszlását hasonlítja össze. Például: különbözik-e a férfiak és nők TV nézési szokásai.
• A függetlenségvizsgálat esetén azt vizsgáljuk, hogy két változó megfigyelései, kontingenciatáblázatban kifejezve, függetlenek-e egymástól. Például van-e kapcsolat a nemek (férfi v. nő) és empátia között (alacsony v. magas).

Hipotézisek

$H_0 = $ A megfigyelt gyakoriságok és az elméleti eloszlás között nincs különbség
$H_1 = $ Van különbség az eloszlások között

Próbastatisztika általánosabban

A próbastatisztika lényegében a megfigyelt relatív gyakoriságok és elméleti értékek (valószínűségek) normalizált négyzetes eltérésének összege (normalized sum of squared deviations), aminek értéke a következő:

\[ \chi^2 = \sum_{i=1}^r \frac{(O_i - E_i)^2}{E_i}, \] ahol

• \(\chi^2 =\) Pearson-féle kumulatív próbastatisztika, ami a \(H_0\) fennállása esetén aszimptotikusan \(\chi^2\) eloszlású.
• \(O_i =\) A megfigyelések száma az \(i\) csoportban.
• \(E_i =\) A várható (elméleti) száma az \(i\) csoportbeli elemeknek.
• \(r =\) a “táblázat” celláinak száma, a lehetséges értékek

Ebből vezethető le a háromféle fenti összehasonlításnak megfelelő egy-egy próbastatisztika.

Próba számítása

Mindhárom próbához/teszthez a próba számítása a következő lépésekből áll:

1. Számoljuk ki a \(\chi^2\) próbastatisztikát
2. Határozzuk meg a szabadságfokot:
   • Illeszkedésvizsgálat esetén a szabadságfok \(df = k - 1\), ahol \(k\) a megfigyelt változó értékkészletének mérete. (Ha az elméleti valószínűségeket is becsüljük, akkor a becsült paraméterek számával csökken a szabadsági fok)
   • Homogenitásvizsgálat esetén a szabadságfok \(df = (rows - 1) \times (cols - 1)\), ahol \(rows\) a kontingencia táblázat sorainak (kategóriák), a \(cols\) az oszlopainak (csoportok) számát jelöli.
   • Függetlenségvizsgálat esetén hasonlóan \(df = (rows - 1) \times (cols - 1)\), ahol \(rows\) az egyik változó értékeinek száma, \(cols\) a másik változó értékeinek száma.
3. Határozzuk meg a konfidenciaszintet.
4. Hasonlítsuk össze a \(\chi^2\) próbastatisztikát a \(\chi^2\)-eloszlástáblázatban található kritikus értékkel (a megfelelő szabadságfoknál és konfidenciaszintnél). A próbastatisztika nagyobb-e mint a kritikus érték?
5. Ha a próbastatisztika értéke nagyobb mint a kritikus érték, akkor elutasíthatjuk a \(H_0\) nullhipotézist a választott konfidenciaszinten.

R példák

A khí-négyzet próbák elvégzéséhez a chisq.test függvényt fogjuk használni.

Illeszkedésvizsgálat

• \(H_0:\) A megfigyelt gyakoriságok és az adott elméleti eloszlás (\(\{p_i\}_{i=1}^r\)) között nincs különbség
• \(H_1:\) Van különbség a gyakoriság és eloszlások között. vagy:
• \(H_0:\) A minta háttérváltozójának eloszlása az adott \(\{p_i\}_{i=1}^r\)
• \(H_1:\) Az eloszlása nem \(\{p_i\}_{i=1}^r\)

A próbastatisztika itt:

\[ \chi^2 = \sum_{i=1}^r \frac{(O_i - n p_i)^2}{n p_i} = n \sum_{i=1}^r \frac{(O_i/n - p_i)^2}{p_i}, \] ahol

• \(n =\) Az összes megfigyelés száma.
• \(O_i =\) A megfigyelések száma az \(i\) csoportban.
• \(p_i =\) Az \(i\) csoport aránya (valószínűsége) a populációban \(H_0\) esetén.
• \(r =\) A tesztelt változó lehetséges értékeinek száma.

Nézzünk egy példát.

Tegyük fel, hogy szedtünk tulipánokat, és összesen 81 piros, 50 sárga és 27 fehér tulipánt gyűjtöttünk.

1. kérdés: Ugyanolyan az gyakorisága/eloszlása a különböző színű tulipánoknak?

Ha ugyanolyan gyakoriak lennének, akkor a várható arányuk 1/3-1/3 lenne.

• \(H_0\): Nincs szignifikáns különbség a színek megfigyelt és az elméleti/várható értékük között
• \(H_1\): Szignifikánsan különbözik a megfigyelt és az elméleti érték.

Ebben az esetben a chisq.test függvénynek a következő inputokat kell megadnunk:

• x: numerikus vektor
• p: valószínűségeket tartalmazó vektor, ami ugyanolyan hosszú mint \(x\)

tulip <- c(81, 50, 27) # megfigyelt gyakoriságok

result <- chisq.test(tulip, p=c(1/3, 1/3, 1/3))

result

## 
##  Chi-squared test for given probabilities
## 
## data:  tulip
## X-squared = 27.886, df = 2, p-value = 8.803e-07

Látjuk, hogy a \(p\)-érték \(8.803 \cdot 10^{-7}\), ami azt jelenti hogy a nullhipotézist nemhogy \(\alpha = 0.05\) szignifikanciaszinten de még \(\alpha = 0.001\) szignifikanciaszinten is elutasítjuk, azaz nem igaz, hogy egyenletes az eloszlása a tulipán színének.

# A nullhipotézis szerint a várhatóértéke a megfigyelt színeknek:
result$expected

## [1] 52.66667 52.66667 52.66667

names(result)

## [1] "statistic" "parameter" "p.value"   "method"    "data.name" "observed" 
## [7] "expected"  "residuals" "stdres"

2. kérdés: Tegyük fel, hogy abban a régióban, ahol szedtük a tulipánokat ott a piros, sárga és fehér tulipánok aránya 3:2:1, azaz a várható gyakoriságok:

• 3/6 (=1/2) a pirosé
• 2/6 (=1/3) a sárgáé
• 1/6 a fehéré

result <- chisq.test(tulip, p=c(1/2, 1/3, 1/6))

result

## 
##  Chi-squared test for given probabilities
## 
## data:  tulip
## X-squared = 0.20253, df = 2, p-value = 0.9037

Látjuk hogy a \(p\)-érték 0.9037, ami nagyobb mint az általában használt szignifikanciaszint \(\alpha=0.05\). Tehát a megfigyelt arányok nem különböznek szignifikánsan a várható arányoktól.

Függetlenségvizsgálat

A próbastatisztika itt, ha a két változó valódi eloszlásait külön-külön (a marginálisokat) nem ismerjük, hanem a mintából becsüljük:

\[ \begin{aligned} \chi^2 &= \sum_{i=1}^{q\cdot r} \frac{(O_i - E_i)^2}{E_i} = \sum_{i=1}^{q} \sum_{j=1}^{r} \frac{(O_{ij} - E_{ij})^2}{E_{ij}}\\ &=\sum_{i=1}^{q} \sum_{j=1}^{r}\frac{\left(\nu_{ij} - n\ \frac{\nu_{i\cdot}}{n}\frac{\nu_{\cdot j}}{n}\right)^2}{n \frac{\nu_{i\cdot}}{n}\frac{\nu_{\cdot j}}{n}}, \end{aligned} \] ahol

• \(v_{ij}\) a kontingenciatáblázat \(i\)-edik sorának \(j\)-edik eleme, azaz azon megfigyelések száma, akik az első változó szerint az \(i\) csoportba, a második változó szerint a \(j\) csoportba tartoznak
• \(n\) az összes megfigyelés száma
• \(q\) és \(r\) az első és a második változó értékkészletének mérete.
• \(v_{i \cdot}\) és \(v_{\cdot j}\) a marginálisokat jelölik, azaz \(v_{i \cdot}\) az \(i\)-edik sorösszeg, azaz \(v_{i \cdot} = \sum_{j=1}^r v_{ij}\), továbbá \(v_{\cdot j}\) a \(j\)-edik oszlopösszeg, azaz \(v_{\cdot j} = \sum_{i=1}^q v_{ij}\)

Hipotézisek

• \(H_0\): a két változó független
• \(H_1\): a két változó nem független

precízebben, ha az első változó értékkészlete \(A_1, \ldots, A_q\) és a második változó értékkészlete \(B_1, \ldots, B_r\), akkor

• \(H_0\): \(\mathbb{P}(A_i \text{ és } B_j) = \mathbb{P}(A_i)\cdot \mathbb{P}(B_j)\) \(\quad \quad \quad \quad (1 \leq i \leq q\) és \(1\leq j \leq r)\)
• \(H_1\): van olyan \(i\) és \(j\), hogy \(\mathbb{P}(A_i \text{ és } B_j) \neq \mathbb{P}(A_i)\cdot \mathbb{P}(B_j)\)

A függetlenségvizsgálatot végezzük el a híres Iris flower dataset-en. Ebben az adathalmaz három növény fajról (setosa, virginica, versicolor) tartalmaz 4-féle mérési eredményeket (szirmok hossza, széle).

Adatok betöltése:

library(datasets)

data(iris)

summary(iris)

##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 

A summary függvénnyel kérhetünk le leíróstatisztikákat az adatokról.

summary(iris$Sepal.Length)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   4.300   5.100   5.800   5.843   6.400   7.900

Nézzük meg az adatokat a View függvénnyel:

View(iris)

Vagy ki is írhatjuk néhány sorát az adatoknak:

iris[1:10, ]

##    Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1           5.1         3.5          1.4         0.2  setosa
## 2           4.9         3.0          1.4         0.2  setosa
## 3           4.7         3.2          1.3         0.2  setosa
## 4           4.6         3.1          1.5         0.2  setosa
## 5           5.0         3.6          1.4         0.2  setosa
## 6           5.4         3.9          1.7         0.4  setosa
## 7           4.6         3.4          1.4         0.3  setosa
## 8           5.0         3.4          1.5         0.2  setosa
## 9           4.4         2.9          1.4         0.2  setosa
## 10          4.9         3.1          1.5         0.1  setosa

Csak egy kategorikus változó van, a Species, de a \(\chi^2\)-próba két kategorikus változót kell megadni, így hozzunk létre még egy kategorikus változót, ami azt mutatja majd, hogy a Sepal.Length nagy vagy kicsi, pontosabban ha ez a hosszúság kisebb mint a medián, akkor legyen az értéke small, különben legyen az értéke big:

#Először mentsük el az iris adatokat egy másik változóba, hogy ne az eredetit módosítsuk:

dat <- iris

dat$size <- ifelse(dat$Sepal.Length < median(dat$Sepal.Length), "small", "big")

dat[1:5, ]

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species  size
## 1          5.1         3.5          1.4         0.2  setosa small
## 2          4.9         3.0          1.4         0.2  setosa small
## 3          4.7         3.2          1.3         0.2  setosa small
## 4          4.6         3.1          1.5         0.2  setosa small
## 5          5.0         3.6          1.4         0.2  setosa small

Most készítsünk egy kontingenciatáblát a Species és size változókat használva, ezt a table függvény segítségével tudjuk megtenni:

con_table <- table(dat$Species, dat$size)

A táblázatból látjuk, hogy például csak 1 nagy és 49 kicsi setosa virág van. Gyakorlásképpen nézzük meg ezeket az adatokat egy oszlopdiagramon:

library(ggplot2)

ggplot(dat) +
  aes(x = Species, fill = size) +
  geom_bar() +
  scale_fill_hue() +
  theme_minimal()

Most a \(\chi^2\)-próba segítségével nézzük meg, hogy van e összefüggés a faj (Species) és a méret (size) változók között:

result <- chisq.test(con_table)

result

## 
##  Pearson's Chi-squared test
## 
## data:  con_table
## X-squared = 86.035, df = 2, p-value < 2.2e-16

Már az oszlopdiagramból is egyből látszott, de ahogy a \(p\)-value is mutatja, elutasítjuk a nullhipotézist, még \(\alpha=0.0001\) szignifikanciaszinten is, azaz összefüggnek ezek a változók.

Homogenitásvizsgálat

A homogenitásvizsgálat tekinthető a függetlenségvizsgálat egy speciális változatának.

Két független (akár különböző méretű: \(n,m\)) mintában vizsgáljuk ugyanazt a változót. A két minta egy másik változó szempontjából lehet akár teljesen különböző.

• \(H_0:\) A két mintát ugyanabból az eloszlásból (\(\{p_i\}_{i=1}^r\)) mintavételeztük.
• \(H_1:\) Különbözik a háttérváltozók eloszlása.

Pongyolán az is mondható, hogy az van tesztelve, hogy a vizsgált változó független-e attól, hogy a két minta közül melyikben nézzük. A próbastatisztika itt megkapható a függetlenséget (becsült paraméterekkel) tesztelő változatból hosszas átalakítgatásokkal. Érdemes beazonosítania lenti képletben, mi minek felel meg a függetlenséget tesztelő változatban.

\[ \begin{aligned} \chi^2 &= \sum_{i=1}^{r\cdot2} \frac{(O_i - E_i)^2}{E_i} = \sum_{i=1}^r \left[ \frac{(O_{i1} - E_{i1})^2}{E_{i1}} + \frac{(O_{i2} - E_{i2})^2}{E_{i2}} \right] \\ &= \sum_{i=1}^r \left[ \frac{\left(\alpha_i - (n+m)\frac{\alpha_i+\beta_i}{n+m}\frac{n}{n+m}\right)^2}{(n+m)\frac{\alpha_i+\beta_i}{n+m}\frac{n}{n+m}} + \frac{\left(\beta_i - (n+m)\frac{\beta_i+\alpha_i}{n+m}\frac{m}{n+m}\right)^2}{(n+m)\frac{\beta_i+\alpha_i}{n+m}\frac{m}{n+m}} \right] \\ &= nm \sum_{i=1}^r \frac{(\alpha_i/n - \beta_i/m)^2}{\alpha_i+\beta_i}, \end{aligned} \] ahol

• \(n\) az első mintában, míg \(m\) a második mintában a megfigyelések száma.
• \(O_{ij} =\) A megfigyelések száma \(j\) minta \(i\) csoportjában.
• \(E_{ij} =\) A megfigyelések várható értéke \(j\) minta \(i\) csoportjában.
• \(\alpha_i =\) az első mintában, míg \(\beta_i\) a megfigyelések száma az \(i\) csoportra.
• \(r =\) A tesztelt változó lehetséges értékeinek száma.

Nézzünk egy példát.

Tegyük fel, hogy megkérdezünk \(n_1 = 240\) nőt (1. csoport) és \(n_2 = 170\) férfit (2. csoport), hogy hogyan készítik el a vasárnapi tojásos reggelit: tükörtojás (sunny side up egg), kétoldalas tükörtojás (over easy egg) vagy rántotta (scrambled eggs). A megfigyeléseket az alábbi táblázat tartalmazza.

Kérdés: Különbözik-e a tojás preferenciájuk a nőknek és a férfiaknak?

egg_preferences <- matrix(c(105,93,63,31,72,46), ncol=3)

colnames(egg_preferences) <- c("Sunny","Over Easy","Scrambled")

rownames(egg_preferences) <- c("Females","Males")


egg_preferences

##         Sunny Over Easy Scrambled
## Females   105        63        72
## Males      93        31        46

Homogenitásvizsgálatnál egy ilyen jellegű táblázatot kell megadni a chisq.test függvénynek.

result <- chisq.test(egg_preferences)

result

## 
##  Pearson's Chi-squared test
## 
## data:  egg_preferences
## X-squared = 5.5606, df = 2, p-value = 0.06202

Látjuk hogy \(p\)-érték 0.06202, ami nem kisebb mint az általában használt \(\alpha=0.05\) szignifikancia szint, így elfogadjuk a \(H_0\)-t azaz a férfiak és nők tojás sütési preferenciái nem különböznek szignifikánsan.

Appendix: A binom.test p-értéke"

A binomiális teszt \(p\)-értékére a következő “anomália” léphet fel: a kétoldali esetre az érték nem lesz duplája az egyoldali esetnek.

Egy oldali eset:

binom.test(4,22,0.075, alternative = "greater")

## 
##  Exact binomial test
## 
## data:  4 and 22
## number of successes = 4, number of trials = 22, p-value = 0.07814
## alternative hypothesis: true probability of success is greater than 0.075
## 95 percent confidence interval:
##  0.06459562 1.00000000
## sample estimates:
## probability of success 
##              0.1818182

varhert <- 22*0.075
varhert

## [1] 1.65

A tesztelt 4 a \(H_0\)-beli várható érték felett van. Próbáljuk meg reprodukálni a kapott \(p\)-értéket. Az alábbiak közül a második értéket kapjuk vissza p-értékként, ahogy várjuk:

pbinom(4,22, 0.075, lower.tail = TRUE)  # <=4

## [1] 0.9787511

pbinom(3,22, 0.075, lower.tail = FALSE) # 3< vagy 4<=

## [1] 0.0781356

pbinom(4,22, 0.075, lower.tail = FALSE) # 4< vagy 5<=

## [1] 0.02124892

Két-oldali eset:

binom.test(4,22,0.075, alternative = "two.sided")

## 
##  Exact binomial test
## 
## data:  4 and 22
## number of successes = 4, number of trials = 22, p-value = 0.07814
## alternative hypothesis: true probability of success is not equal to 0.075
## 95 percent confidence interval:
##  0.0518673 0.4028458
## sample estimates:
## probability of success 
##              0.1818182

Azt látjuk, hogy a kétoldali \(p\)-érték nem lett a duplája az előzőnek, sőt ugyanazt kaptuk vissza.

Miért nem duplája fenn a p-érték?

A binomiális eloszlás nem szimmetrikus (a várható értéke körül), ha a \(p\) paramétere nem \(0.5\). A tanult duplázós szabály szimmetrikus eloszlásokra igaz. Aszimmetrikus eloszlásokra nincs kőbevésett mód a kétoldali p-érték számolására.

Sajnos a binom.test fv. dokumentációja nem írja konkrétan, milyen módszerrel számolja a két-oldali p-értéket. Interneten utána lehet nézni a citált források alapján vagy a forráskódból kihámozni.

Egy lehetséges módszer lehet pl. a \(H_0\)-beli várható értékhez (példánkban \(1.65\)) képesti ellentétes oldali kilengést nézni. Példánkban ez a kilengés \(4-1.65=2.35\), így egy 2-oldali p-érték lehet itt: \[ p_{2,a}:=\mathbb{P}(X>1.65+2.35) + \mathbb{P}(X<1.65-2.35) = \mathbb{P}(X>4) + \mathbb{P}(X<-0.7)= \mathbb{P}(X>4) +0 \] Tehát pl. ezzel a módszerrel az egyoldali p-értéket kapjuk vissza \(4\)-re, mert a másik irányba ekkora kilengés nem lehet, hiszen a Binom eloszlás nemnegatív. Természetesen lehet, hogy az R fv.-e valamilyen teljesen más módszert használ, ami ugyanezt eredményezte a példánkban. Mindegyik módszernek megvannak a maga erősségei, gyengeségei többnyire.

Binom-eloszl. különböző \(p\) paraméterekre:

x<-seq(0,22,length.out = 23)
plot(x, dbinom(x,22,.5))

plot(x, dbinom(x,22,.75))

plot(x, dbinom(x,22,.075))

Fenn a harmadik ábrán láthatjuk, hogy a duplázás életszerűtlen lenne, hiszen a másik irányba a kilengés nem is lenne lehetséges.

Nagyobb elemszámokra hasonló eredmény tesztelése:

Legyen most 10-szer ekkora a minta és a sikerek száma.

binom.test(40,220,0.075, alternative = "greater")

## 
##  Exact binomial test
## 
## data:  40 and 220
## number of successes = 40, number of trials = 220, p-value = 1.736e-07
## alternative hypothesis: true probability of success is greater than 0.075
## 95 percent confidence interval:
##  0.1401849 1.0000000
## sample estimates:
## probability of success 
##              0.1818182

binom.test(40,220,0.075, alternative = "two.sided")

## 
##  Exact binomial test
## 
## data:  40 and 220
## number of successes = 40, number of trials = 220, p-value = 2.092e-07
## alternative hypothesis: true probability of success is not equal to 0.075
## 95 percent confidence interval:
##  0.1331849 0.2392463
## sample estimates:
## probability of success 
##              0.1818182

Ugyan itt sincs duplázás, de már van különbség a két p-érték között.

Adjunk a fenn vázolt módszer segítségével egy kétoldali \(p\)-érték szerű dolgot.

jobbo<-pbinom(39,220, 0.075, lower.tail = FALSE)

diff<- 40-220*0.075

balo<-pbinom(220*0.075-diff,220, 0.075, lower.tail = TRUE)

ketp<-jobbo+balo

A fenn vázolt módszer ebben az extrém példában is pont az egyoldali értéket adná vissza kétoldali esetben is, tehát annyit már tudunk, hogy a beépített fv. nem ezt használja. Érdemes ennek utánézni, hogy akkor mit számolhat valójában.

Próbáljuk meg kihámozni a forráskódból. (Az F1 nyomvatartása mellett kattintsunk a fv. nevére a forráskód eléréséhez.) Ott ez szerepel:

# PVAL <- switch(alternative,
#                less = pbinom(x, n, p),
#                greater = pbinom(x - 1, n, p, lower.tail = FALSE),
#                two.sided = {
#     if (p == 0) (x == 0) else if (p == 1) (x == n) else {
#         relErr <- 1 + 1e-07
#         d <- dbinom(x, n, p)
#         m <- n * p
#         if (x == m) 1 else if (x < m) {
#             i <- seq.int(from = ceiling(m), to = n)
#             y <- sum(dbinom(i, n, p) <= d * relErr)
#             pbinom(x, n, p) + pbinom(n - y, n, p, lower.tail = FALSE)
#         } else {
#             i <- seq.int(from = 0, to = floor(m))
#             y <- sum(dbinom(i, n, p) <= d * relErr)
#             pbinom(y - 1, n, p) + pbinom(x - 1, n, p, lower.tail = FALSE)
#         }
#     }
# })

Ez alapján is láthatjuk, hogy az egyoldali eset az, amit várunk, a kétoldali eset a trükkös, amit a nagy else ág ír le. Ez alapján a fenti példában (\(x=4, n=22, p_0=0.075\)): \[ p_{2,b}:=\mathbb{P}(X>4) + \mathbb{P}\left( X< \sum_{i=0}^{\lfloor1.65\rfloor} \mathbb{I}_{\left\{ \mathbb{P}(X=i) \le \mathbb{P}(X=4)\right\}}-1 \right) \] Tehát nagyjából az egyoldali \(p\)-értékhez (\(\mathbb{P}(X>4)\)) azokat a valószínűségeket adta még hozzá a \(H_0\)-beli várható érték másik oldaláról, amelyek a megadott értéknél is kevésbé valószínűek, ha van ilyen, kivéve a legutolsót. Ha a megadott érték a \(H_0\)-beli várható érték másik oldalára esett volna, a fentivel analóg módon történne a számolás. A dokumentációban citált irodalmak valamelyikében vállalkozó kedvűek utánajárhatnak, kihez fűződik ez a módszer.

Binom-eloszl. nagyobb elemszámmal:

A binomiális eloszlás CHT miatt, ha \(n\) elég nagy, akkor elég közel lesz egy normális eloszláshoz, ami már szimmetrikus, így a duplázós szabály nagy \(n\)-re már elég jó közelítés lehet a p-értékre, ha \(p\) paraméter sem extrém kicsi vagy nagy.

Pl. már 220-as elemszámra is elég normálisnak tűnik az eloszlás:

x<-seq(0,220,length.out = 221)
plot(x, dbinom(x,220,.5))

plot(x, dbinom(x,220,.75))

plot(x, dbinom(x,220,.075))

Az asszimmetrikus eloszlás problémája nem okozott gondot az Előjelpróbánál, mert ott a háttérben \(p=0.5\)-re történik a binomiális teszt, amire szimmetrikus a binomiális eloszlás is.