5 Eloszlások vizsgálata
5.1 Grafikus vizsgálatok a normalitásra
A normalitás ellenőrzésére tipikusan a hisztogramot, a simított hisztogramot és a QQ ábrát használhatjuk.
5.1.1 Példa: a szisztolés vérnyomás normalitása
Vizsgáljuk meg a szisztolés vérnyomás normalitását a teljes mintában, valamint az egészséges és magas vérnyomásban szenvedők körében.
Forrás: (Taeger and Kuhnt 2014)
Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.
d <- read.table(file = textConnection("
no,status,mmhg
1,0,120
2,0,115
3,0,94
4,0,118
5,0,111
6,0,102
7,0,102
8,0,131
9,0,104
10,0,107
11,0,115
12,0,139
13,0,115
14,0,113
15,0,114
16,0,105
17,0,115
18,0,134
19,0,109
20,0,109
21,0,93
22,0,118
23,0,109
24,0,106
25,0,125
26,1,150
27,1,142
28,1,119
29,1,127
30,1,141
31,1,149
32,1,144
33,1,142
34,1,149
35,1,161
36,1,143
37,1,140
38,1,148
39,1,149
40,1,141
41,1,146
42,1,159
43,1,152
44,1,135
45,1,134
46,1,161
47,1,130
48,1,125
49,1,141
50,1,148
51,1,153
52,1,145
53,1,137
54,1,147
55,1,169
"), header=T, sep=",")
d$no <- factor(d$no)
d$status <- factor(d$status, labels=c("egészséges", "magas vérnyomás"))
str(d) # az adattábla szerkezete## 'data.frame': 55 obs. of 3 variables:
## $ no : Factor w/ 55 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ status: Factor w/ 2 levels "egészséges","magas vérnyomás": 1 1 1 1 1 1 1 1 1 1 ...
## $ mmhg : int 120 115 94 118 111 102 102 131 104 107 ...d # a beolvasott adattábla## no status mmhg
## 1 1 egészséges 120
## 2 2 egészséges 115
## 3 3 egészséges 94
## 4 4 egészséges 118
## 5 5 egészséges 111
## 6 6 egészséges 102
## 7 7 egészséges 102
## 8 8 egészséges 131
## 9 9 egészséges 104
## 10 10 egészséges 107
## 11 11 egészséges 115
## 12 12 egészséges 139
## 13 13 egészséges 115
## 14 14 egészséges 113
## 15 15 egészséges 114
## 16 16 egészséges 105
## 17 17 egészséges 115
## 18 18 egészséges 134
## 19 19 egészséges 109
## 20 20 egészséges 109
## 21 21 egészséges 93
## 22 22 egészséges 118
## 23 23 egészséges 109
## 24 24 egészséges 106
## 25 25 egészséges 125
## 26 26 magas vérnyomás 150
## 27 27 magas vérnyomás 142
## 28 28 magas vérnyomás 119
## 29 29 magas vérnyomás 127
## 30 30 magas vérnyomás 141
## 31 31 magas vérnyomás 149
## 32 32 magas vérnyomás 144
## 33 33 magas vérnyomás 142
## 34 34 magas vérnyomás 149
## 35 35 magas vérnyomás 161
## 36 36 magas vérnyomás 143
## 37 37 magas vérnyomás 140
## 38 38 magas vérnyomás 148
## 39 39 magas vérnyomás 149
## 40 40 magas vérnyomás 141
## 41 41 magas vérnyomás 146
## 42 42 magas vérnyomás 159
## 43 43 magas vérnyomás 152
## 44 44 magas vérnyomás 135
## 45 45 magas vérnyomás 134
## 46 46 magas vérnyomás 161
## 47 47 magas vérnyomás 130
## 48 48 magas vérnyomás 125
## 49 49 magas vérnyomás 141
## 50 50 magas vérnyomás 148
## 51 51 magas vérnyomás 153
## 52 52 magas vérnyomás 145
## 53 53 magas vérnyomás 137
## 54 54 magas vérnyomás 147
## 55 55 magas vérnyomás 169Első lépésként végezzük el a normalitás ellenőrzését hisztogrammal, és QQ-ábrával. A középső képen a hisztogram mellett a simított hisztogramot és a mintából számolt paraméterekhez tartozó normális eloszlás görbéjét is berajzoltuk (a piros görbe).
par(mfrow=c(1,3)) # a grafikus eszköz felosztása 1 sorra és 3 oszlopra
hist(d$mmhg) # normalitás ellenőrzése hisztogrammal
hist(d$mmhg, freq = F); lines(density(d$mmhg)) # normalitás ellenőrzése hisztogrammal, és simított gyakorásgi görbével
curve(dnorm(x,mean = mean(d$mmhg), sd = sd(d$mmhg)), add = T, col = "red")
qqnorm(d$mmhg);qqline(d$mmhg) # normalitás ellenőrzése QQ-ábrával 
par(mfrow=c(1,1)) # az alapértelmezett felosztás visszaállítása a grafikus eszközönAz egészségesek és a magas vérnyomásban szenvedők csoportjában is végezzük el a fenti grafikus vizsgálatokat.
library(ggplot2)
library(PASWR2)
range.x <- range(d$mmhg) # a közös x tengely miatt
plots <- lapply(X = split(d, d$status), FUN = function(x) {
ggplot(x, aes(x=mmhg)) + geom_histogram(aes(y = ..density..), binwidth = 10, right = F, fill="white", colour = "black") + geom_density() + facet_grid(status ~ .) + stat_function( fun = dnorm, args = list(mean = mean(x$mmhg), sd = sd(x$mmhg)), colour = "red") + theme_bw() + coord_cartesian(xlim = range.x)
})
multiplot(plotlist = plots, cols = 1) # ábrák rajzolása
library(ggplot2)
library(PASWR2)
p1 <- ggplot(d, aes(sample = mmhg)) + stat_qq() + facet_grid(status ~ .) + theme_bw() # QQ-ábrák
p2 <- ggplot(d, aes(x=mmhg)) + geom_density(aes(fill=status), size=1, alpha=.5) + theme_bw() # simított hisztogramok
multiplot(p1, p2, cols = 2) # ábrák rajzolása
A fenti ábrák alapján azt állapíthatjuk meg, hogy a teljes mintára vonatkozó normalitás sérülhet, a két csoportban azonban a normális eloszlást nagy biztonsággal feltételezhetjük.
5.2 Normalitás ellenőrzése ferdeségi és csúcsossági együtthatókkal
A normalitás ellenőrzése alapozható a két alakmutatóra is, a ferdeségi és a csúcsossági együtthatóra, amelyek értéke normális eloszlás esetén nullával egyenlő. Kiszámításukra több módszer is ismert, ezeket lentebb ismertetjük.
A ferdeségi együttható néhány értelmezése:
- \(\gamma_1 = \frac{m_3}{m_2^{3/2}}\), ahol \(m_1=\frac{\sum_{j=1}^{n}X_j}{n}\), \(m_i=\frac{\sum_{j=1}^{n}(X_j-m_1)^i}{n}\), \(i={2, 3, 4}\)
- \(G_1 = \gamma_1 \frac{\sqrt{n(n-1)}}{n-2}\)
- \(b_1 = \frac{m_3}{s^3} = \gamma_1 \left(\frac{n-1}{n}\right)^{3/2}\)
A csúcsossági együtthatóra is több számolási mód adható:
- \(\gamma_2 = \frac{m_4}{m_2^2} - 3\)
- \(G_2 = \frac{((n+1) \gamma_2 + 6)(n-1) }{ (n-2)(n-3)}\)
- \(b_2 = \frac{m_4}{s^4} - 3 = (\gamma_2 + 3) (1 - \frac{1}{n})^2 - 3\)
A standard hiba aszimptotikus kiszámítása a ferdeségi és csúcsossági együttható esetén:
\[ASE.skew = \sqrt{ \frac{6n(n-1)}{(n-2)(n+1)(n+3)} }\] \[ASE.kurt = \sqrt{ \frac{n^2 - 1}{(n-3)(n+5)} }\]
A standard hiba egyszerűsített számítása:
\[SE.skew = \sqrt{ \frac{6}{n} }\] \[SE.kurt = \sqrt{ \frac{24}{n} }\]
5.2.1 Példa: a szisztolés vérnyomás normalitása
Vizsgáljuk meg a szisztolés vérnyomás normalitását a teljes mintában, valamint az egészséges és magas vérnyomásban szenvedők körében.
Forrás: (Taeger and Kuhnt 2014)
Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.
d <- read.table(file = textConnection("
no,status,mmhg
1,0,120
2,0,115
3,0,94
4,0,118
5,0,111
6,0,102
7,0,102
8,0,131
9,0,104
10,0,107
11,0,115
12,0,139
13,0,115
14,0,113
15,0,114
16,0,105
17,0,115
18,0,134
19,0,109
20,0,109
21,0,93
22,0,118
23,0,109
24,0,106
25,0,125
26,1,150
27,1,142
28,1,119
29,1,127
30,1,141
31,1,149
32,1,144
33,1,142
34,1,149
35,1,161
36,1,143
37,1,140
38,1,148
39,1,149
40,1,141
41,1,146
42,1,159
43,1,152
44,1,135
45,1,134
46,1,161
47,1,130
48,1,125
49,1,141
50,1,148
51,1,153
52,1,145
53,1,137
54,1,147
55,1,169
"), header=T, sep=",")
d$no <- factor(d$no)
d$status <- factor(d$status, labels=c("egészséges", "magas vérnyomás"))
str(d) # az adattábla szerkezete## 'data.frame': 55 obs. of 3 variables:
## $ no : Factor w/ 55 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ status: Factor w/ 2 levels "egészséges","magas vérnyomás": 1 1 1 1 1 1 1 1 1 1 ...
## $ mmhg : int 120 115 94 118 111 102 102 131 104 107 ...d # a beolvasott adattábla## no status mmhg
## 1 1 egészséges 120
## 2 2 egészséges 115
## 3 3 egészséges 94
## 4 4 egészséges 118
## 5 5 egészséges 111
## 6 6 egészséges 102
## 7 7 egészséges 102
## 8 8 egészséges 131
## 9 9 egészséges 104
## 10 10 egészséges 107
## 11 11 egészséges 115
## 12 12 egészséges 139
## 13 13 egészséges 115
## 14 14 egészséges 113
## 15 15 egészséges 114
## 16 16 egészséges 105
## 17 17 egészséges 115
## 18 18 egészséges 134
## 19 19 egészséges 109
## 20 20 egészséges 109
## 21 21 egészséges 93
## 22 22 egészséges 118
## 23 23 egészséges 109
## 24 24 egészséges 106
## 25 25 egészséges 125
## 26 26 magas vérnyomás 150
## 27 27 magas vérnyomás 142
## 28 28 magas vérnyomás 119
## 29 29 magas vérnyomás 127
## 30 30 magas vérnyomás 141
## 31 31 magas vérnyomás 149
## 32 32 magas vérnyomás 144
## 33 33 magas vérnyomás 142
## 34 34 magas vérnyomás 149
## 35 35 magas vérnyomás 161
## 36 36 magas vérnyomás 143
## 37 37 magas vérnyomás 140
## 38 38 magas vérnyomás 148
## 39 39 magas vérnyomás 149
## 40 40 magas vérnyomás 141
## 41 41 magas vérnyomás 146
## 42 42 magas vérnyomás 159
## 43 43 magas vérnyomás 152
## 44 44 magas vérnyomás 135
## 45 45 magas vérnyomás 134
## 46 46 magas vérnyomás 161
## 47 47 magas vérnyomás 130
## 48 48 magas vérnyomás 125
## 49 49 magas vérnyomás 141
## 50 50 magas vérnyomás 148
## 51 51 magas vérnyomás 153
## 52 52 magas vérnyomás 145
## 53 53 magas vérnyomás 137
## 54 54 magas vérnyomás 147
## 55 55 magas vérnyomás 169Alakmutatók kiszámolása. A ferdeségi és csúcsossági együtthatókat a DescTools csomag Skew() és Kurt() függvényével is kiszámolhatjuk. A különböző mutatók értékét a method= argumentum megadásával számolhatjuk ki (method=1 esetén a \(\gamma_1\) és \(\gamma_2\), a method=2 a \(G_1\) és \(G_2\), a method=3 pedig a \(b_1\) és \(b_2\) mutatókat számolja). A conf.level= argumentum megadásával adott megbízhatósági szinten konfidenciaintervallumot készíthetünk a populációbeli paraméterekre. A konfidenciaintervallumok határai a két függvényben alapértelmezés szerint a bootstrap eljárással készülnek.
Első lépésben számoljuk ki a teljes mintában az alakmutatókat:
library(DescTools)
DescTools::Skew(d$mmhg, conf.level = 0.95, method = 1) # a ferdeség hagyományos, tankönyvi definíciója## skew lwr.ci upr.ci
## -0.08008799 -0.49612203 0.35835211DescTools::Skew(d$mmhg, conf.level = 0.95, method = 2) # ferdeség a SAS és SPSS programcsomagokban## skew lwr.ci upr.ci
## -0.08235117 -0.60347390 0.32393525DescTools::Skew(d$mmhg, conf.level = 0.95, method = 3) # ferdeség a MINITAB és BMDP programcsomagokban## skew lwr.ci upr.ci
## -0.07791373 -0.49382762 0.42055716DescTools::Kurt(d$mmhg, conf.level = 0.95, method = 1) # a csúcsosság hagyományos, tankönyvi definíciója## kurt lwr.ci upr.ci
## -1.0588047 -1.3778092 -0.5580666DescTools::Kurt(d$mmhg, conf.level = 0.95, method = 2) # csúcsosság a SAS és SPSS programcsomagokban## kurt lwr.ci upr.ci
## -1.044204 -1.399506 -0.452023DescTools::Kurt(d$mmhg, conf.level = 0.95, method = 3) # csúcsosság a MINITAB és BMDP programcsomagokban## kurt lwr.ci upr.ci
## -1.1287519 -1.4527139 -0.6185631Alakmutatók értékelése. A fenti outputból kiolvashatjuk a különböző ferdeségi és csúcsossági együtthatók értékét. Az értékelést a 3. típus (method = 3) alapján végzünk el. A ferdeségi együttható értéke \(b_1=-0,078\), a 95%-os megbízhatóságú konfidenciaintervallum határai tartalmazzák a nullát, így a ferdeségre vonatkozó ellenhipotézist nem tudjuk elfogadni. A csúcsossági együttható \(b_2=-1.129\), a 95%-os konfidenciaintervallum pedig nem tartalmazza a nullát, amely a csúcsosságra vonatkozó ellenhipotézis elfogadását jelenti. Ezek alapján azt mondhatjuk, hogy a változó eloszlása a csúcsossági paraméterre vonatkozó 95%-os konnfidenciaintervallum alapján eltér a normális eloszlástól.
Az egyes csoportokra is meghatározzuk a ferdeségi és csúcsossági együtthatókat. Most csak a \(b_1\) és \(b_2\) mutatókat számoljuk ki:
by(d$mmhg, d$status, DescTools::Skew, conf.level = 0.95, method = 3)## d$status: egészséges
## skew lwr.ci upr.ci
## 0.4675099 -0.1677403 1.2485703
## --------------------------------------------------------
## d$status: magas vérnyomás
## skew lwr.ci upr.ci
## -0.1011255 -0.7803134 0.7719545by(d$mmhg, d$status, DescTools::Kurt, conf.level = 0.95, method = 3)## d$status: egészséges
## kurt lwr.ci upr.ci
## -0.1037284 -1.0393472 1.7619865
## --------------------------------------------------------
## d$status: magas vérnyomás
## kurt lwr.ci upr.ci
## -0.02156674 -0.84692168 1.66162622Az egyes mutatók leolvashatók a fenti outputból, a 95%-os konfidenciaintervallumok mindegyike tartalmazza a nulla értéket, amelyek így a normalitást támogatják.
Alternatívák alakmutatók számítására. Az R-ben számos függvény létezik a fenti mutatók kiszámítására. A pastecs csomag stat.desc() függvénye más mutatók mellett a ferdeségi és a csúcsossági mutatót is szolgáltatja a norm=T argumentum megadása után. A stat.desc() függvény a \(b_1\) és \(b_2\) értékeket számolja ki. Az outputban a skew.2SE és kurt.2SE címkével jelzett értékek is szerepelnek, melyek számítási módja: \(skew.2SE = \frac{b_1}{2ASE.skew}\), \(kurt.2SE = \frac{b_2}{2ASE.kurt}\). Ha ezek a mutatók -1 és +1 között vannak, akkor támogatják a változó normalitását.
library(pastecs)
stat.desc(d$mmhg, norm = T, basic = F, desc = F)## skewness skew.2SE kurtosis kurt.2SE normtest.W normtest.p
## -0.07791373 -0.12108103 -1.12875193 -0.89087549 0.96077467 0.07011906A QuantPsyc csomag Skew(), Kurt() és norm() függvényei a \(G_1\) és \(G_2\) mutatókat szolgáltatják. Az SE oszlopban az SE.skew, illetve SE.kurt értékek olvashatók, a t oszlopban a két előző érték hányadosa, és a hozzá tartozó p érték a standard normális eloszlás alapján.
library(QuantPsyc)
QuantPsyc::Skew(d$mmhg)## sk se.sk t.sk p.sk
## [1,] -0.08235117 0.3302891 -0.2493305 0.4015526QuantPsyc::Kurt(d$mmhg)## ku se.ku t.ku p.ku
## [1,] -1.044204 0.6605783 -1.580742 0.05696853QuantPsyc::norm(d$mmhg)## Statistic SE t-val p
## Skewness -0.08235117 0.3302891 -0.2493305 0.40155256
## Kurtosis -1.04420373 0.6605783 -1.5807419 0.056968535.3 D’Agostino–Pearson-próba
A D’Agostino–Pearson-próba a ferdeségi és a csúcsossági együtthatók alapján a változó normalitására vonatkozó nullhipotézist teszteli.
Alkalmazási feltételek:
- folytonos eloszlású \(X\) változó
Null- és ellenhipotézisek:
- Nullhipotézis:
- \(H_0: X\) normális eloszlású
- Ellenhipotézis:
- \(H_1: X\) nem normális eloszlású
Próbastatisztika:
A ferdeségre vonatkozó Z3 próbastatisztika a nullhipotézis igaz volta esetén standard normális eloszlást követ. Kiszámítása:
\[Y= \gamma_1 \sqrt{\frac{(n + 1)(n + 3)}{(6 (n - 2))}}, \beta_2(\gamma_1) = \frac{3 (n^2 + 27 n - 70) (n + 1) (n + 3)}{(n - 2) (n + 5) (n + 7) (n + 9)}\] \[W= \sqrt{-1 + \sqrt{2 (\beta_2(\gamma_1) - 1)}}, \delta = \frac{1}{\sqrt{\ln(W)}}, \alpha = \sqrt{\frac{2}{W^2 - 1}}\] \[Z3= \delta \ln\left(\frac{Y}{\alpha} + \sqrt{\left(\frac{Y}{\alpha}\right)^2 + 1}\right) \sim N(0,1)\]
A csúcsosságra vonatkozó Z4 próbastatisztika szintén standard normális eloszlású, ha \(H_0\) igaz. Ennek is megadjuk a kiszámítását:
\[\gamma_2 = \gamma_2 + 3, \mu(\gamma_2) = \frac{3 (n - 1)}{n + 1}, \sigma(\gamma_2) = \frac{24 n (n - 2) (n - 3)}{(n + 1)^2 (n + 3) (n + 5)}, x = \frac{\gamma_2 - \mu(\gamma_2)}{\sqrt(\sigma(\gamma_2))}\] \[\beta_1(\gamma_2) = \frac{6 (n^2 - 5n + 2)}{(n + 7) (n + 9)} \sqrt{\frac{6 (n + 3) (n + 5)}{n (n - 2) (n - 3)}}, A = 6 + \frac{8}{\beta_1(\gamma_2)} \left(\frac{2}{\beta_1(\gamma_2)} + \sqrt{1 + \frac{4}{\beta_1(\gamma_2)^2}}\right)\] \[Z4 = \frac{(1 - \frac{2}{9 * A}) - \left(\frac{1 - \frac{2}{A}}{1 + x\sqrt{\frac{2}{A - 4}}}\right)^\frac{1}{3}}{\sqrt{\frac{2}{9A}}} \sim N(0,1)\]
A D’Agostino–Pearson-próba \(\chi^2\) próbastatisztikája a fenti \(Z3\) és \(Z4\) értékeken alapul:
\[\chi^2 = Z3^2 + Z4^2 \sim \chi^2(2)\]
5.3.1 Példa: a szisztolés vérnyomás normalitása
Vizsgáljuk meg a szisztolés vérnyomás normalitását a teljes mintában, és az egészséges és magas vérnyomásban szenvedők körében is.
Forrás: (Taeger and Kuhnt 2014)
Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.
d <- read.table(file = textConnection("
no,status,mmhg
1,0,120
2,0,115
3,0,94
4,0,118
5,0,111
6,0,102
7,0,102
8,0,131
9,0,104
10,0,107
11,0,115
12,0,139
13,0,115
14,0,113
15,0,114
16,0,105
17,0,115
18,0,134
19,0,109
20,0,109
21,0,93
22,0,118
23,0,109
24,0,106
25,0,125
26,1,150
27,1,142
28,1,119
29,1,127
30,1,141
31,1,149
32,1,144
33,1,142
34,1,149
35,1,161
36,1,143
37,1,140
38,1,148
39,1,149
40,1,141
41,1,146
42,1,159
43,1,152
44,1,135
45,1,134
46,1,161
47,1,130
48,1,125
49,1,141
50,1,148
51,1,153
52,1,145
53,1,137
54,1,147
55,1,169
"), header=T, sep=",")
d$no <- factor(d$no)
d$status <- factor(d$status, labels=c("egészséges", "magas vérnyomás"))
str(d) # az adattábla szerkezete## 'data.frame': 55 obs. of 3 variables:
## $ no : Factor w/ 55 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ status: Factor w/ 2 levels "egészséges","magas vérnyomás": 1 1 1 1 1 1 1 1 1 1 ...
## $ mmhg : int 120 115 94 118 111 102 102 131 104 107 ...d # a beolvasott adattábla## no status mmhg
## 1 1 egészséges 120
## 2 2 egészséges 115
## 3 3 egészséges 94
## 4 4 egészséges 118
## 5 5 egészséges 111
## 6 6 egészséges 102
## 7 7 egészséges 102
## 8 8 egészséges 131
## 9 9 egészséges 104
## 10 10 egészséges 107
## 11 11 egészséges 115
## 12 12 egészséges 139
## 13 13 egészséges 115
## 14 14 egészséges 113
## 15 15 egészséges 114
## 16 16 egészséges 105
## 17 17 egészséges 115
## 18 18 egészséges 134
## 19 19 egészséges 109
## 20 20 egészséges 109
## 21 21 egészséges 93
## 22 22 egészséges 118
## 23 23 egészséges 109
## 24 24 egészséges 106
## 25 25 egészséges 125
## 26 26 magas vérnyomás 150
## 27 27 magas vérnyomás 142
## 28 28 magas vérnyomás 119
## 29 29 magas vérnyomás 127
## 30 30 magas vérnyomás 141
## 31 31 magas vérnyomás 149
## 32 32 magas vérnyomás 144
## 33 33 magas vérnyomás 142
## 34 34 magas vérnyomás 149
## 35 35 magas vérnyomás 161
## 36 36 magas vérnyomás 143
## 37 37 magas vérnyomás 140
## 38 38 magas vérnyomás 148
## 39 39 magas vérnyomás 149
## 40 40 magas vérnyomás 141
## 41 41 magas vérnyomás 146
## 42 42 magas vérnyomás 159
## 43 43 magas vérnyomás 152
## 44 44 magas vérnyomás 135
## 45 45 magas vérnyomás 134
## 46 46 magas vérnyomás 161
## 47 47 magas vérnyomás 130
## 48 48 magas vérnyomás 125
## 49 49 magas vérnyomás 141
## 50 50 magas vérnyomás 148
## 51 51 magas vérnyomás 153
## 52 52 magas vérnyomás 145
## 53 53 magas vérnyomás 137
## 54 54 magas vérnyomás 147
## 55 55 magas vérnyomás 169A D’Agostino–Pearson-próba végrehajtása. A D’Agostino–Pearson-próba végrehajtását az fBasics csomag dagoTest() függvényével kezdeményezhetjük. Egyetlen paramétere a numerikus adatminta. A próbát elvégezzük a teljes mintára és a két csoportra is (egészségesek és magas vérnyomásúak csoportja).
library(fBasics)
dagoTest(x = d$mmhg) # D'Agostino--Pearson-próba a teljes mintára##
## Title:
## D'Agostino Normality Test
##
## Test Results:
## STATISTIC:
## Chi2 | Omnibus: 7.802
## Z3 | Skewness: -0.268
## Z4 | Kurtosis: -2.7803
## P VALUE:
## Omnibus Test: 0.02022
## Skewness Test: 0.7887
## Kurtosis Test: 0.00543
##
## Description:
## Wed May 18 23:11:26 2016 by user: kali2library(fBasics)
by(d$mmhg, d$status, dagoTest) # D'Agostino--Pearson-próba a két csoportra##
## Title:
## D'Agostino Normality Test
##
## Test Results:
## STATISTIC:
## Chi2 | Omnibus: 1.9236
## Z3 | Skewness: 1.1804
## Z4 | Kurtosis: 0.7282
## P VALUE:
## Omnibus Test: 0.3822
## Skewness Test: 0.2378
## Kurtosis Test: 0.4665
##
## Description:
## Wed May 18 23:11:26 2016 by user: kali2
##
## --------------------------------------------------------
##
## Title:
## D'Agostino Normality Test
##
## Test Results:
## STATISTIC:
## Chi2 | Omnibus: 0.6472
## Z3 | Skewness: -0.2782
## Z4 | Kurtosis: 0.7548
## P VALUE:
## Omnibus Test: 0.7236
## Skewness Test: 0.7809
## Kurtosis Test: 0.4504
##
## Description:
## Wed May 18 23:11:26 2016 by user: kali2A D’Agostino–Pearson-próba értékelése. A D’Agostino–Pearson-próba értékelését most a teljes mintára nézve végezzük el. A normalitás ellenőrzésére végrehajtott D’Agostino–Pearson-próba eredménye alapján, a normális eloszlásra vonatkozó nullhipotézist visszautasítjuk (\(\chi^2(2)=7,802\); \(p=0,0202\)). A normális eloszlástól való eltérés az alakmutatók alapján nem a szimmetrikus tulajdonságon (\(Z3=-0,268\); \(p=0,7887\)), hanem inkábba a csúcsosságon múlik (\(Z4=-2,7803\); \(p=0,0054\)).
A D’Agostino–Pearson-próba kézi számolása. A D’Agostino–Pearson-próba outputjában szereplő próbastatisztika értékeket magunk is kiszámolhatjuk:
n <- sum(!is.na(d$mmhg)) # mintaelemszám
m.1 <- sum(d$mmhg)/n # 1. momentum, mintaátlag
m.2 <- sum((d$mmhg-m.1)^2)/n # 2. momentum, variancia
m.3 <- sum((d$mmhg-m.1)^3)/n # 3. momentum
m.4 <- sum((d$mmhg-m.1)^4)/n # 4. momentum
gamma.1 <- m.3/m.2^(3/2) # ferdeségi mutató, gamma1
gamma.2 <- (m.4/m.2^2)-3 # csúcsossági mutató, gamma2
gamma.1 # ferdeségi mutató, gamma1 kiíratása## [1] -0.08008799gamma.2 # csúcsossági mutató, gamma2 kiíratása## [1] -1.058805# Z3 számolása
Y <- gamma.1 * sqrt((n + 1) * (n + 3) / (6 * (n - 2)))
b2ng1 <- (3 * (n^2 + 27 * n - 70) * (n + 1) * (n + 3)) / ((n - 2) * (n + 5) * (n + 7) * (n + 9))
W <- sqrt(-1 + sqrt(2 * (b2ng1 - 1)))
delta <- 1 / sqrt(log(W))
alfa <- sqrt(2 / (W^2 - 1))
Z3 <- delta * log(Y / alfa + sqrt((Y / alfa)^2 + 1)) # Z3 próbastatisztika a ferdeségi próbához
Z3 # Z3 próbastatisztika kiíratása## [1] -0.2680036# Z4 számolása
gamma.2 <- gamma.2 + 3
Eg2 <- 3 * (n - 1) / (n + 1)
varg2 <- (24 * n * (n - 2) * (n - 3)) / ((n + 1)^2 * (n + 3) * (n + 5))
x <- (gamma.2 - Eg2) / sqrt(varg2)
ngbg2 <- (6 * (n^2 - 5 * n + 2)) / ((n + 7) * (n + 9)) * sqrt((6 * (n + 3) * (n + 5)) / (n * (n - 2) * (n - 3)))
A <- 6 + 8/ngbg2 * (2/ngbg2 + sqrt(1 + 4/ngbg2^2))
# Z4 próbastatisztika a csúcsossági próbához
Z4 <- ((1 - 2/(9 * A)) - ((1 - 2/A)/(1 + x*sqrt(2/(A - 4))))^(1/3))/sqrt(2/(9 * A))
Z4 # Z4 próbastatisztika kiírása## [1] -2.780328# Chi2 számolása
Chi2 <- Z3^2 + Z4^2 # Chi2 próbastatisztika a D'Agostino--Pearson-próbához
Chi2 # Chi2 próbastatisztika kiírása## [1] 7.8020482 * pnorm (q = Z3, lower.tail = TRUE) # p érték a Z3 ferdeségi próbastatisztikához## [1] 0.78869652 * pnorm (q = Z4, lower.tail = TRUE) # p érték a Z4 csúcsossági próbastatisztikához## [1] 0.005430407pchisq(q = Chi2, df = 2, lower.tail = FALSE) # p érték a Chi2 normalitási próbastatisztikához## [1] 0.02022119Alternatívák a D’Agostino–Pearson-próbára. A D’Agostino–Pearson-próba outputjában a ferdeségi és csúcsossági próbastatisztika értékek (Z3 és Z4), és a hozzájuk tartozó p értékek is szerepelnek. A moments csomag agostino.test() és anscombe.test() függvényeivel a ferdeségre vonatkozó D’Agostino-próba és a csúcsosságra vonatkozó Anscombe–Glynn-próba külön-külön is végrehajtható:
library(moments)
agostino.test(x = d$mmhg) # D'Agostino-próba a ferdeségre##
## D'Agostino skewness test
##
## data: d$mmhg
## skew = -0.080088, z = -0.268000, p-value = 0.7887
## alternative hypothesis: data have a skewnessanscombe.test(x = d$mmhg) # Anscombe-Glynn-próba a csúcsosságra##
## Anscombe-Glynn kurtosis test
##
## data: d$mmhg
## kurt = 1.9412, z = -2.7803, p-value = 0.00543
## alternative hypothesis: kurtosis is not equal to 3A két alakmutatót szintén együttesen veszi figyelembe a Jarque–Bera-próba, így alkalmas a \(H_0:\) a változó normális eloszlású a populációban és a \(H_1:\) a változó nem normális eloszlású a populációban hipotézisek tesztelésére. Ez utóbbi próba a nullhipotézis igaz volta mellett \(\chi^2(n-2)\) eloszlású próbastatisztikára vezet, így nagyobb mintaelemszám esetén érdemes használni. A próba alkalmazásának feltétele, hogy a vizsgált \(X\) változó folytonos eloszlású legyen. A Jarque–Bera-próba próbastatisztikája: \(LM=n\left(\frac{\gamma_1^2}{6}+\frac{(\gamma_2-3)^2}{24}\right)\)
library(DescTools)
JarqueBeraTest(d$mmhg, robust = FALSE) # Jarque-Bera-próba végrehajtása##
## Jarque Bera Test
##
## data: d$mmhg
## X-squared = 2.6279, df = 2, p-value = 0.2688jbTest(d$mmhg)##
## Attaching package: 'akima'## The following object is masked from 'package:tolerance':
##
## interp##
## Title:
## Jarque - Bera Normality Test
##
## Test Results:
## PARAMETER:
## Sample Size: 55
## STATISTIC:
## LM: 2.628
## ALM: 2.782
## P VALUE:
## LM p-value: 0.145
## ALM p-value: 0.171
## Asymptotic: 0.269
##
## Description:
## Wed May 18 23:11:28 2016 by user: kali2# A Jarque--Bera-próba kézi számítása
n <- sum(!is.na(d$mmhg)) # mintaelemszám
m.1 <- sum(d$mmhg)/n # 1. momentum, mintaátlag
m.2 <- sum((d$mmhg-m.1)^2)/n # 2. momentum, variancia
m.3 <- sum((d$mmhg-m.1)^3)/n # 3. momentum
m.4 <- sum((d$mmhg-m.1)^4)/n # 4. momentum
gamma.1 <- m.3/m.2^(3/2) # ferdeségi mutató, gamma1
gamma.2 <- (m.4/m.2^2) # csúcsossági mutató, gamma2
LM <- n * (gamma.1^2 / 6 + (gamma.2 - 3)^2 / 24) # Jarque--Bera-próba LM próbastatisztikája
LM # próbastatisztika kiírása## [1] 2.627909pchisq(q = LM, df = 2, lower.tail = FALSE) # p érték az LM próbastatisztikához## [1] 0.2687552A Jarque–Bera-próba outputjából kiolvasható, hogy kevésbé érzékeny, mint a D’Agostino–Pearson-próba, hiszen a p érték alapján a nullhipotézist megtartjuk, azaz a változó normalitása a Jarque–Bera-próba eredménye szerint nem sérül (\(\chi^2=2,6279\); \(p=0,2688\)).
5.4 Shapiro–Wilk-próba
Ha azt szeretnénk vizsgálni, hogy a minta származhat-e normális eloszlású változóból, akkor a Shapiro–Wilk-próbát használhatjuk. (A normalitás ellenőrzésére a leginkább javasolt próba, ha kevés azonos érték fordul elő a mintában. Az SPSS hagyományosan ezt és a Kolmogorov–Szmirnov-próba Lilliefors változatát outputolja, és az a hagyomány alakult ki, hogy \(n<50\) esetén a Shapiro–Wilk-próba, ellenkező esetben a Lilliefors-próba eredményét tekintették.)
Alkalmazási feltételek:
- legyen \(X\) intervallum vagy arány skálájú változó.
Null- és ellenhipotézisek:
- Nullhipotézis:
- \(H_0:\) a változó normális eloszlású a populációban
- Ellenhipotézis:
- \(H_1:\) a változó nem normális eloszlású a populációban
Próbastatisztika:
A próbastatisztikát \(W\)-vel jelöljük, amelyre teljesül, hogy \(0<W \leq 1\). A \(W\) eloszlása függ a minta méretétől (\(3 \leq n \leq 5000\)).
5.4.1 Példa: a szisztolés vérnyomás normalitása
Vizsgáljuk meg a szisztolés vérnyomás normalitását a teljes mintában, és az egészséges és magas vérnyomásban szenvedők körében is.
Forrás: (Taeger and Kuhnt 2014)
Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.
d <- read.table(file = textConnection("
no,status,mmhg
1,0,120
2,0,115
3,0,94
4,0,118
5,0,111
6,0,102
7,0,102
8,0,131
9,0,104
10,0,107
11,0,115
12,0,139
13,0,115
14,0,113
15,0,114
16,0,105
17,0,115
18,0,134
19,0,109
20,0,109
21,0,93
22,0,118
23,0,109
24,0,106
25,0,125
26,1,150
27,1,142
28,1,119
29,1,127
30,1,141
31,1,149
32,1,144
33,1,142
34,1,149
35,1,161
36,1,143
37,1,140
38,1,148
39,1,149
40,1,141
41,1,146
42,1,159
43,1,152
44,1,135
45,1,134
46,1,161
47,1,130
48,1,125
49,1,141
50,1,148
51,1,153
52,1,145
53,1,137
54,1,147
55,1,169
"), header=T, sep=",")
d$no <- factor(d$no)
d$status <- factor(d$status, labels=c("egészséges", "magas vérnyomás"))
str(d) # az adattábla szerkezete## 'data.frame': 55 obs. of 3 variables:
## $ no : Factor w/ 55 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ status: Factor w/ 2 levels "egészséges","magas vérnyomás": 1 1 1 1 1 1 1 1 1 1 ...
## $ mmhg : int 120 115 94 118 111 102 102 131 104 107 ...d # a beolvasott adattábla## no status mmhg
## 1 1 egészséges 120
## 2 2 egészséges 115
## 3 3 egészséges 94
## 4 4 egészséges 118
## 5 5 egészséges 111
## 6 6 egészséges 102
## 7 7 egészséges 102
## 8 8 egészséges 131
## 9 9 egészséges 104
## 10 10 egészséges 107
## 11 11 egészséges 115
## 12 12 egészséges 139
## 13 13 egészséges 115
## 14 14 egészséges 113
## 15 15 egészséges 114
## 16 16 egészséges 105
## 17 17 egészséges 115
## 18 18 egészséges 134
## 19 19 egészséges 109
## 20 20 egészséges 109
## 21 21 egészséges 93
## 22 22 egészséges 118
## 23 23 egészséges 109
## 24 24 egészséges 106
## 25 25 egészséges 125
## 26 26 magas vérnyomás 150
## 27 27 magas vérnyomás 142
## 28 28 magas vérnyomás 119
## 29 29 magas vérnyomás 127
## 30 30 magas vérnyomás 141
## 31 31 magas vérnyomás 149
## 32 32 magas vérnyomás 144
## 33 33 magas vérnyomás 142
## 34 34 magas vérnyomás 149
## 35 35 magas vérnyomás 161
## 36 36 magas vérnyomás 143
## 37 37 magas vérnyomás 140
## 38 38 magas vérnyomás 148
## 39 39 magas vérnyomás 149
## 40 40 magas vérnyomás 141
## 41 41 magas vérnyomás 146
## 42 42 magas vérnyomás 159
## 43 43 magas vérnyomás 152
## 44 44 magas vérnyomás 135
## 45 45 magas vérnyomás 134
## 46 46 magas vérnyomás 161
## 47 47 magas vérnyomás 130
## 48 48 magas vérnyomás 125
## 49 49 magas vérnyomás 141
## 50 50 magas vérnyomás 148
## 51 51 magas vérnyomás 153
## 52 52 magas vérnyomás 145
## 53 53 magas vérnyomás 137
## 54 54 magas vérnyomás 147
## 55 55 magas vérnyomás 169A Shapiro–Wilk-próba végrehajtása. Shapiro–Wilk-próbát a shapiro.test() függvénnyel hajthatunk végre. Ha csoportokra szeretnénk kérni normalitás-vizsgálatot, akkor használjuk a by() függvényt.
shapiro.test(x = d$mmhg) # Shapiro-Wilk-próba a teljes mintára##
## Shapiro-Wilk normality test
##
## data: d$mmhg
## W = 0.96077, p-value = 0.07012by(d$mmhg, d$status, shapiro.test) # Shapiro-Wilk-próba az egyes csoportokra## d$status: egészséges
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.96052, p-value = 0.4249
##
## --------------------------------------------------------
## d$status: magas vérnyomás
##
## Shapiro-Wilk normality test
##
## data: dd[x, ]
## W = 0.97971, p-value = 0.8179A Shapiro–Wilk-próba értékelése. A Shapiro–Wilk-próba eredménye alapján azt mondhatjuk, hogy sem a teljes minta (\(W=0,96077; p=0,07\)), sem az egyes csoportok mintái (egészséges csoport: \(W=0,96052, p=0,4249\); magas vérnyomás csoport: \(W=0,979717, p=0,8179\)) nem sértik meg a normalitásra vonatkozó nullhipotézist.
Alternatívák a Shapiro–Wilk-próbára. A Shapiro–Wilk-próba egyszerűsített változata a Shapiro–Francia-próba, amelyet nagyobb mintaelemszámok esetén használhatunk.
library(DescTools)
ShapiroFranciaTest(d$mmhg) # Shapiro--Francia-próba a teljes mintára##
## Shapiro-Francia normality test
##
## data: d$mmhg
## W = 0.96724, p-value = 0.1231by(d$mmhg, d$status, ShapiroFranciaTest) # Shapiro--Francia-próba az egyes csoportokra## d$status: egészséges
##
## Shapiro-Francia normality test
##
## data: dd[x, ]
## W = 0.95938, p-value = 0.3412
##
## --------------------------------------------------------
## d$status: magas vérnyomás
##
## Shapiro-Francia normality test
##
## data: dd[x, ]
## W = 0.97552, p-value = 0.60615.5 Kolmogorov–Szmirnov-próba
Eloszlás vizsgálatára az egymintás Kolmogorov–Szmirnov-próbát is használhatjuk, amely azt vizsgálja, hogy a minta vajon adott (pl. normális, Poisson, egyenletes, exponenciális stb.) eloszlású-e. Használata nagyobb mintaelemszám (\(n \geq 50\)) esetén javasolt.
Alkalmazási feltételek:
- legyen \(X\) intervallum vagy arány skálájú változó.
Null- és ellenhipotézisek:
- Nullhipotézis:
- \(H_0:\) a változó normális eloszlású a populációban
- Ellenhipotézis:
- \(H_1:\) a változó nem normális eloszlású a populációban
Próbastatisztika:
A \(D\) próbastatisztika az elméleti és az empirikus eloszlásfüggvény közötti maximális eltérésén alapul.
5.5.1 Példa: a szisztolés vérnyomás normalitása
Vizsgáljuk meg a szisztolés vérnyomás normalitását a teljes mintában, és az egészséges és magas vérnyomásban szenvedők körében is.
Forrás: (Taeger and Kuhnt 2014)
Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.
d <- read.table(file = textConnection("
no,status,mmhg
1,0,120
2,0,115
3,0,94
4,0,118
5,0,111
6,0,102
7,0,102
8,0,131
9,0,104
10,0,107
11,0,115
12,0,139
13,0,115
14,0,113
15,0,114
16,0,105
17,0,115
18,0,134
19,0,109
20,0,109
21,0,93
22,0,118
23,0,109
24,0,106
25,0,125
26,1,150
27,1,142
28,1,119
29,1,127
30,1,141
31,1,149
32,1,144
33,1,142
34,1,149
35,1,161
36,1,143
37,1,140
38,1,148
39,1,149
40,1,141
41,1,146
42,1,159
43,1,152
44,1,135
45,1,134
46,1,161
47,1,130
48,1,125
49,1,141
50,1,148
51,1,153
52,1,145
53,1,137
54,1,147
55,1,169
"), header=T, sep=",")
d$no <- factor(d$no)
d$status <- factor(d$status, labels=c("egészséges", "magas vérnyomás"))
str(d) # az adattábla szerkezete## 'data.frame': 55 obs. of 3 variables:
## $ no : Factor w/ 55 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ status: Factor w/ 2 levels "egészséges","magas vérnyomás": 1 1 1 1 1 1 1 1 1 1 ...
## $ mmhg : int 120 115 94 118 111 102 102 131 104 107 ...d # a beolvasott adattábla## no status mmhg
## 1 1 egészséges 120
## 2 2 egészséges 115
## 3 3 egészséges 94
## 4 4 egészséges 118
## 5 5 egészséges 111
## 6 6 egészséges 102
## 7 7 egészséges 102
## 8 8 egészséges 131
## 9 9 egészséges 104
## 10 10 egészséges 107
## 11 11 egészséges 115
## 12 12 egészséges 139
## 13 13 egészséges 115
## 14 14 egészséges 113
## 15 15 egészséges 114
## 16 16 egészséges 105
## 17 17 egészséges 115
## 18 18 egészséges 134
## 19 19 egészséges 109
## 20 20 egészséges 109
## 21 21 egészséges 93
## 22 22 egészséges 118
## 23 23 egészséges 109
## 24 24 egészséges 106
## 25 25 egészséges 125
## 26 26 magas vérnyomás 150
## 27 27 magas vérnyomás 142
## 28 28 magas vérnyomás 119
## 29 29 magas vérnyomás 127
## 30 30 magas vérnyomás 141
## 31 31 magas vérnyomás 149
## 32 32 magas vérnyomás 144
## 33 33 magas vérnyomás 142
## 34 34 magas vérnyomás 149
## 35 35 magas vérnyomás 161
## 36 36 magas vérnyomás 143
## 37 37 magas vérnyomás 140
## 38 38 magas vérnyomás 148
## 39 39 magas vérnyomás 149
## 40 40 magas vérnyomás 141
## 41 41 magas vérnyomás 146
## 42 42 magas vérnyomás 159
## 43 43 magas vérnyomás 152
## 44 44 magas vérnyomás 135
## 45 45 magas vérnyomás 134
## 46 46 magas vérnyomás 161
## 47 47 magas vérnyomás 130
## 48 48 magas vérnyomás 125
## 49 49 magas vérnyomás 141
## 50 50 magas vérnyomás 148
## 51 51 magas vérnyomás 153
## 52 52 magas vérnyomás 145
## 53 53 magas vérnyomás 137
## 54 54 magas vérnyomás 147
## 55 55 magas vérnyomás 169A Kolmogorov–Szmirnov-próba végrehajtása. Kolmogorov–Szmirnov-próbát a ks.test() függvénnyel hajthatunk végre. Ha csoportokra szeretnénk kérni normalitás-vizsgálatot, akkor használjuk a by() függvényt. Ha a ks.test() függvényt a normalitás ellenőrzésére használjuk, akkor az x = d$mmhg minta után az elméleti eloszlásfüggvényt ("pnorm") és azok paramétereit (mean= és sd=) is meg kell adnunk.
# Kolmogorov--Szmirnov-próba a teljes mintára
ks.test(x = d$mmhg, "pnorm", mean = mean(d$mmhg), sd = sd(d$mmhg))##
## One-sample Kolmogorov-Smirnov test
##
## data: d$mmhg
## D = 0.11725, p-value = 0.4361
## alternative hypothesis: two-sided# Kolmogorov--Szmirnov-próba az egyes csoportokra
by(d$mmhg, d$status, function(x) ks.test(x, "pnorm", mean = mean(x), sd = sd(x)))## d$status: egészséges
##
## One-sample Kolmogorov-Smirnov test
##
## data: x
## D = 0.14603, p-value = 0.6606
## alternative hypothesis: two-sided
##
## --------------------------------------------------------
## d$status: magas vérnyomás
##
## One-sample Kolmogorov-Smirnov test
##
## data: x
## D = 0.11729, p-value = 0.8036
## alternative hypothesis: two-sidedAlternatívák a Kolmogorov–Szmirnov-próbára. A Kolmogorov–Szmirnov-próbát módosított változata a Lilliefors-próba, amelyet a DescTools csomag LillieTest() függvényével tudjuk végrehajtani. A függvény egyetlen argumentuma a minta.
library(DescTools)
LillieTest(x = d$mmhg) # Lilliefors-próba a teljes mintára##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: d$mmhg
## D = 0.11725, p-value = 0.05723by(d$mmhg, d$status, LillieTest) # Lilliefors-próba az egyes csoportokra## d$status: egészséges
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: dd[x, ]
## D = 0.14603, p-value = 0.1837
##
## --------------------------------------------------------
## d$status: magas vérnyomás
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: dd[x, ]
## D = 0.11729, p-value = 0.3645