4 Varianciára vonatkozó próbák

Ebben a fejezetben normális eloszlású változók varianciájára (és egyben szórására) vonatkozó hipotézisek tesztelését ismertetjük. Az egy változó varianciájára vonatkozó \(\sigma^2=\sigma_0^2\) hipozézissel kezdjük, amely khi-négyzet eloszlású próbastatisztikára vezet. A \(\sigma_1^2=\sigma_2^2\) hipotézis vizsgálatát páros és független mintás esetben is ismertetjük, majd bemutatjuk a 3 vagy annál több független változóhoz tartozó \(\sigma_1^2=\sigma_2^2= \dots = \sigma_k^2\) hipotézis vizsgálatát.

4.1 Khi-négyzet próba egy variancia vizsgálatára

A normális eloszlású változó varianciájára vonatkozó hipotézist teszteljük. Azt a kérdést tesszük fel, hogy a variancia (vagy szórás) eltér-e egy előre definiált (hipotetikus) értéktől.

Alkalmazási feltételek:

  • \(X\) normális eloszlású a populációban

Null- és ellenhipotézisek:

  • Nullhipotézis:
    • \(H_0: \sigma^2=\sigma^2_0\)
  • Ellenhipotézis:
    • \(H_1: \sigma^2\neq\sigma^2_0\) (kétoldali ellenhipotézis) vagy
    • \(H_1: \sigma^2<\sigma^2_0\) (egyoldali, kisebb (bal oldali) ellenhipotézis) vagy
    • \(H_1: \sigma^2>\sigma^2_0\) (egyoldali, nagyobb (jobb oldali) ellenhipotézis)

Próbastatisztika:

\[\chi^2=\frac{(n-1)s^2}{\sigma_0^2} \sim \chi^2(n-1)\]

4.1.1 Példa: szisztolés vérnyomás varianciája

Vizsgáljuk meg, hogy a szisztolés vérnyomás varianciája lehet-e 400 a populációban (vagyis a szórás lehet-e 20 a populációban) a minta alapján!
Forrás: (Taeger and Kuhnt 2014)

A feladat szövegéből a szisztolés vérnyomás varianciájára (és szórására) a következő null- és ellenhipotézis olvasható ki:

  • \(H_0: \sigma^2=400\) (\(H_0: \sigma=20\))
  • \(H_1: \sigma^2\neq400\) (\(H_1: \sigma\neq20\))

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  169

Az egy variancia vizsgálatára vonatkozó khi-négyzet próba végrehajtása. Az egy varianciára vonatkozó próbát a TeachingDemos csomag sigma.test() függvényével hajtjuk végre.

library(TeachingDemos)
sigma.test(x = d$mmhg, sigma = 20)
## 
##  One sample Chi-squared test for variance
## 
## data:  d$mmhg
## X-squared = 49.595, df = 54, p-value = 0.7104
## alternative hypothesis: true variance is not equal to 400
## 95 percent confidence interval:
##  260.3684 557.4611
## sample estimates:
## var of d$mmhg 
##      367.3704

Az egy variancia vizsgálatára vonatkozó khi-négyzet próba értékelése. Az egy variancia vizsgálatára vonatkozó khi-négyzet próba eredménye alapján azt mondhatjuk, hogy nincs elegendő bizonyítékunk arra, hogy a szisztolés vérnyomás varianciája eltér 400-tól (vagy a szórás a 20-tól) (\(\chi^2=49,595; p=0,7104\)).

Az egy variancia vizsgálatára vonatkozó khi-négyzet próba kézi számolása. A próba outputjában szereplő értékeket magunk is kiszámolhatjuk:

n        <- sum(!is.na(d$mmhg))                          # mintaeleszám
s        <- sd(x = d$mmhg, na.rm = T)                    # minta szórása
s.0      <- 20                                           # hipotetikus szórás
khi.2    <- (n-1)*s^2/s.0^2                              # próbastatisztika értéke
khi.2                                                    # próbastatisztika értékének kiírása
## [1] 49.595
2*pchisq(q = khi.2, df = n - 1)                          # p érték kétoldali próba esetén
## [1] 0.7103942

