Testy četností

Porovnáváme obdržené počty v diskrétních kategoriích s očekávanými hodnotami.

Použití testu \({\chi ^2}\) je jen aproximací problému. Funguje dobře, pokud jsou očekávané četnosti větší nebo rovny 5 v 80 % kategorií a žádná nemá očekávaný počet 0. Jinak je dobré kategorie s malým zastoupením sloučit do společných skupin. Pro malé četnosti je také možné použít Fisher’s exact test – funkce fisher.test().

Porovnávání četností – test dobré shody \({\chi ^2}\) (goodness of fit)

Jsou odlišnosti od H0 jen dílem náhody?

H0: ano

H1: není to náhoda

Příklad: Na otázku „váš nejoblíbenější předmět“ odpovědělo \(n = 70\) žáků následovně:

  1. češtinu: 20
  2. matematiku: 21
  3. tělocvik: 29

Na hladině významnosti \(\alpha = 0.05\) rozhodněte mezi H0: předměty jsou stejně oblíbené, H1: mezi předměty jsou výrazné rozdíly.

Vzhledem k H0, která předpokládá rovnoměrné zastoupení, jsou očekávané hodnoty počtů všech skupin stejné. Můžeme tedy použít defaultní nastavení funkce chisq.test.

chisq.test(c(20, 21, 29))
## 
##  Chi-squared test for given probabilities
## 
## data:  c(20, 21, 29)
## X-squared = 2.0857, df = 2, p-value = 0.3524

P-hodnota není menší než 0.05, takže nezamítáme H0 na hladině významnosti \(\alpha = 0.05\). Pokud máme pocit, že tělocvik má navrch, musíme sehnat více dat a uvidíme, jak test dopadne. chisq.test(c(77, 37), p = c(1147/12, 1145/12), rescale.p = TRUE)

Specifikace očekávaných počtů

Očekávané počty můžeme specifikovat buď formou očekávaných pravděpodobností (parametr p), a nebo stejným parametrem zadat vektor očekávaných četností pro H0 a přidat parametr rescale.p = TRUE, který řekne funkci, aby si uvnitř tyto hodnoty přepočetla na pravděpodobnosti s celkovým součtem 1.

Příklad: V období duben až říjen proběhlo v městě 78 svateb, zatímco od listopadu do března jen 36. H0: svatby jsou napříč těmito obdobími rovnoměrně rozděleny, H1: nejsou rovnoměrně rozděleny.

Rok má 12 měsíců a pro očekávané četnosti je tedy dobré zohlednit nestejnoměrnou velikost sledovaných období. Duben až říjen je 7 měsíců, listopad až březen je 5 měsíců.

svatby <- c(78, 36)
celkem <- sum(svatby)
chisq.test(svatby, p = c(celkem*7/12, celkem*5/12), rescale.p = TRUE)
## 
##  Chi-squared test for given probabilities
## 
## data:  svatby
## X-squared = 4.7729, df = 1, p-value = 0.02891

P-hodnota je nižší než 0.05, proto přijímáme hypotézu H1 na hladině významnosti \(\alpha = 0.05\).

Získání souhrnné tabulky z data.frame jednotlivých dotázaných

Výše jsme pracovaly s vektory souhrnných hodnot. Pokud máme ale data.frame (získaný např. načtením tabulky ze souboru) s odpověďmi jednotlivých dotázaných, můžeme souhrn počtů podle požadovaných kategorií vyrobit funkcí table().

data(survey, package="MASS")  # interní příklad v R, dotazník 237 studentů statistiky na University of Adelaide

head(survey)
##      Sex Wr.Hnd NW.Hnd W.Hnd    Fold Pulse    Clap Exer Smoke Height
## 1 Female   18.5   18.0 Right  R on L    92    Left Some Never 173.00
## 2   Male   19.5   20.5  Left  R on L   104    Left None Regul 177.80
## 3   Male   18.0   13.3 Right  L on R    87 Neither None Occas     NA
## 4   Male   18.8   18.9 Right  R on L    NA Neither None Never 160.00
## 5   Male   20.0   20.0 Right Neither    35   Right Some Never 165.00
## 6 Female   18.0   17.7 Right  L on R    64   Right Some Never 172.72
##        M.I    Age
## 1   Metric 18.250
## 2 Imperial 17.583
## 3     <NA> 16.917
## 4   Metric 20.333
## 5   Metric 23.667
## 6 Imperial 21.000

