NMFM301
Praktické cvičení 3 – 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/cv3.RData"))
load(url("http://www.karlin.mff.cuni.cz/~pesta/NMFM301/functions.RData"))Nyní máte k dispozici datové soubory miry, tlak, sex, deti.
Data si prohlédněte poklepáním na příslušný řádek v okénku “Workspace”. Vypište si strukturu dat
str(miry)## 'data.frame': 500 obs. of 6 variables:
## $ id : num 62765 69815 71901 70078 68737 ...
## $ pohl : Factor w/ 2 levels "Male","Female": 2 2 2 2 1 1 2 1 2 1 ...
## $ vzdelani: Factor w/ 5 levels "Elementary unfinished",..: 5 4 4 2 2 4 4 1 2 1 ...
## $ stav : Factor w/ 6 levels "Married","Widowed",..: 1 2 1 1 1 1 1 5 1 5 ...
## $ vaha : num 51.9 90.6 140.4 55.7 101.4 ...
## $ vyska : num 166 164 153 150 178 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:385] 13 32 48 75 83 93 103 105 119 145 ...
## .. ..- attr(*, "names")= chr [1:385] "30" "62" "87" "130" ...str(tlak)## 'data.frame': 2309 obs. of 7 variables:
## $ id : num 62177 62178 62179 62191 62199 ...
## $ vek : num 51 80 55 70 57 62 65 77 54 80 ...
## $ pohl : Factor w/ 2 levels "Male","Female": 1 1 1 1 1 2 2 2 1 1 ...
## $ sys.tl.1: num 152 124 126 166 112 116 132 136 108 126 ...
## $ dia.tl.1: num 68 72 78 68 70 60 72 56 68 64 ...
## $ sys.tl.2: num 144 118 124 156 112 106 130 130 108 126 ...
## $ dia.tl.2: num 70 68 74 58 66 58 70 56 72 68 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:257] 1 44 47 51 52 53 54 60 76 93 ...
## .. ..- attr(*, "names")= chr [1:257] "14" "154" "170" "198" ...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(deti)## 'data.frame': 2695 obs. of 5 variables:
## $ id : num 62164 62172 62176 62189 62206 ...
## $ vek : num 44 43 34 30 35 62 22 65 77 38 ...
## $ vzdelani: Factor w/ 5 levels "Elementary unfinished",..: 4 3 5 4 4 1 3 4 1 4 ...
## $ stav : Factor w/ 6 levels "Married","Widowed",..: 1 5 1 1 1 4 4 1 5 5 ...
## $ poc.deti: num 0 2 2 2 0 5 1 3 0 1 ...
## - attr(*, "na.action")=Class 'omit' Named int [1:3364] 1 3 4 6 8 9 10 11 12 13 ...
## .. ..- attr(*, "names")= chr [1:3364] "1" "9" "10" "14" ...a základní popisné statistiky
summary(miry)## id pohl vzdelani
## Min. :62176 Male :264 Elementary unfinished: 55
## 1st Qu.:64654 Female:236 Elementary : 70
## Median :66958 High school :113
## Mean :67024 College unfinished :136
## 3rd Qu.:69495 College :126
## Max. :71909
## stav vaha vyska
## Married :249 Min. : 43.8 Min. :144
## Widowed : 46 1st Qu.: 67.7 1st Qu.:160
## Divorced : 47 Median : 79.5 Median :168
## Separated : 16 Mean : 82.5 Mean :168
## Single :102 3rd Qu.: 92.6 3rd Qu.:175
## Living with partner: 40 Max. :216.1 Max. :203summary(tlak)## id vek pohl sys.tl.1
## Min. :62177 Min. :50.0 Male :1152 Min. : 84
## 1st Qu.:64646 1st Qu.:56.0 Female:1157 1st Qu.:118
## Median :67275 Median :63.0 Median :130
## Mean :67118 Mean :64.6 Mean :132
## 3rd Qu.:69528 3rd Qu.:72.0 3rd Qu.:144
## Max. :71915 Max. :80.0 Max. :238
## dia.tl.1 sys.tl.2 dia.tl.2
## Min. : 0.0 Min. : 80 Min. : 0.0
## 1st Qu.: 64.0 1st Qu.:118 1st Qu.: 62.0
## Median : 72.0 Median :128 Median : 70.0
## Mean : 71.1 Mean :130 Mean : 69.5
## 3rd Qu.: 80.0 3rd Qu.:140 3rd Qu.: 78.0
## Max. :120.0 Max. :232 Max. :114.0summary(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(deti)## id vek vzdelani
## Min. :62164 Min. :20.0 Elementary unfinished:247
## 1st Qu.:64666 1st Qu.:34.0 Elementary :360
## Median :67117 Median :49.0 High school :530
## Mean :67099 Mean :48.9 College unfinished :866
## 3rd Qu.:69602 3rd Qu.:63.0 College :692
## Max. :71908 Max. :80.0
## stav poc.deti
## Married :1203 Min. : 0.00
## Widowed : 320 1st Qu.: 0.00
## Divorced : 323 Median : 2.00
## Separated : 120 Mean : 1.75
## Single : 534 3rd Qu.: 3.00
## Living with partner: 195 Max. :12.00Jednovýběrový Kolmogorovův-Smirnovův test
Kolmogorovův-Smirnovův test testuje hypotézu \(H_0: F_X(x)=F_0(x) \quad\forall x\in\mathbb{R}\) proti alternativě \(H_1: \exists x\in\mathbb{R}: F_X(x)\neq F_0(x)\), kde \(F_0\) je nějaká pevně specifikovaná spojitá distribuční funkce (bez neznámých parametrů).
Jeho testová statistika je \[ K_n=\sup_{x\in\mathbb{R}}\left|\widehat{F}_{n}(x)-F_0(x)\right|, \] kde \(\widehat{F}_{n}\) je empirická distribuční funkce náhodného výběru \(X_1,\ldots,X_n\).
Hypotézu zamítáme pro velké hodnoty \(K_n\) nebo \(\sqrt{n}K_n\). Používáme asymptotické kritické hodnoty vypočtené numericky.
Otestujeme hypotézu, že výška obyvatel USA zaznamenaná v datovém souboru miry se řídí rozdělením N(168,100). Nejprve spočítáme a namalujeme empirickou distribuční funkci a hypotetickou distribuční funkci. Empirická distribuční funkce se počítá funkcí ecdf.
plot(ecdf(miry$vyska),cex.points=0.5,
main="Empiricka distribucni funkce vysky")
xpts=seq(135,200,length=500)
lines(xpts,pnorm(xpts,mean=168,sd=10))
Funkce lines domaluje na předchozí obrázek propojené hodnoty hypotetické distribuční funkce spočítané v bodech xpts.
Provedeme Kolmogorovův-Smirnovův test pomocí funkce ks.test. Jako první argument má testovaný náhodný výběr, další tři argumenty specifikují hypotetickou distribuční funkci, poslední argument říká, že máme počítat asymptotický test (nikoli přesný).
ks.test(miry$vyska,"pnorm",mean=168,sd=10,exact=FALSE)## Warning: ties should not be present for the Kolmogorov-Smirnov test##
## One-sample Kolmogorov-Smirnov test
##
## data: miry$vyska
## D = 0.037, p-value = 0.4993
## alternative hypothesis: two-sidedR nás varuje, že v datech jsou shodné hodnoty (ties), čímž porušujeme předpoklady. Výsledek testu tedy nemusí být spolehlivý. Nyní budeme tento problém ignorovat.
Hodnota testové statistiky \(K_n\) byla 0.037, p-hodnota vyšla 0.4993. Na hladině \(\alpha=0.05\), nemůžeme zamítnout hypotézu, že výška má rozdělení N(168,100).
Nyní opakujte celý postup pro výšku mužů (miry$vyska[miry$pohl=="Male"]) a otestujte, zdali má výška mužů normální rozdělení se střední hodnotou 174 a směrodatnou odchylkou 7.5.
Znaménkový test
Znaménkový test testuje hypotézu \(H_0: m_X=m_0\) proti alternativě \(H_1: m_X\neq m_0\). Je založen na statistice \[ Y_n=\sum_{i=1}^n \mathbb{I}_{(0,\infty)}(X_i-m_0). \] Můžeme jej provést přesně (\(Y_n\) má za hypotézy binomické rozdělení) nebo asymptoticky pomocí statistiky \[ \frac{2}{\sqrt{n}}Y_n-\sqrt{n}, \] která má za \(H_0\) přibližně normované normální rozdělení.
Asymptotický test je implementován funkcí sign.test. Jejím prvním argumentem je náhodný výběr, druhým argumentem je hypotetický medián \(m_0\). Např. test hypotézy, že medián hmotnosti je 78 kg:
sign.test(miry$vaha,78)## n Y_n Test. statistika Krit. hodnota
## 500.0000 267.0000 1.5205 1.9600
## P-hodnota
## 0.1284Funkce sign.test počítá celkový počet pozorování n, hodnotu \(Y_n\), hodnotu asymptotické testové statistiky, kritickou hodnotu pro testování a p-hodnotu.
Otestujte nyní asymptotickým znaménkovým testem hypotézy, že medián hmotnosti mužů (žen) je 78 kg.
Přesný asymptotický test lze získat pomocí funkce binom.test (jde vlastně o test na parametr binomického rozdělení). Do této funkce musíme zadat hodnotu \(Y_n\) a hodnotu n.
binom.test(sum(miry$vaha>78),length(miry$vaha))##
## Exact binomial test
##
## data: sum(miry$vaha > 78) and length(miry$vaha)
## number of successes = 267, number of trials = 500, p-value =
## 0.1399
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.4892 0.5784
## sample estimates:
## probability of success
## 0.534Otestujte nyní přesným znaménkovým testem hypotézy, že medián hmotnosti mužů (žen) je 78 kg. Liší se nějak zásadně výsledky přesného testu od testu asymptotického? Proč asi?
Párové testy
Nyní budeme pracovat s datovým souborem tlak, který obsahuje 2309 pozorování. Nalezneme v něm dvě měření systolického tlaku v mm Hg (sys.tl.1 a sys.tl.2) vykonaná na téže osobě v určitém časovém intervalu (několik minut) a dvě měření diastolického tlaku (dia.tl.1 a dia.tl.2). Dále máme informace o věku (veličina vek, od 50 let výš) a pohlaví (veličina pohl) pacientů.
Nejprve spočítejme průměrný systolický tlak při prvním a druhém měření
mean(tlak$sys.tl.1)## [1] 132.1mean(tlak$sys.tl.2)## [1] 130a namalujme krabicový diagram obou měření. Takto se malují krabicové diagramy dvou veličin (nebo dvou skupin) vedle sebe:
boxplot(c(tlak$sys.tl.1,tlak$sys.tl.2)~rep(1:2,nrow(tlak)),main="Systolicky tlak",
names=c("1. mereni","2. mereni"))
Prohlédněme si ještě histogram rozdílů mezi prvním a druhým měřením (s 24 intervaly o délce 2 mm Hg):
hist(tlak$sys.tl.1-tlak$sys.tl.2,breaks=24)
Řekli byste podle obrázků, že mezi prvním a druhým měřením tlaku je nějaký systematický rozdíl?
Otestujme si, zdali se systematicky liší střední hodnota prvního a druhého měření systolického tlaku. Problém je třeba řešit párovým testem, protože se jedná o náhodný výběr dvojic pozorování \((X_i,Y_i)\) učiněných na tom samém pacientovi \(i\in\{1,\ldots,n\}\). Párové testy jsou vlastně jednovýběrové testy aplikované na rozdíly \(Z_i=X_i-Y_i\). Ukážeme si na tomto příkladě několik různých párových testů.
Párový t-test
Párový t-test spočítáme příkazem
t.test(tlak$sys.tl.1,tlak$sys.tl.2,paired=TRUE)##
## Paired t-test
##
## data: tlak$sys.tl.1 and tlak$sys.tl.2
## t = 14.27, df = 2308, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 1.842 2.429
## sample estimates:
## mean of the differences
## 2.135T-statistika nabyla hodnoty více než 14, proto s velikou rezervou zamítáme hypotézu o rovnosti středních hodnot. Můžeme tedy tvrdit, že střední hodnota systolického tlaku mezi prvním a druhým měřením klesla.
Párový znaménkový test
Párový znaménkový test testuje rovnost středních hodnot, jestliže rozdělení rozdílů mezi prvním a druhým měřením tlaku je symetrické. Pomocí histogramu rozdílů (viz výše) můžeme neformálně posoudit, zdali je tento předpoklad splněn. Myslíte, že ano?
Párový znaménkový test provedeme jako jednovýběrový znaménkový test aplikovaný na rozdíly mezi prvním a druhým měřením tlaku.
sign.test(tlak$sys.tl.1-tlak$sys.tl.2,m0=0)## n Y_n Test. statistika Krit. hodnota
## 2.309e+03 1.272e+03 4.891e+00 1.960e+00
## P-hodnota
## 1.006e-06Vidíme, že i znaménkový test zamítá hypotézu o rovnosti středních hodnot. Počet lidí, jimž se mezi oběma měřeními tlak snížil, byl příliš velký. Poznáte z výstupu funkce sign.test, kolik jich bylo?
Párový Wilcoxonův test
I párový Wilcoxonův test testuje rovnost středních hodnot za předpokladu symetrie rozdělení rozdílů mezi prvním a druhým měřením tlaku. Počítá se funkcí wilcox.test, argumenty vypadají podobně jako u párového t-testu.
wilcox.test(tlak$sys.tl.1,tlak$sys.tl.2,paired=TRUE,correct=FALSE)##
## Wilcoxon signed rank test
##
## data: tlak$sys.tl.1 and tlak$sys.tl.2
## V = 1371104, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Wilcoxonův test také zamítá hypotézu o rovnosti středních hodnot.
Dvouvýběrové testy
Nyní budeme pracovat s datovým souborem sex. Tento soubor obsahuje informace o sexuálním životě 4467 Američanů. 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.
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.