Konfidenciaintervallum kézi számolása az egy variancia vizsgálatára vonatkozó khi-négyzet próbához. A próba outputjában szereplő konfidenciaintervallum határait kézi számolással is meghatározzuk. A varianciára vonatkozó \(\left(1-\alpha\right)\) megbízhatóságú konfidenciaintervallum:

\[\left]\frac{(n-1)s^2}{\chi_{\alpha/2}^2},\frac{(n-1)s^2}{\chi_{1-\alpha/2}^2}\right[\]

A konfidenciaintervallum határait a DescTools csomag VarCI() függvényével és kézzel is kiszámolhatjuk:

library(DescTools)
VarCI(d$mmhg)
##      var   lwr.ci   upr.ci 
## 367.3704 260.3684 557.4611
n        <- sum(!is.na(d$mmhg))                                 # mintaeleszám
s        <- sd(x = d$mmhg, na.rm = T)                           # minta szórása
alfa     <- 0.05                                                # alfa=1-0.95  
chi.alfa.1 <- qchisq(p = alfa/2, df = n - 1, lower.tail = F)    # a bal oldali határhoz
chi.alfa.2 <- qchisq(p = 1-alfa/2, df = n - 1, lower.tail = F)  # a jobb oldali határhoz
(n-1)*s^2*1/c(chi.alfa.1, chi.alfa.2)                           # az intervallum
## [1] 260.3684 557.4611

A hatás nagysága az egy variancia vizsgálatára vonatkozó khi-négyzet próbához.

s          <- sd(x = d$mmhg, na.rm = T)                         # minta szórása
s.0        <- 20                                                # hipotetikus szórás
delta      <- s.0^2/s^2                                         # delta
alfa       <- 0.05                                              # alfa=1-0.95  
chi.alfa.1 <- qchisq(p = alfa/2, df = n - 1, lower.tail = F)    # a bal oldali határhoz
chi.alfa.2 <- qchisq(p = 1-alfa/2, df = n - 1, lower.tail = F)  # a jobb oldali határhoz
pchisq(q = delta*chi.alfa.2, df = n - 1, lower.tail = F) -      # beta
  pchisq(q = delta*chi.alfa.1, df = n - 1, lower.tail = F)
## [1] 0.9347617

4.2 Pitman–Morgan-próba

A Pitman–Morgan-próbát varianciák egyezésének vizsgálatára használjuk összetartozó minták esetén.

Alkalmazási feltételek:

  • véletlen, összetartozó (páros) minta: \((X_1, X_2)\)
  • az \(X_1\) és \(X_2\) normális eloszlású (elegendő, ha a \(X_1-X_2\) különbségváltozó normális eloszlású)

Null- és ellenhipotézisek:

  • Nullhipotézis:
    • \(H_0: \sigma_1^2=\sigma_2^2\)
  • Ellenhipotézis:
    • \(H_1: \sigma_1^2\neq\sigma_2^2\) (kétoldali ellenhipotézis) vagy
    • \(H_1: \sigma_1^2<\sigma_2^2\) (egyoldali, kisebb (bal oldali) ellenhipotézis) vagy
    • \(H_1: \sigma_1^2>\sigma_2^2\) (egyoldali, nagyobb (jobb oldali) ellenhipotézis)

Próbastatisztika:

\[T=\frac{\sqrt{n-2}(s_1^2-s_2^2)}{\sqrt{4(1-r^2)s_1^2s_2^2}} \sim t(n-2)\]

4.2.1 Példa: intelligencia varianciája

Vizsgáljuk meg azt a hipotézist, hogy a tréning előtt és tréning után az intelligencia varianciája megegyezik egymással!
Forrás: (Taeger and Kuhnt 2014)