Zaměřme se na otázku kouření (sloupec Smoke).

kuraci <- table(survey$Smoke)
kuraci
## 
## Heavy Never Occas Regul 
##    11   189    19    17

Řešme otázku, že H0: mezi studenty statistiky je 85 % nekuřáků a zbylé kategorie jsou rovnoměrně zastoupeny, každá 5 %. H1: rozdělení je jiné.

chisq.test(kuraci, p=c(.05, .85, .05, .05)) # pravděpodobnosti musejí odpovídat pořadí hodnot v tabulce
## 
##  Chi-squared test for given probabilities
## 
## data:  kuraci
## X-squared = 7.4098, df = 3, p-value = 0.05992

Na hladině významnosti \(\alpha = 0.05\) nezamítáme H0.

Očekávané četnosti můžeme zjistit přidáním expected a pomocí residuals vypíšeme hodnoty \(\chi\) bez umocnění, takže můžeme sledovat, kde jsou největší odchylky od očekávaných hodnot podle H0.

vysledek <- chisq.test(kuraci, p=c(.05, .85, .05, .05))
vysledek$expected
## Heavy Never Occas Regul 
##  11.8 200.6  11.8  11.8
vysledek$residuals
## 
##      Heavy      Never      Occas      Regul 
## -0.2328890 -0.8190163  2.0960010  1.5137785

Test nezávislosti (Chi-Squared test of independence)

Kontingenční tabulka – tabulka „se dvěma vstupy“. Je souvislost mezi oběma jevy?

Příklad: Dotazníku se zúčastnilo \(n = 400\) studentů. Mimo jiné odpovídali na otázky 1. „Byl(a) jste loni ubytován(a) na kolejích?“ (byl(a)/nebyl(a)) a 2. „Jaký byl váš loňský průměrný prospěch?“ (A: < 1.6, B: 1.6 – 2.1, C: > 2.1).

H0: ubytování na koleji s prospěchem nesouvisí, H1: souvisí.

Na \({\chi ^2}\) testu je dobré, že není nutné, aby byly kategorie rovnoměrně zastoupeny (např. kolej ano/ne). Stačí, když je v každé kategorii dostatečné zastoupení.

Jednou z možností, jak počty do R zadat, je vytvořit si dílčí vektory po řádcích, a ty pak spojit to jednoho data.frame. Pojmenování sloupců je sice hezké, ale pro test nemá význam.

kolejAno = c(39, 107, 93)
kolejNe = c(41, 73, 47)

tabulka <- as.data.frame(rbind(kolejAno, kolejNe))  # spojování vektorů jako řádky, z jmen vektorů se stanou jména řádků
tabulka
##          V1  V2 V3
## kolejAno 39 107 93
## kolejNe  41  73 47
names(tabulka) <- c("A", "B", "C")  # pojmenování sloupců - není nutné, ale je to hezké :-)
tabulka
##           A   B  C
## kolejAno 39 107 93
## kolejNe  41  73 47
chisq.test(tabulka)
## 
##  Pearson's Chi-squared test
## 
## data:  tabulka
## X-squared = 6.6286, df = 2, p-value = 0.03636

Na hladině významnosti \(\alpha = 0.05\) přijímáme H1.

Jak vytvořit kontingenční tabulku počtů (contingency table nebo cross tab) z data.frame s jednotlivými položkami?

Využijme interní data.frame mtcars a vytáhneme souhrnné počty testovaných automobilů podle počtu válců (cyl) a převodovky (am: 0 = automatická / 1 = manuální), pro více info o datech viz ?mtcars.

head(mtcars)
##                    mpg cyl disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1
tab <- table(mtcars$cyl, mtcars$am)
tab
##    
##      0  1
##   4  3  8
##   6  4  3
##   8 12  2

