NMFM301
Praktické cvičení 4 – Instruktáž
Příprava k práci
Přesvědčte se, že pracujete se síťovou verzí R (na disku L:) a váš pracovní direktorář je umístěn na disku H: (nikoli C:).
Načtěte data pro dnešní cvičení a potřebné funkce:
load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/cv4.RData"))
load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/functions.RData"))Nyní máte k dispozici datové soubory sex a zam. Data sex pocházejí ze zdravotního průzkumu provedeného v USA v roce 2012. Obsahují údaje o respondentech starších 20 let. Data zam pocházejí ze sociálního průzkumu provedeného v ČR v roce 2004.
Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace”. Vypište si strukturu dat
str(sex)## 'data.frame': 3786 obs. of 9 variables:
## $ id : num 62161 62169 62172 62176 62179 ...
## $ active : Factor w/ 2 levels "Yes","No": 2 2 1 1 1 1 1 1 1 1 ...
## $ vek : num 22 21 43 34 55 35 26 30 57 42 ...
## $ pohl : Factor w/ 2 levels "Male","Female": 1 1 2 2 1 1 1 2 1 1 ...
## $ vzdelani : Factor w/ 5 levels "Elementary unfinished",..: 3 3 3 5 5 5 3 4 5 5 ...
## $ stav : Factor w/ 6 levels "Married","Widowed",..: 5 5 5 1 1 1 5 1 6 1 ...
## $ veksex : num NA NA 17 16 28 17 14 17 19 18 ...
## $ pocpart : num 0 0 4 15 3 4 10 1 36 2 ...
## $ pocpart.rok: num 0 0 2 1 1 1 1 1 1 1 ...str(zam)## 'data.frame': 2799 obs. of 8 variables:
## $ idno : num 1 2 3 4 5 6 7 8 9 10 ...
## $ tvtot : Factor w/ 8 levels "No time at all",..: 3 6 6 8 7 6 7 7 8 8 ...
## $ gndr : Factor w/ 2 levels "Male","Female": 2 2 1 2 2 2 1 2 2 2 ...
## $ agea : int 34 39 16 33 40 15 55 64 65 61 ...
## $ domicil : Factor w/ 5 levels "A big city","Suburbs or outskirts of big city",..: 3 4 4 4 4 4 4 4 4 4 ...
## $ eduyrs : int 17 11 10 13 18 9 12 16 16 12 ...
## $ unempl : num 0 0 0 0 0 0 0 0 0 0 ...
## $ regioncz: Factor w/ 14 levels "Prague","Central Bohemia",..: 14 14 12 12 12 12 12 12 12 12 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:227] 17 24 58 60 67 91 94 122 125 159 ...
## .. ..- attr(*, "names")= chr [1:227] "10316" "10323" "10357" "10359" ...a základní popisné statistiky
summary(sex)## id active vek pohl
## Min. :62161 Yes:3598 Min. :20.0 Male :1922
## 1st Qu.:64688 No : 188 1st Qu.:31.0 Female:1864
## Median :67094 Median :44.0
## Mean :67104 Mean :43.8
## 3rd Qu.:69574 3rd Qu.:56.0
## Max. :71915 Max. :69.0
##
## vzdelani stav veksex
## Elementary unfinished: 224 Married :1774 Min. : 9.0
## Elementary : 508 Widowed : 134 1st Qu.:15.0
## High school : 771 Divorced : 412 Median :17.0
## College unfinished :1236 Separated : 150 Mean :17.6
## College :1047 Single : 973 3rd Qu.:19.0
## Living with partner: 343 Max. :51.0
## NA's :193
## pocpart pocpart.rok
## Min. : 0.0 Min. : 0.00
## 1st Qu.: 2.0 1st Qu.: 0.00
## Median : 6.0 Median : 1.00
## Mean : 17.2 Mean : 1.13
## 3rd Qu.: 12.8 3rd Qu.: 1.00
## Max. :2000.0 Max. :69.00
## summary(zam)## idno tvtot gndr
## Min. : 1 More than 3 hours :761 Male :1295
## 1st Qu.: 802 More than 1,5 hours, up to 2 hours:459 Female:1504
## Median :3183 More than 2 hours, up to 2,5 hours:457
## Mean :2761 More than 2,5 hours, up to 3 hours:387
## 3rd Qu.:5038 More than 1 hour, up to 1,5 hours :297
## Max. :5833 0,5 hour to 1 hour :282
## (Other) :156
## agea domicil eduyrs
## Min. :14.0 A big city : 519 Min. : 6.0
## 1st Qu.:33.0 Suburbs or outskirts of big city: 138 1st Qu.:11.0
## Median :49.0 Town or small city :1360 Median :12.0
## Mean :48.2 Country village : 776 Mean :12.2
## 3rd Qu.:62.0 Farm or home in countryside : 6 3rd Qu.:13.0
## Max. :95.0 Max. :24.0
##
## unempl regioncz
## Min. :0.000 South Moravia : 335
## 1st Qu.:0.000 Moravian Silesia Reg.: 309
## Median :0.000 Prague : 262
## Mean :0.055 Central Bohemia : 229
## 3rd Qu.:0.000 Usti Reg. : 217
## Max. :1.000 Vysocina : 212
## (Other) :1235Dvouvýběrové testy
Nyní budeme pracovat s datovým souborem sex. Tento soubor obsahuje informace o sexuálním životě 4467 Američanů starších 20 let. Obsahuje veličiny active (byl/a někdy sexuálně aktivní?), veksex (věk sexuální iniciace - pouze u aktivních), pocpart (počet partnerů za celý život), pocpart.rok (počet partnerů za uplynulý rok) a dále demografická data (věk, pohlaví, vzdělání, rodinný status).
Porovnáme nejprve věk při sexuální iniciaci (prvním pohlavním styku) mužů a žen.
Začneme deskriptivní statistikou. Nakreslíme si krabicový diagram věku iniciace pro obě pohlaví
boxplot(veksex~pohl,data=sex)
a spočítáme průměry, mediány a rozptyly pro obě pohlaví. Použijeme k tomu funkci tapply (argument na.rm=TRUE říká, že se mají vynechat chybějící hodnoty veličiny veksex).
tapply(sex$veksex,sex$pohl,mean,na.rm=TRUE)## Male Female
## 17.35 17.90tapply(sex$veksex,sex$pohl,median,na.rm=TRUE)## Male Female
## 17 17tapply(sex$veksex,sex$pohl,var,na.rm=TRUE)## Male Female
## 18.24 15.54Ještě můžeme namalovat průměry věku iniciace a intervaly spolehlivosti pro střední věk iniciace pro obě pohlaví.
plotmeans(veksex~pohl,data=sex,xlab="Pohlavi",ylab="Prumerny vek iniciace")
Nakonec porovnáme empirické distribuční funkce věku iniciace obou pohlaví.
plot(ecdf(sex$veksex[sex$pohl=="Male"]),col="blue",cex.points=0.5,verticals=TRUE,
main="Sexualni iniciace")
lines(ecdf(sex$veksex[sex$pohl=="Female"]),col="red",cex.points=0.5,verticals=TRUE)
legend(40,0.4,lty=1,col=c("blue","red"),legend=c("Muzi","Zeny"))
Nasvědčují podle vás popisné statistiky a grafy tomu, že se věk iniciace mužů a žen může lišit?
Nyní přikročíme k testování několika různými metodami.
Dvouvýběrový Kolmogorovův-Smirnovův test
Dvouvýběrový Kolmogorovův-Smirnovův test testuje rovnost distribučních funkcí věku iniciace mužů a žen. Jeho testová statistika je \[
K_{n,m}=\sup_{x\in\mathbb{R}}\left|\widehat{F}_{X}(x)-\widehat{F}_{Y}(x)\right|.
\] Výsledky nám spočítá funkce ks.test.
ks.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])## Warning: p-value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## D = 0.0845, p-value = 5.431e-06
## alternative hypothesis: two-sidedOpět dostáváme varování o výskytu shodných hodnot. Věk iniciace zaokrouhlený na kalendářní rok fakticky není spojitá veličina, proto bychom Kolmogorovův-Smirnovův test raději neměli používat. Nicméně, p-hodnota je velice malá, kdybychom tomuto testu věřili navzdory porušeným předpokladům, zamítli bychom hypotézu o rovnosti distribučních funkcí s velkou rezervou.
Dvouvýběrový t-test
Proveďme nyní dvouvýběrový t-test na rovnost středních hodnot. Testová statistika pro hypotézu \(H_0: \mu_X=\mu_Y\) je \[ T_{n,m}=\sqrt{\frac{mn}{n+m}}\,\frac{\overline{X}_{n}-\overline{Y}_{m}}{S_{n,m}}, \] kde \[ S^2_{n,m}=\frac{n-1}{n+m-2}S^2_{X}+\frac{m-1}{n+m-2}S^2_{Y} \] je nestranný odhad společného rozptylu spočítaný z obou výběrů. Tento test předpokládá normalitu a shodné rozptyly. Ani jeden předpoklad není nejspíš splněn, i když porušení normality by při takhle velkém počtu pozorování vadit nemělo.
t.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"],var.equal=TRUE)##
## Two Sample t-test
##
## data: sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## t = -3.996, df = 3591, p-value = 6.583e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.8171 -0.2792
## sample estimates:
## mean of x mean of y
## 17.35 17.90I t-test zamítá rovnost středních hodnot u mužů a u žen s velkou rezervou.
Všimněte si, že výstup funkce t.test zahrnuje 95%-ní interval spolehlivosti pro rozdíl středního věku iniciace mužů a žen.
Dvouvýběrový z-test
Lépe by bylo použít dvouvýběrový z-test, který nepožaduje normalitu ani shodnost rozptylů. Jeho testová statistika je \[
Z_{n,m}=\frac{\overline{X}_{n}-\overline{Y}_{m}}{\sqrt{S^2_X/n+S^2_Y/m}}.
\] Tento test dostaneme funkcí t.test, v níž vynecháme argument var.equal=TRUE (požadující shodnost rozptylů). Kritické hodnoty a p-hodnoty se počítají Welchovou aproximací.
t.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"])##
## Welch Two Sample t-test
##
## data: sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## t = -3.998, df = 3578, p-value = 6.503e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.8169 -0.2794
## sample estimates:
## mean of x mean of y
## 17.35 17.90Welchův test dává prakticky stejné výsledky jako dvouvýběrový t-test se shodnými rozptyly. I zde máme k dispozici 95%-ní interval spolehlivosti pro rozdíl středních hodnot.
Dvouvýběrový Wilcoxonův test
Dvouvýběrový Wilcoxonův test testuje shodnost rozdělení (za předpokladu, že platí model posunutí v poloze), obecně jde o test hypotézy \(H_0: P[X_i\lt Y_j]=1/2\). Tento test počítá funkce wilcox.test. Testová statistika je \[
W^*_{n,m}=\sum_{i=1}^n\sum_{j=1}^m \mathbb{I}_{\{X_i\lt Y_j\}},
\] tj. odpovídá Mann-Whitneyho formulaci Wilcoxonova testu. Kritické hodnoty a p-hodnoty se počítají normální aproximací.
wilcox.test(sex$veksex[sex$pohl=="Male"],sex$veksex[sex$pohl=="Female"],correct=FALSE)##
## Wilcoxon rank sum test
##
## data: sex$veksex[sex$pohl == "Male"] and sex$veksex[sex$pohl == "Female"]
## W = 1453028, p-value = 2.06e-07
## alternative hypothesis: true location shift is not equal to 0Argument correct=FALSE zabraňuje provedení tzv. opravy na spojitost, kterou jsme do testové statistiky nezařazovali, ale R ji normálně provádí. Také Wilcoxonův test zamítá nulovou hypotézu. Výstup Wilcoxonova testu však neudává žádný odhad (bodový ani intervalový) rozdílu mezi porovnávanými populacemi.
Odhadování a testování pravděpodobností a šancí
Nyní budeme pracovat s datovým souborem zam. Tento soubor obsahuje data ze sociologického průzkumu provedeného v ČR v roce 2004. Zahrnuje veličiny idno (identifikační kód), tvtot (jak dlouho se obvykle denně díváte na televizi?), gndr (pohlaví), agea (věk), domicil (velikost sídla bydliště), eduyrs (počet let vzdělání), unempl (je nyní nezaměstnaný/á), regioncz (kraj bydliště).
Budeme se zabývat nezaměstnaností. Vypišme si, kolik lidí je nezaměstnaných:
table(zam$unempl)##
## 0 1
## 2645 154Je jich 154 z 2799 respondentů. Odhad pravděpodobnosti, že respondent je nezaměstnaný je \(\widehat{p}_n=154/2799=0.055\). Odhad šance na nezaměstnanost je \(\widehat{p}_n/(1-\widehat{p}_n)=0.0582\). Protože pravděpodobnost nezaměstnanosti je malá, šance je blízko hodnoty pravděpodobnosti.
Budeme testovat, zdali je míra (pravděpodobnost) nezaměstnanosti rovna 0.05. Použijeme funkci propor (není součástí R), která počítá přibližný interval spolehlivosti pro pravděpodobnost a testuje hypotézu \(H_0: p_X=p_0\) dvěma způsoby: jednak přímo pomocí věty 7.1, jednak přes šanci pomocí věty 7.2. Vstupní argumenty jsou tři: počet úspěchů \(X_n\), počet pokusů (velikost populace) \(n\), a hypotetická pravděpodobnost \(p_0\).
x <- sum(zam$unempl)
n <- nrow(zam)
propor(x,n,p0=0.05)##
## Testovani pravdepodobnosti
##
## 154 uspechu z 2799 pokusu
##
## Odhadnuta pravdepodobnost: 0.055
## Odhadnuta sance : 0.0582
## Odhadnuta log-sance : -2.84
##
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( 0.0466 , 0.0635 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika 1.165, p-hodnota 0.244
##
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.0472 , 0.0641 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika 1.218, p-hodnota 0.223Hypotézu, že pravděpodobnost nezaměstnanosti je 0.05, nemůžeme zamítnout. Obě metody dávají velmi podobné výsledky jak co se týče výsledku testu tak co se týče intervalu spolehlivosti pro parametr \(p_X\). R nám dává i možnost otestovat hypotézu přesným testem založeným na binomickém rozdělení a získat přesný interval spolehlivosti. Toto provádí funkce binom.test. Její vstupní argumenty jsou stejné, jako u funkce propor (jen poslední argument se jmenuje p místo p0).
binom.test(x,n,p=0.05)##
## Exact binomial test
##
## data: x and n
## number of successes = 154, number of trials = 2799, p-value =
## 0.2245
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
## 0.04686 0.06412
## sample estimates:
## probability of success
## 0.05502Její výsledky jsou velmi podobné oběma metodám z funkce propor, neboť pracujeme s dosti velkým počtem pozorování.
Součástí základní distribuce R je funkce prop.test, která provádí asymptotický test metodou 2, ale uvádí čtverec testové statistiky a používá asymptotické rozdělení \(\chi^2_1\):
prop.test(x,n,p=0.05,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: x out of n, null probability 0.05
## X-squared = 1.485, df = 1, p-value = 0.223
## alternative hypothesis: true p is not equal to 0.05
## 95 percent confidence interval:
## 0.04717 0.06409
## sample estimates:
## p
## 0.05502Nyní budeme zkoumat nezaměstnanost mezi některými subpopulacemi (což povede ke zmenšování rozsahu výběru).
Nezaměstnanost ve Zlínském kraji: [Nesmíme zapomenout správně měnit rozsah celé zkoumané populace \(n\).]
x <- sum(zam$unempl[zam$regioncz=="Zlin Reg."])
n <- sum(zam$regioncz=="Zlin Reg.")
propor(x,n,p0=0.05)##
## Testovani pravdepodobnosti
##
## 4 uspechu z 113 pokusu
##
## Odhadnuta pravdepodobnost: 0.0354
## Odhadnuta sance : 0.0367
## Odhadnuta log-sance : -3.31
##
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( 0.00133 , 0.0695 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika -0.84, p-hodnota 0.401
##
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.0133 , 0.0905 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika -0.7083, p-hodnota 0.479binom.test(x,n,0.05)##
## Exact binomial test
##
## data: x and n
## number of successes = 4, number of trials = 113, p-value = 0.6645
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
## 0.009728 0.088156
## sample estimates:
## probability of success
## 0.0354Nezaměstnanost na vesnicích ve Zlínském kraji:
x <- sum(zam$unempl[zam$regioncz=="Zlin Reg." & zam$domicil=="Country village"])
n <- sum(zam$regioncz=="Zlin Reg." & zam$domicil=="Country village")
propor(x,n,p0=0.05)##
## Testovani pravdepodobnosti
##
## 4 uspechu z 75 pokusu
##
## Odhadnuta pravdepodobnost: 0.0533
## Odhadnuta sance : 0.0563
## Odhadnuta log-sance : -2.88
##
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( 0.00248 , 0.104 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika 0.1285, p-hodnota 0.898
##
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.0202 , 0.134 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika 0.1324, p-hodnota 0.895Tady už je rozsah výběru tak malý (pouze 4 ,,úspěchy’’), že bychom rozhodně neměli používat asymptotické metody. Ale, jak ukazuje přesný test (dole), metoda 2 stále dává použitelné výsledky.
binom.test(x,n,0.05)##
## Exact binomial test
##
## data: x and n
## number of successes = 4, number of trials = 75, p-value = 0.7899
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
## 0.01472 0.13096
## sample estimates:
## probability of success
## 0.05333Nezaměstnanost mezi ženami na vesnicích ve Zlínském kraji:
x <- sum(zam$unempl[zam$regioncz=="Zlin Reg." & zam$domicil=="Country village" & zam$gndr=="Female"])
n <- sum(zam$regioncz=="Zlin Reg." & zam$domicil=="Country village" & zam$gndr=="Female")
propor(x,n,p0=0.05)##
## Testovani pravdepodobnosti
##
## 1 uspechu z 32 pokusu
##
## Odhadnuta pravdepodobnost: 0.0312
## Odhadnuta sance : 0.0323
## Odhadnuta log-sance : -3.43
##
## 95%-ni as. interval spolehlivosti pro pravdapodobnost - metoda 1:
## ( -0.029 , 0.0915 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 1:
## Testova statistika -0.6096, p-hodnota 0.542
##
## 95%-ni as. interval spolehlivosti pro pravdepodobnost - metoda 2:
## ( 0.00438 , 0.191 )
##
## As. test hypotezy p=0.05 proti p!=0.05 - typ 2:
## Testova statistika -0.4818, p-hodnota 0.63binom.test(x,n,0.05)##
## Exact binomial test
##
## data: x and n
## number of successes = 1, number of trials = 32, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.05
## 95 percent confidence interval:
## 0.0007909 0.1621710
## sample estimates:
## probability of success
## 0.03125V tomto extrémním případě už metoda 1 zcela selhává - interval spolehlivosti zahrnuje i záporná čísla a jeho horní okraj je špatný ve srovnání s přesným intervalem.
Máte-li čas, otestujte, zdali pravděpodobnost, že občan má aspoň 12 let školní docházky je rovna jedné polovině. Proveďte to v celé populaci, mezi ženami, mezi ženami staršími 50 let, a mezi ženami staršími 50 let žijícími na vesnici.