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.   :203
summary(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.0
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(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.00

Jednový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))

plot of chunk unnamed-chunk-4 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-sided

R 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.1284

Funkce 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.534

Otestujte 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.1
mean(tlak$sys.tl.2)
## [1] 130

a 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"))

plot of chunk unnamed-chunk-9

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)

plot of chunk unnamed-chunk-10 Ř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.135

T-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-06

Vidí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 0

Wilcoxonů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)

plot of chunk unnamed-chunk-14 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.90
tapply(sex$veksex,sex$pohl,median,na.rm=TRUE)
##   Male Female 
##     17     17
tapply(sex$veksex,sex$pohl,var,na.rm=TRUE)
##   Male Female 
##  18.24  15.54

Ješ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")

plot of chunk unnamed-chunk-16 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"))

plot of chunk unnamed-chunk-17 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-sided

Opě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.90

I 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.90

Welchův test dává prakticky stejné výsledky jako dvouvýběrový t-test se shodnými rozptyly.