Nová kontingenční tabulka podle vybraných faktorů z již existující větší kontingenční tabulky (s více faktory)

data(UCBAdmissions)     # data o přijatých a nepřijatých uchazečích podle pohlaví a oddělení
DF <- as.data.frame(UCBAdmissions)
head(DF)
##      Admit Gender Dept Freq
## 1 Admitted   Male    A  512
## 2 Rejected   Male    A  313
## 3 Admitted Female    A   89
## 4 Rejected Female    A   19
## 5 Admitted   Male    B  353
## 6 Rejected   Male    B  207
tab <- xtabs(Freq ~ Gender + Admit, data = DF) # kontingenčí tabulka hodnot Freq dle Gender a Admit (Department nás tedy v tuto chvíli nezajímá)
tab
##         Admit
## Gender   Admitted Rejected
##   Male       1198     1493
##   Female      557     1278

Sloučení málo zastoupených kategorií

(Převzato z http://www.r-tutor.com/elementary-statistics/goodness-fit/chi-squared-test-independence)

Vraťme se k příkladu se studenty statistiky na University of Adelaide. Zajímá nás, zda je vztah mezi kuřáctvím a frekvencí sportování.

data(survey, package="MASS")  # interní příklad v R, dotazník 237 studentů statistiky na University of
tbl <- table(survey$Smoke, survey$Exer)  # vytvoření kontingenční tabulky závislosti 2 věcí
tbl
##        
##         Freq None Some
##   Heavy    7    1    3
##   Never   87   18   84
##   Occas   12    3    4
##   Regul    9    1    7

H0: kouření nemá souvislost s frekvencí sportování, H1: souvislost je.

chisq.test(tbl)
## Warning in chisq.test(tbl): Chi-squared approximation may be incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  tbl
## X-squared = 5.4885, df = 6, p-value = 0.4828

P-hodnota sice ukazuje na nezamítnutí H0, ale zároveň dostáváme varování o malém zastoupení kategorií.

Zkusme tedy některé sloučit kategorie None a Some do jedné společné.

ctbl <- cbind(tbl[,"Freq"], tbl[,"None"] + tbl[,"Some"])
ctbl
##       [,1] [,2]
## Heavy    7    4
## Never   87  102
## Occas   12    7
## Regul    9    8
chisq.test(ctbl)
## 
##  Pearson's Chi-squared test
## 
## data:  ctbl
## X-squared = 3.2328, df = 3, p-value = 0.3571

P-hodnota je sice nižší, ale stále nezamítáme H0 na hladině významnosti \(\alpha = 0.5\). A co je hlavní, nyní máme mnohem lepší pocit, že se nedopouštíme statistického zločinu, protože funkce již žádné varování nehlásí.

McNemarův \({\chi ^2}\) test – pro spárovaná měření

Situace: máme „metodu“, zkoumáme úspěch před a po.

H0: metoda nemá vliv, neboli úspěch->neúspěch == neúspěch->úspěch, H1: metoda má vliv (ať kladný či záporný).

Poznámka: McNemarův test je pouze aproximací pro tuto úlohu a zvláště pro nízké počty může být nepřesný. Zde ho tedy uvádíme spíše z historických důvodů. Vhodnější je pro tyto úlohy použít tzv. oboustranný znaménkový test popsaný v kapitole 16. Testy středních hodnot.

mat <- as.table(rbind(c(38, 14), 
                     c( 22, 36) ))
colnames(mat) <- rownames(mat) <- c("úspěch", "neúspěch")
names(dimnames(mat)) <- c("Před", "Po")
mat
##           Po
## Před       úspěch neúspěch
##   úspěch       38       14
##   neúspěch     22       36
mcnemar.test(mat, correct = TRUE)
## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  mat
## McNemar's chi-squared = 1.3611, df = 1, p-value = 0.2433

Na hladině významnosti \(\alpha = 0.05\) nezamítáme H0.

Poznámka: pro testy s více situacemi (více než 2) se používá Cochran’s Q Test.

Úlohy

  1. Viz domácí úkol.

© 16. 4. 2015 Tomáš Bořil, borilt@gmail.com