A példa alapján az összetartozó minták mögötti változók varianciájának egyezését vizsgáljuk:

  • \(H_0: \sigma_1^2=\sigma_2^2\)
  • \(H_1: \sigma_1^2\neq\sigma_2^2\)

Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.

d <- read.table(file = textConnection("
no,IQ1,IQ2
1,127,137
2,98,108
3,105,115
4,83,93
5,133,143
6,90,100
7,107,117
8,98,108
9,91,101
10,100,110
11,88,98
12,96,106
13,110,120
14,87,97
15,88,98
16,88,100
17,105,115
18,95,111
19,79,89
20,106,116
"), header=T, sep=",")
d$no <- factor(d$no)
str(d)            # az adattábla szerkezete
## 'data.frame':    20 obs. of  3 variables:
##  $ no : Factor w/ 20 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ IQ1: int  127 98 105 83 133 90 107 98 91 100 ...
##  $ IQ2: int  137 108 115 93 143 100 117 108 101 110 ...
d                 # a beolvasott adattábla
##    no IQ1 IQ2
## 1   1 127 137
## 2   2  98 108
## 3   3 105 115
## 4   4  83  93
## 5   5 133 143
## 6   6  90 100
## 7   7 107 117
## 8   8  98 108
## 9   9  91 101
## 10 10 100 110
## 11 11  88  98
## 12 12  96 106
## 13 13 110 120
## 14 14  87  97
## 15 15  88  98
## 16 16  88 100
## 17 17 105 115
## 18 18  95 111
## 19 19  79  89
## 20 20 106 116

Az adatokat készítsük elő a páros mintás elemzésre is:

library(PairedData)
d.paired <- with(d, paired(IQ1, IQ2))    # adatelőkészítés a PairedData csomag használatához
d.paired                                 # a páros minta
## Object of class "paired"
##    IQ1 IQ2
## 1  127 137
## 2   98 108
## 3  105 115
## 4   83  93
## 5  133 143
## 6   90 100
## 7  107 117
## 8   98 108
## 9   91 101
## 10 100 110
## 11  88  98
## 12  96 106
## 13 110 120
## 14  87  97
## 15  88  98
## 16  88 100
## 17 105 115
## 18  95 111
## 19  79  89
## 20 106 116

Páros Pitman–Morgan-próba végrehajtása. Páros Pitman-Morgan próbát a PairedData csomag var.test() függvényével végzünk:

library(PairedData)
var.test(d.paired)                                    # Páros Pitman--Morgan-próba
## 
##  Paired Pitman-Morgan test
## 
## data:  IQ1 and IQ2
## t = 0.2967, df = 18, p-value = 0.7701
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.9174365 1.1213427
## sample estimates:
## variance of x variance of y 
##      188.4316      185.7789

Páros Pitman–Morgan-próba értékelése. Az páros mintán alapuló Pitman–Morgan-próba eredménye alapján azt mondhatjuk, hogy nincs elegendő bizonyítékunk arra, hogy az intelligencia pontszámok varianciája eltér a tréning előtt és után (\(t(18)=0,2967; p=0,7701\)).

Páros Páros Pitman–Morgan-próba kézi számolása.

szoras.1 <- sd(d$IQ1)                                 # 1. minta szórása
szoras.2 <- sd(d$IQ2)                                 # 2. minta szórása 
r        <- cor(x = d$IQ1, y = d$IQ2)                 # korrelációs együttható
n        <- sum(!is.na(d$IQ1))                        # mintaelemszám  
t <- (sqrt(n-2)*(szoras.1^2-szoras.2^2))/sqrt(4*(1-r^2)*szoras.1^2*szoras.2^2) # próbastatisztika értéke
t                                                     # próbastatisztika értékének kiírása
## [1] 0.2966984
2*(1-pt(q = abs(t),df = n - 2))                       # p érték kétoldali próbára
## [1] 0.7700932

4.3 F-próba

Az F-próbát varianciák egyezésének vizsgálatára használjuk független minták esetén.

Alkalmazási feltételek:

  • két véletlen, független minta: \(X_1\), \(X_2\)
  • az \(X_1\) és \(X_2\) normális eloszlású

Null- és ellenhipotézisek:

  • Nullhipotézis:
    • \(H_0: \sigma_1^2=\sigma_2^2\)
  • Ellenhipotézis:
    • \(H_1: \sigma_1^2\neq\sigma_2^2\) (kétoldali ellenhipotézis) vagy
    • \(H_1: \sigma_1^2<\sigma_2^2\) (egyoldali, kisebb (bal oldali) ellenhipotézis) vagy
    • \(H_1: \sigma_1^2>\sigma_2^2\) (egyoldali, nagyobb (jobb oldali) ellenhipotézis)

Próbastatisztika:

\[F=\frac{s_1^2}{s_2^2} \sim F(n_1-1,n_2-1)\]

4.3.1 Példa: szisztolés vérnyomás varianciája két csoportban

Vizsgáljuk meg azt a hipotézist, hogy a szisztolés vérnyomás varianciája megegyezik egészségesek és magas vérnyomásban szenvedők körében!
Forrás: (Taeger and Kuhnt 2014)

A példa alapján független minták mögötti változók populációbeli varianciájának egyezését vizsgáljuk:

  • \(H_0: \sigma_1^2=\sigma_2^2\)
  • \(H_1: \sigma_1^2\neq\sigma_2^2\)

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  169

F-próba végrehajtása. Az F-próba végrehajtásához az alapértelmezetten betöltődő stats csomag var.test() függvénnyel végezhető el. Az első argumentumában megadjuk a formulát (formula=mmhg~status), amely definiálja független változójában a két csoportot (status) és a függő változót (mmhg).

var.test(formula = mmhg ~ status, data = d)
## 
##  F test to compare two variances
## 
## data:  mmhg by status
## F = 1.0363, num df = 24, denom df = 29, p-value = 0.918
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4811235 2.2980313
## sample estimates:
## ratio of variances 
##           1.036343

Az F-próba értékelése. Az F-próba eredménye alapján azt mondhatjuk, hogy nincs elegendő bizonyítékunk arra, hogy a szisztolés vérnyomás varianciája eltér az egészségesek és a magas vérnyomásban szenvedők két csoportjában (\(F(24, 29)=1,0363; p=0,918\)).

Az F-próba kézi számolása.

s.1      <- sd(d$mmhg[d$status=="egészséges"])                         # 1. minta szórása
s.2      <- sd(d$mmhg[d$status=="magas vérnyomás"])                    # 2. minta szórása
n.1      <- sum(!is.na(d$mmhg[d$status=="egészséges"]))                # 1. minta elemszáma
n.2      <- sum(!is.na(d$mmhg[d$status=="magas vérnyomás"]))           # 2. minta elemszáma
F.stat   <- s.1^2/s.2^2                                                # próbastatisztika értéke
F.stat                                                                 # próbastatisztika értékének kiírása
## [1] 1.036343
2*(1-pf(q = F.stat, df1 = n.1 - 1, df2 = n.2 - 1))                     # p érték kétoldali próbára
## [1] 0.9179638

4.4 Levene-próba

Ha két vagy több ismeretlen populációbeli varianciát szeretnénk összehasonlítani \(k\) független minta alapján, úgy hogy nem vagyunk biztosak a változók normális eloszlásában, akkor a Levene-próbát használhatjuk.

Alkalmazási feltételek:

  • \(k\) független minta, melyek intervallum vagy arány skálájú változóktól származnak

Null- és ellenhipotézisek:

  • Nullhipotézis:
    • \(H_0: \sigma_1^2=\sigma_2^2= \dots = \sigma_k^2\)
  • Ellenhipotézis:
    • \(H_1:\) nem mindegyik variancia azonos

Próbastatisztika:

A Levene-próba a csoportátlagoktól (\(\bar{Y_i}\)) mért abszolút eltéréseken (\(Z_{i,j}\)) alapuló varianciaelemzésen alapul. Amennyiben nem a mintaátlagot vesszük alapul, hanem a mediánt (\(\tilde{Y_i}\)), akkor a varianciaelemzés a Brown–Forsyth-próbára vezet. A próbastatisztika tehát a nullhipotézis igaz volta esetén \(F(k-1, N-k)\) eloszlású, ahol \(k\) a csoportok száma, az \(N\) pedig az összmintaelemszám.

\[Z_{i,j}= \begin{cases} \left|Y_{i,j}-\bar{Y_i}\right|, & \text{ahol}\ \bar{Y_i}\ \text{az}\ i. \text{mintaátlag} \\ \left|Y_{i,j}-\tilde{Y_i}\right|, & \text{ahol}\ \tilde{Y_i}\ \text{az}\ i. \text{minta mediánja} \end{cases}\]

\[F=\frac{(N-k)\sum_{i=1}^kn_i(\bar{Z_i}-\bar{Z})^2}{(k-1)\sum_{i=1}^k\sum_{j=1}^{n_i}(Z_{i,j}-\bar{Z_i})^2} \sim F(k-1, N-k)\]

4.4.1 Példa: hangulatpontszámok varianciája több csoportban.

Vizsgáljuk meg, hogy a hangulatpontszámok varianciája eltér-e egymástól a populációban az alkalmazott három kezelés esetén!
Forrás: http://health.adelaide.edu.au/psychology/ccs/teaching/lsr/

Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.

d <- read.table(file = textConnection("
     drug    therapy mood.gain
  placebo no.therapy       0.5
  placebo no.therapy       0.3
  placebo no.therapy       0.1
 anxifree no.therapy       0.6
 anxifree no.therapy       0.4
 anxifree no.therapy       0.2
 joyzepam no.therapy       1.4
 joyzepam no.therapy       1.7
 joyzepam no.therapy       1.3
  placebo        CBT       0.6
  placebo        CBT       0.9
  placebo        CBT       0.3
 anxifree        CBT       1.1
 anxifree        CBT       0.8
 anxifree        CBT       1.2
 joyzepam        CBT       1.8
 joyzepam        CBT       1.3
 joyzepam        CBT       1.4
"), header=T, sep="")
str(d)                       # az adattábla szerkezete
## 'data.frame':    18 obs. of  3 variables:
##  $ drug     : Factor w/ 3 levels "anxifree","joyzepam",..: 3 3 3 1 1 1 2 2 2 3 ...
##  $ therapy  : Factor w/ 2 levels "CBT","no.therapy": 2 2 2 2 2 2 2 2 2 1 ...
##  $ mood.gain: num  0.5 0.3 0.1 0.6 0.4 0.2 1.4 1.7 1.3 0.6 ...
d                            # a beolvasott adattábla
##        drug    therapy mood.gain
## 1   placebo no.therapy       0.5
## 2   placebo no.therapy       0.3
## 3   placebo no.therapy       0.1
## 4  anxifree no.therapy       0.6
## 5  anxifree no.therapy       0.4
## 6  anxifree no.therapy       0.2
## 7  joyzepam no.therapy       1.4
## 8  joyzepam no.therapy       1.7
## 9  joyzepam no.therapy       1.3
## 10  placebo        CBT       0.6
## 11  placebo        CBT       0.9
## 12  placebo        CBT       0.3
## 13 anxifree        CBT       1.1
## 14 anxifree        CBT       0.8
## 15 anxifree        CBT       1.2
## 16 joyzepam        CBT       1.8
## 17 joyzepam        CBT       1.3
## 18 joyzepam        CBT       1.4

A Levene-próba végrehajtása. Levene-próbát a DescTools csomag LeveneTest() függvényével is végrehajthatjuk. A formula argumentum a függő változó (mood.gain) és a faktor (drug) megadásával definiálja a csoportokat. A Levene-próba a csoportközepek kiszámításához a mintaátlagot használja, így a center=mean argumentumot is megadjuk. Amennyiben ezt nem tesszük meg, akkor a center=median alapértelmezett értékkel hívódik a függvény, így Brown–Forsyth kerül végrehajtásra.

library(DescTools)
LeveneTest(mood.gain ~ drug, data = d, center = mean)                    # Levene-próba végrehajtása
## Levene's Test for Homogeneity of Variance (center = mean)
##       Df F value Pr(>F)
## group  2  1.4497 0.2657
##       15

A Levene-próba értékelése. A Levene-próba a hangulatpontszámok 3 csoportjában az azonos varianciára vonatkozó nullhipotézist támogatja (\(F(2,15)=1,4497\); \(p=0,2657\)).

A Levene-próba kézi kiszámítása. A fenti outputban szereplő numerikus értékeket magunk is meghatározhatjuk:

Y.i    <- tapply(d$mood.gain, d$drug, mean, na.rm = T)       # csoportok mintaátlaga
Zd     <- abs(d$mood.gain-Y.i[d$drug])                       # a mintaelemek abszolút eltérései a saját csoportátlaguktól
Z.i    <- tapply(Zd, d$drug, mean, na.rm = T)                # Zd csoportátlagai
Z      <- mean(Zd)                                           # összátlag
N      <- length(d$mood.gain)                                # összmintaelemszám
k      <- nlevels(d$drug)                                    # csoportok száma
n.i    <- tapply(d$mood.gain, d$drug, length)                # csoportok mintaelemszáma
F.test <- (N - k) * sum(n.i*(Z.i-Z)^2) / ((k - 1) * sum((Zd-Z.i[d$drug])^2))
F.test                                                       # próbastatisztika kiírása
## [1] 1.449739
pf(q = F.test, df1 = k - 1, df2 = N - k, lower.tail = FALSE) # p érték
## [1] 0.2656941
# elvégezhetjük a variancielemzés módszerével is
Y.i    <- tapply(d$mood.gain, d$drug, mean, na.rm = T)       # csoportok mintaátlaga
Zd     <- abs(d$mood.gain-Y.i[d$drug])                       # a mintaelemek abszolút eltérései a saját csoportátlaguktól
summary(aov(Zd ~ d$drug))
##             Df Sum Sq Mean Sq F value Pr(>F)
## d$drug       2 0.0616 0.03080    1.45  0.266
## Residuals   15 0.3187 0.02125
# Brown-Forsyth próba a variancielemzés módszerével
Y.i    <- tapply(d$mood.gain, d$drug, median, na.rm = T)     # csoportok mediánja
Zd     <- abs(d$mood.gain-Y.i[d$drug])                       # a mintaelemek abszolút eltérései a saját csoportmediánjuktól
summary(aov(Zd ~ d$drug))
##             Df Sum Sq Mean Sq F value Pr(>F)
## d$drug       2 0.0844 0.04222   1.467  0.262
## Residuals   15 0.4317 0.02878

Alternatívák Levene-próbára: a Brown–Forsyth-próba. Ha a LeveneTest() függvényt a center=median argumentummal hívjuk, akkor Brown–Forsyth-próbát hajtunk végre, amely a Levene-próba robosztus változatának is tekinthető.

library(DescTools)
LeveneTest(mood.gain ~ drug, data = d, center = median)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.4672 0.2618
##       15

Brown–Forsyth-próbát a HH csomag hov() függvényével is végrehajthatunk. A hovPlot() függvénnyel grafikus ábrázolást is rendelhetünk a próbához.

library(HH)
hov(mood.gain ~ drug, data = d)
## 
##  hov: Brown-Forsyth
## 
## data:  mood.gain
## F = 1.4672, df:drug = 2, df:Residuals = 15, p-value = 0.2618
## alternative hypothesis: variances are not identical
hovPlot(mood.gain ~ drug, data = d)

4.5 Bartlett-próba

Ha két vagy több ismeretlen populációbeli varianciát szeretnénk összehasonlítani \(k\) független minta alapján, úgy hogy a változók normális eloszlásúak, akkor Bartlett-próbát használunk.

Alkalmazási feltételek:

  • \(k\) független minta
  • normális eloszlásúak a változók

Null- és ellenhipotézisek:

  • Nullhipotézis:
    • \(H_0: \sigma_1^2=\sigma_2^2= \dots = \sigma_k^2\)
  • Ellenhipotézis:
    • \(H_1:\) nem mindegyik variancia azonos

Próbastatisztika:

A próbastatisztika \(\chi^2(k-1)\) eloszlású a nullhipotézis igaz volta esetén. Ha a \(k\) független minta \(n_i\) elemszámmal és \(s_i\) szórással rendelkezik és a pooled variancia \(s_p^2\), valamint az összeminatelemszám \(N\), akkor kiszámítása:

\[s_p^2=\frac{1}{N-k}\sum_{i=1}^k (n_i-1)s_i^2\] \[\chi^2=\frac{(N-k)\ln(s_p^2)-\sum_{i=1}^k(n_i-1)\ln(s_i^2)}{1+\frac{1}{3(k-1)}\left(\sum_{i=1}^k\left(\frac{1}{n_i-1}\right)-\frac{1}{N-k}\right)} \sim \chi^2(k-1),\]

4.5.1 Példa: hangulatpontszámok varianciája több csoportban.

Vizsgáljuk meg, hogy a hangulatpontszámok varianciája eltér-e egymástól a populációban az alkalmazott három kezelés esetén!
Forrás: http://health.adelaide.edu.au/psychology/ccs/teaching/lsr/

Adatok beolvasása. Végezzük el az adatok beolvasását egy d adattáblába.

d <- read.table(file = textConnection("
     drug    therapy mood.gain
  placebo no.therapy       0.5
  placebo no.therapy       0.3
  placebo no.therapy       0.1
 anxifree no.therapy       0.6
 anxifree no.therapy       0.4
 anxifree no.therapy       0.2
 joyzepam no.therapy       1.4
 joyzepam no.therapy       1.7
 joyzepam no.therapy       1.3
  placebo        CBT       0.6
  placebo        CBT       0.9
  placebo        CBT       0.3
 anxifree        CBT       1.1
 anxifree        CBT       0.8
 anxifree        CBT       1.2
 joyzepam        CBT       1.8
 joyzepam        CBT       1.3
 joyzepam        CBT       1.4
"), header=T, sep="")
str(d)                       # az adattábla szerkezete
## 'data.frame':    18 obs. of  3 variables:
##  $ drug     : Factor w/ 3 levels "anxifree","joyzepam",..: 3 3 3 1 1 1 2 2 2 3 ...
##  $ therapy  : Factor w/ 2 levels "CBT","no.therapy": 2 2 2 2 2 2 2 2 2 1 ...
##  $ mood.gain: num  0.5 0.3 0.1 0.6 0.4 0.2 1.4 1.7 1.3 0.6 ...
d                            # a beolvasott adattábla
##        drug    therapy mood.gain
## 1   placebo no.therapy       0.5
## 2   placebo no.therapy       0.3
## 3   placebo no.therapy       0.1
## 4  anxifree no.therapy       0.6
## 5  anxifree no.therapy       0.4
## 6  anxifree no.therapy       0.2
## 7  joyzepam no.therapy       1.4
## 8  joyzepam no.therapy       1.7
## 9  joyzepam no.therapy       1.3
## 10  placebo        CBT       0.6
## 11  placebo        CBT       0.9
## 12  placebo        CBT       0.3
## 13 anxifree        CBT       1.1
## 14 anxifree        CBT       0.8
## 15 anxifree        CBT       1.2
## 16 joyzepam        CBT       1.8
## 17 joyzepam        CBT       1.3
## 18 joyzepam        CBT       1.4

A Bartlett-próba végrehajtása. Bartlett-próba a bartlett.test() függvénnyel hajtható végre. A formula=mood.gain~drug argumentum definiálja a függő változót (mood.gain) és a csoportokat a drug faktor segítségével.

bartlett.test(formula = mood.gain ~ drug, data = d)                             # Bartlett-próba végrehajtása
## 
##  Bartlett test of homogeneity of variances
## 
## data:  mood.gain by drug
## Bartlett's K-squared = 1.6761, df = 2, p-value = 0.4326

A Bartlett-próba értékelése. A varianciák (szórások) egyezésére vonatkozó nullhipotézist megtartjuk (\(\chi^2(2)=1,6761; p=0,4326\)).

A Bartlett-próba kézi számítása. A fenti outputban szereplő numerikus értékeket magunk is kiszámíthatjuk:

k     <- nlevels(d$drug)                                                         # csoportok száma
n.i   <- tapply(d$mood.gain, d$drug, function(x) sum(!is.na(x)))                 # csoportok mintaelemszáma
N     <- sum(n.i)                                                                # összmintaelemszám
s.i.2 <- tapply(d$mood.gain, d$drug, var)                                        # csoportok szórásnégyzete
s.p.2 <- (1 / (N - k)) * sum((n.i - 1) *  s.i.2)                                 # pooled variancia

Chi2 <- ((N - k) * log(s.p.2) - sum((n.i-1)*log(s.i.2)))/ (1 + 1/(3*(k-1))*(sum(1/(n.i-1))-(1/(N - k))) )
Chi2                                                                             # próbastatisztika értéke
## [1] 1.676109
pchisq(q = Chi2, df = k - 1, lower.tail = FALSE)                                 # p érték
## [1] 0.4325513

Statisztikai mutatók a varianciára és a szórásra. A csoportokban a hangulatpontszám varianciáját és szórását magunk is meghatározhatjuk. A DescTools csomag VarCI() függvénye a populációbeli varianciára ad pont- és intervallumbecslést. A csoportokra vonatkozó mutatókat a by() segítségével íratjuk ki. A psych csomag describeBy() függvénye a szórások meghatározásának egyik módját jelenti a három csoportban.

library(DescTools)
by(d$mood.gain, d$drug, VarCI, conf.level = 0.95)
## d$drug: anxifree
##        var     lwr.ci     upr.ci 
## 0.15366667 0.05987401 0.92435346 
## -------------------------------------------------------- 
## d$drug: joyzepam
##        var     lwr.ci     upr.ci 
## 0.04566667 0.01779336 0.27469940 
## -------------------------------------------------------- 
## d$drug: placebo
##        var     lwr.ci     upr.ci 
## 0.07900000 0.03078121 0.47520991
library(psych)
describeBy(x = d$mood.gain, group = d$drug)
## group: anxifree
##   vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## 1    1 6 0.72 0.39    0.7    0.72 0.52 0.2 1.2     1    0    -1.89 0.16
## -------------------------------------------------------- 
## group: joyzepam
##   vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## 1    1 6 1.48 0.21    1.4    1.48 0.15 1.3 1.8   0.5 0.49    -1.83 0.09
## -------------------------------------------------------- 
## group: placebo
##   vars n mean   sd median trimmed  mad min max range skew kurtosis   se
## 1    1 6 0.45 0.28    0.4    0.45 0.22 0.1 0.9   0.8 0.34    -1.46 0.11

A Bartlett-próba alkalmazási feltételének ellenőrzése. A Bartlett-próba alkalmazási feltétele a normális eloszlású változók használata. Ehhez a Shapiro–Wilk-próbát használjuk.

by(d$mood.gain, d$drug, shapiro.test)                           # Normalitás ellenőrzése a 3 csoportban
## d$drug: anxifree
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.95617, p-value = 0.7898
## 
## -------------------------------------------------------- 
## d$drug: joyzepam
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.82404, p-value = 0.09563
## 
## -------------------------------------------------------- 
## d$drug: placebo
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.9614, p-value = 0.8305