Luku 4 Ryhmien odotusarvojen vertailu

Kurssilla Tilastotieteen perusteet käytiin jakauman sijaintiin liittyvä t-testi jotakuinkin seuraavanlaisessa tilanteessa. Oletetaan, että satunnaismuuttuja $X\sim N(\mu, \sigma^2)$ noudattaa normaalijakaumaa odotusarvolla $\mu$ ja varianssilla $\sigma^2$. Odotusarvo $\mu$ ja varianssi $\sigma^2$ ovat siis deterministisiä mutta tuntemattomia lukuja, jotka voimme estimoida otoksen avulla. Olkoon $\overrightarrow X = \{X_1, \ldots, X_n\}$ satunnaisotos riippumattomia havaintoja satunnaismuuttujasta $X$. Merkitään satunnaisotosta $\overrightarrow X$ vastaavaa keskiarvoa ja otosvarianssia merkinnöillä $\bar X = \frac{1}{n}\sum_{i = 1}^n X_i$ ja $S^2 = \frac{1}{n-1} \sum_{i = 1}^n (X_i - \bar X)^2$.

Voimme testata seuraavaa nollahypoteesia $H_0$ vaihtoehtoista hypoteesia $H_1$ vastaan \[\begin{align*} &H_0: \mu = \mu_0, \\ &H_1: \mu \neq \mu_0, \end{align*}\] (kaksisuuntaisella) yhden otoksen t-testillä. Tällöin satunnaisotoksesta laskettu testisuure noudattaa $t$-jakaumaa vapausasteella $n - 1$, \[\begin{equation} \tag{4.1} T = \frac{\bar X - \mu_0}{S/\sqrt{n}} \sim t(n-1). \end{equation}\] Jos havaitusta otoksesta laskettu testisuure $t$ on tarpeeksi poikkeava suhteessa valittuun merkitsevyystasoon, niin hylkäämme nollahypoteesin $H_0$. Lämmittelynä tämän luvun aiheisiin käymme läpi esitellyn t-testin suorittamisen esimerkin kautta. Kyseinen esimerkki on itseasiassa suoraan kurssilta Tilastotieteen perusteet.

Esimerkki 4.1 (Yhden otoksen t-testi odotusarvolle) Vilin makeistehtaat tuottavat suklaalevyjä. Vili omistaa yhteensä $4$ tehdasta. Johtoporras on asettanut laatukriteeriksi, että suklaalevyjen keskimääräisen painon pitäisi olla $\mu = 100\mathrm{g}$ joka tehtaassa. Nyt laaduntarkastajat nappaavat jokaisen tehtaan tuotantolinjalta $20$ suklaalevyn otoksen. Ensimmäistä tehdasta vastaava otos näkyy alla.

choco[, 1]

## # A tibble: 20 × 1
##    factory1
##       <dbl>
##  1    103. 
##  2     98.3
##  3    111. 
##  4    101. 
##  5    107. 
##  6    102. 
##  7    107. 
##  8    101. 
##  9     97.8
## 10    107. 
## 11     93.4
## 12    103. 
## 13     91.8
## 14     88.9
## 15     97.4
## 16     98.8
## 17    108. 
## 18    101. 
## 19     99.3
## 20    104.

Voimme laskea ensimmäisen tehtaan otosta vastaavan keskipainon ja otoskeskihajonnan.

x_bar <- mean(choco$factory1)
s <- sqrt(var(choco$factory1))
x_bar

## [1] 101.0809

## [1] 5.629763

Keskipainoksi saatiin $\bar x \approx 101.08\mathrm{g}$ ja otoskeskihajonnaksi $s \approx 5.63\mathrm{g}$.

Oletetaan, että suklaalevyjen painot ovat normaalijakautuneita.

Onko ensimmäisen tehtaan tuotantolinjassa jotain vikaa vai selittyykö havaittu poikkeama (tavoite $\mu = 100$ vs havaittu $\bar x \approx 101.08$) satunnaisvaihtelulla?

Analysoidaan kysymystä yhden otoksen t-testin avulla käyttäen merkitsevyystasoa $\alpha = 0.05$. R:llä voimme laskea testisuureen arvon $t$ poimittujen $20$ suklaalevyn pohjalta.

n <- 20
mu_0 <- 100

t <- sqrt(n) * (x_bar - mu_0) / s
t

## [1] 0.8586106

Koska $t$-jakauma on symmetrinen ja tässä tapauksessa $t\approx 0.86 > 0$, niin $p$-arvon kaavaksi saadaan \[\begin{equation*} \begin{split} \mathbb{P}\left(|T|\geq |t|\right) &= \mathbb{P}\left(|T|\geq t\right) = 2\mathbb{P}\left(T\geq t\right) \\ &= 2\mathbb{P}\left(T\leq -t\right) = 2F_{n-1}(-t), \end{split} \end{equation*}\] jossa $F_{n-1}$ on jakauman $t(n-1)$ kertymäfunktio. Nyt p-arvo voidaan laskea R:n avulla.

p_value <- 2 * pt(-t, df = n - 1)
p_value

## [1] 0.401258

Koska $p$-arvo ($\approx 0.40$) on suurempi kuin asetettu merkitsevyystaso $\alpha = 0.05$, niin emme hylkää nollahypoteesiä. Testin tulos ei siis vielä riitä vakuuttamaan laadunvalvontaa siitä, että ensimmäisen tehtaan tuotantolinjassa olisi jotain vikaa.

Testiä ei tarvitse suorittaa R:ssä manuaalisesti. Sama tulos saadaan käyttämällä t.test() funktiota.

t.test(choco$factory1, mu = 100)

## 
##  One Sample t-test
## 
## data:  choco$factory1
## t = 0.85861, df = 19, p-value = 0.4013
## alternative hypothesis: true mean is not equal to 100
## 95 percent confidence interval:
##   98.44605 103.71568
## sample estimates:
## mean of x 
##  101.0809

Tässä luvussa emme kuitenkaan ole niinkään kiinnostuneita testaamaan, onko odotusarvon $\mu$ arvo jokin tietty $\mu_0$, vaan tavoiteenamme on vertailla useampaa ryhmää keskenään. Tämä osio huipentuukin luvussa 4.3 esiteltävään yksisuuntaiseen varianssianalyysin (eng: one-way analysis of variance, ANOVA), joka on tilastollinen testi useamman ryhmän odotusarvon vertailuun. Ennen ANOVA:a, luvussa 4.2 käymme läpi kahden ryhmän odotusarvojen vertailua kahden otoksen t-testillä.

Hieman yllättäen ANOVA:n voi ymmärtää intuitiivisesti lineaarisen regression kautta. Tämän vuoksi luvussa 4.1 käsittelemme luokkamuuttujia lineaarisessa regressiomallissa, sillä tilastolliset testit voidaan usein muotoilla lineaarisen regression kautta nimenomaan luokkamuuttujien avulla. Myös yhden otoksen t-testi voidaan muotoilla lineaarisen regression kautta, vaikka tämä ei välttämättä ole kovin hyödyllistä.

Esimerkki 4.2 (Yhden otoksen t-testi ja lineaarinen regressio) Jatketaan lämmittelyesimerkkiä 4.1. Määritellään lineaarinen regressiomalli \[\begin{equation} \tag{4.2} \texttt{factory1}_i - \mu_0 = \beta_0 + \varepsilon_i, \quad \varepsilon_i\sim N(0, \sigma^2), \end{equation}\] jossa $\mu_0 = 100$. Huomaa, että yllä olevassa mallissa on vain vakiotermi $\beta_0$ eikä yhtään selittäjää! Nollahypoteesin $\mu = \mu_0$ alla suklaalevyjen painot noudattavat normaalijakaumaa $N(\mu_0, \sigma^2)$. Toisin sanoen nollahypoteesin alla¹⁵ mallissa (4.2) asetetaan $\beta_0 = 0$. Tämän vuoksi yhden otoksen t-testi odotusarvolle palautuu mallin (4.2) parametrin $\beta_0$ tilastollisen merkitsevyyden testaamiseen luvun 3.5 työkaluilla.

Sovitetaan kaavan (4.2) lineaarinen malli. Lineaarinen regressiomalli, jossa on ainoastaan vakiotermi saadaan käyttämällä syntaksia y ~ 1 funktiossa lm().

factory_lm <- lm(choco$factory1 - mu_0 ~ 1)

Komennon summary(factory_lm) tulosteesta voimme lukea vakiotermiä vastaavan merkitsevyystestin p-arvon.

summary(factory_lm)

## 
## Call:
## lm(formula = choco$factory1 - mu_0 ~ 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.2170  -2.8902   0.1911   3.7369   9.5701 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)    1.081      1.259   0.859    0.401
## 
## Residual standard error: 5.63 on 19 degrees of freedom

Huomataan, että p-arvo ($\approx 0.40$) on sama, joka saatiin esimerkissä 4.1!

Yksinkertaisuuden vuoksi tämän luvun testeissä oletamme datan noudattavan normaalijakaumaa. Huomautamme kuitenkin, että tämän luvun testit ovat usein käyttökelpoisia myös tapauksissa, joissa data ei noudata normaalijakaumaa kunhan otoskoko on tarpeeksi suuri. Tämä voidaan perustella keskeisen raja-arvolauseen avulla. Esimerkiksi kaavan (4.1) testisuure noudattaa likimain standardinormaalijakaumaa $N(0,1)$ suurella otoskoolla. Vastaus kysymykseen, kuinka suuri otoskoon täytyy olla, jotta voimme perustella keskeisen raja-arvolauseen käytön onkin sitten monimutkaisempi…

4.1 Luokkamuuttujat lineaarisessa regressiossa

Tarkastellaan tilannetta, jossa haluamme selittää erään yrityksen työntekijöiden palkkaa (salary, yksikkönä tuhat euroa) muuttujilla työkokemus (experience, yksikkönä vuosi) ja työllisyystyyppi (type, kokoaikainen tai osa-aikainen). Muuttuja experience on intervalliasteikkoinen muuttuja. Eli voimme mitata eroja työkokemuksessa $|\texttt{experience}_i - \texttt{experience}_j|$ ja jokaista muuttujan experience arvoa vastaa yksikäsitteinen luku (esim. $2$ vuotta tai $15.3$ vuotta). Työllisyystyyppi on kuitenkin luokkamuuttuja (eng: categorical variable), ja luokkamuuttujan arvojen koodaus lukuarvoiksi voidaan tehdä monella tavalla. Voimme esimerkiksi määritellä \[\begin{align*} &\text{kokoaikainen} = 1, \\ &\text{osa-aikainen} = 0 \end{align*}\] tai \[\begin{align*} &\text{kokoaikainen} = 1, \\ &\text{osa-aikainen} = -1. \end{align*}\]

Lineaarisessa regressiossa luokkamuuttuja halutaan ottaa tietyllä tavalla malliin mukaan. Tähän käytämme osoitinmuuttujia (eng: dummy variable, indicator variable), jotka ovat 0–1-arvoisia satunnaismuuttujia. Tässä tapauksessa luokkamuuttujalla type on vain kaksi luokkaa (kokoaikainen/osa-aikainen), joten meidän tarvitsee sisällyttää vain yksi osoitinmuuttuja lineaariseen malliin. Alla olevassa esimerkissä havainnollistamme yllä kuvaillun lineaarisen mallin rakennusta ja tulkitsemista.

Esimerkki 4.3 (Luokkamuuttujan tulkinta lineaarisessa regressiossa) Osoitinmuuttujan \[\begin{equation*} \texttt{type}_i = \begin{cases} 0,& \text{jos työntekijä $i$ on osa-aikainen}, \\ 1,& \text{jos työntekijä $i$ on kokoaikainen}, \end{cases} \end{equation*}\] avulla määrittelemme lineaarisen regressiomallin, \[\begin{equation*} \texttt{salary}_i = \beta_0 + \beta_1\texttt{experience}_i + \beta_2\texttt{type}_i + \varepsilon_i. \end{equation*}\] Kuvitellaan, että pienimmän neliösumman menetelmällä saamme sovitteen \[\begin{equation*} \widehat{\texttt{salary}}_i = \hat\beta_0 + \hat\beta_1\texttt{experience}_i + \hat\beta_2\texttt{type}_i, \end{equation*}\] jossa $\hat\beta_0 = 2$, $\hat\beta_1 = 0.1$ ja $\hat\beta_2 = 1$. Toisaalta sijoittamalla osoitinmuuttujan arvot ($0$ tai $1$) voimme kirjoittaa sovitteen muodossa \[\begin{equation*} \widehat{\texttt{salary}}_i = \begin{cases} \hat\beta_0 + \hat\beta_1\texttt{experience}_i,& \text{jos työntekijä $i$ on osa-aikainen}, \\ (\hat\beta_0 + \hat\beta_2) + \hat\beta_1\texttt{experience}_i, & \text{jos työntekijä $i$ on kokoaikainen}. \end{cases} \end{equation*}\] Tulkinta luokkamuuttujaa vastaavalle osoitinmuuttujalle on siis seuraava.

Kokoaikainen työntekijä tienaa keskimäärin $\hat\beta_2 = 1$ (eli tuhat euroa) enemmän kuin saman kokemuksen omaava osa-aikainen työntekijä.

Lineaariseen regressiomalliin lisättävässä luokkamuuttujassa voi myös olla enemmän kuin kaksi luokkaa – kun luokkia on $k$ kappaletta, niin osoitinmuuttujia täytyy olla $k - 1$ kappaletta.

Esimerkki 4.4 (Useampi luokka lineaarisessa regressiossa) Jatketaan esimerkkiä 4.3. Työllisyystyypin type sijasta sisällytetään luokkamuuttuja, joka kuvaa työntekijän toimialaa (rahoitus, markkinointi tai kirjanpito). Tässä tapauksessa voimme sisällyttää kaksi osoitinmuuttujaa \[\begin{equation*} \texttt{r}_i = \begin{cases} 1,& \text{jos työntekijän $i$ toimiala on rahoitus}, \\ 0,& \text{jos työntekijä $i$ toimiala on joku muu}, \end{cases} \end{equation*}\] ja \[\begin{equation*} \texttt{m}_i = \begin{cases} 1,& \text{jos työntekijän $i$ toimiala on markkinointi}, \\ 0,& \text{jos työntekijä $i$ toimiala on joku muu}, \end{cases} \end{equation*}\] jolloin saamme lineaarisen regressiomallin \[\begin{equation*} \texttt{salary}_i = \beta_0 + \beta_1\texttt{experience}_i + \beta_2\texttt{r}_i + \beta_3\texttt{m}_i + \varepsilon_i. \end{equation*}\] Kuvitellaan, että pienimmän neliösumman menetelmällä saamme sovitteen \[\begin{equation} \tag{4.3} \widehat{\texttt{salary}}_i = \hat\beta_0 + \hat\beta_1\texttt{experience}_i + \hat\beta_2\texttt{r}_i + \hat\beta_3\texttt{m}_i. \end{equation}\] Luokkamuuttujia vastaavien osoitinmuuttujien tulkinnat ovat seuraavat.

Estimaatti $\hat\beta_2$ kuvaa sitä, kuinka paljon vähemmän/enemmän rahoituksessa työskentelevä tienaa keskimäärin kuin kirjanpidossa, jos työntekijöillä on sama määrä kokemusta.
Estimaatti $\hat\beta_3$ kuvaa sitä, kuinka paljon vähemmän/enemmän markkinoinnissa työskentelevä tienaa keskimäärin kuin kirjanpidossa, jos työntekijöillä on sama määrä kokemusta.

Huomaa, että olisimme voineet valita osoitinmuuttujat eri tavalla, joka muuttaisi estimaattien $\hat\beta_0$, $\hat\beta_2$ ja $\hat\beta_3$ tulkitsemista. Nimittäin sovitteen (4.3) vakiotermin estimaatti vastaa kirjanpidossa työskentelevän ihmisen keskimääräistä aloituspalkkaa. Toisaalta, jos valitsisimme osoitinmuuttujat vastaamaan luokkia kirjanpito ja markkinointi, niin vakiotermin estimaatti vastaisi rahoituksessa työskentelevän ihmisen keskimääräistä aloituspalkkaa.

Luokkamuuttujia voi siis sisällyttää melko vapaasti lineaariseen malliin. Huomioi kuitenkin, että luokkamuuttujaa vastaavien luokkien täytyy olla

toisensa posisulkevia ja
jokaisen havainnon täytyy kuulua johonkin luokkaan.

Toisin sanoen luokkamuuttujat täytyy olla muodostettu niin, että jokainen havainto kuuluu tasan yhteen luokkaan. Lisäksi harvinaiset luokat tekevät mallin estimaateista epävakaita. Joskus ennen mallin sovittamista luokkia täytyykin manuaalisesti muokata tai yhdistellä.

Kun R ohjelmointikielellä sovitetaan lineaarinen regressiomalli, joka sisältää luokkamuuttujia, niin vastaavia osoitinmuuttujia ei tarvitse tehdä itse manuaalisesti. Luokkamuuttujat täytyy vain antaa tietotyyppinä factor, joka on erityisesti luokkamuuttujille räätälöity tietotyyppi R:ssä.

Esimerkki 4.5 (Luokkamuuttujat lineaarisessa regressiossa ja R) Käytetään R:ssä olevaa valmista havaintoaineistoa mtcars. Kyseinen aineisto sisältää mittauksia autojen eri ominaisuuksista kuten polttoaineen kulutuksesta mpg (mailit per gallona), auton painosta wt (tuhansissa nauloissa) ja sylintereiden määrästä cyl.

mtcars <- tibble::as_tibble(mtcars)
mtcars

## # A tibble: 32 × 11
##      mpg   cyl  disp    hp  drat    wt  qsec    vs    am  gear  carb
##    <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1  21       6  160    110  3.9   2.62  16.5     0     1     4     4
##  2  21       6  160    110  3.9   2.88  17.0     0     1     4     4
##  3  22.8     4  108     93  3.85  2.32  18.6     1     1     4     1
##  4  21.4     6  258    110  3.08  3.22  19.4     1     0     3     1
##  5  18.7     8  360    175  3.15  3.44  17.0     0     0     3     2
##  6  18.1     6  225    105  2.76  3.46  20.2     1     0     3     1
##  7  14.3     8  360    245  3.21  3.57  15.8     0     0     3     4
##  8  24.4     4  147.    62  3.69  3.19  20       1     0     4     2
##  9  22.8     4  141.    95  3.92  3.15  22.9     1     0     4     2
## 10  19.2     6  168.   123  3.92  3.44  18.3     1     0     4     4
## # ℹ 22 more rows

Havaintoaineistossa sarakkeen cyl tietotyyppi on desimaaliluku, vaikka kyseinen muuttuja saa vain arvoja $4$, $6$ ja $8$. Eli cyl on luokkamuuttuja.

unique(mtcars$cyl)

## [1] 6 4 8

Selitämme polttoaineen kulutusta painolla ja sylintereiden määrällä. Muuttamalla muuttujan cyl tietotyyppi funktiolla factor() ennen lineaarisen mallin sovitusta, R ymmärtää, että kyseessä on luokkamuuttuja.

cars_lm <- lm(mpg ~ wt + factor(cyl), data = mtcars)
summary(cars_lm)

## 
## Call:
## lm(formula = mpg ~ wt + factor(cyl), data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5890 -1.2357 -0.5159  1.3845  5.7915 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   33.9908     1.8878  18.006  < 2e-16 ***
## wt            -3.2056     0.7539  -4.252 0.000213 ***
## factor(cyl)6  -4.2556     1.3861  -3.070 0.004718 ** 
## factor(cyl)8  -6.0709     1.6523  -3.674 0.000999 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.557 on 28 degrees of freedom
## Multiple R-squared:  0.8374, Adjusted R-squared:   0.82 
## F-statistic: 48.08 on 3 and 28 DF,  p-value: 3.594e-11

Komennon summary(cars_lm) tulosteesta näemme, että sovitettu malli on muotoa \[\begin{equation*} \widehat{\texttt{mpg}}_i = \hat\beta_0 + \hat\beta_1\texttt{wt}_i + \hat\beta_2\texttt{c_six}_i + \hat\beta_3\texttt{c_eight}_i, \end{equation*}\] jossa \[\begin{equation*} \texttt{c_six}_i = \begin{cases} 1,& \text{jos autossa on 6 sylinteriä}, \\ 0,& \text{muuten}, \end{cases} \end{equation*}\] ja \[\begin{equation*} \texttt{c_eight}_i = \begin{cases} 1,& \text{jos autossa on 8 sylinteriä}, \\ 0,& \text{muuten}. \end{cases} \end{equation*}\]

Kun auton paino pidetään vakiona, niin sylintereiden määrän kasvaessa suhde $\text{miles}/\text{gallona}$ pienenee eli polttoaineen kulutus kasvaa.

4.2 Kahden otoksen t-testi

Olkoon $X_1$ ja $X_2$ satunnaismuuttujia, jotka noudattavat normaalijakaumaa samalla varianssilla $\sigma^2$ mutta mahdollisesti eri odotusarvoilla $\mu_1$ ja $\mu_2$. Toisin sanoen $X_1\sim N(\mu_1, \sigma^2)$ ja $X_2\sim N(\mu_2, \sigma^2)$, ja haluamme testata seuraavaa nollahypoteesia $H_0$ vaihtoehtoista hypoteesia $H_1$ vastaan \[\begin{align*} &H_0: \mu_1 = \mu_2, \\ &H_1: \mu_1 \neq \mu_2. \end{align*}\] Olkoon $\overrightarrow{X_1} = \{X_{11}, \ldots, X_{n_11}\}$ satunnaisotos satunnaismuuttujasta $X_1$ ja $\overrightarrow{X_2} = \{X_{12}, \ldots, X_{n_22}\}$ satunnaisotos satunnaismuuttujasta $X_2$. Meidän tapauksessa oletamme satunnaisotoksista, että

otoskoot $n_1$ ja $n_2$ voivat olla erisuuria ja
yhdistetyn otoksen $\overrightarrow{X_1} \cup \overrightarrow{X_1} = \{X_{11}, \ldots, X_{n_11}, X_{12}, \ldots, X_{n_22}\}$ kaikki havainnot ovat riippumattomia.

Annetussa asetelmassa voimme soveltaa kahden otoksen t-testiä. Katso huomio 4.1 liittyen muuanlaisiin t-testeihin, joissa verrataan kahden ryhmän odotusarvoja. Merkitään satunnaisotoksia vastaavia keskiarvoja notaatioilla \[\begin{equation*} \bar X_1 = \frac{1}{n_1}\sum_{i = 1}^{n_1} X_{i1} \quad\text{ja}\quad \bar X_2 = \frac{1}{n_2}\sum_{i = 1}^{n_2} X_{i2}. \end{equation*}\] Vastaavasti merkitään satunnaisotoksia vastaavia variansseja notaatioilla \[\begin{equation*} S_1^2 = \frac{1}{n_1 - 1}\sum_{i = 1}^{n_1} \left(X_{i1} - \bar X_1\right)^2 \quad\text{ja}\quad S_2^2 = \frac{1}{n_2 - 1}\sum_{i = 1}^{n_2} \left(X_{i2} - \bar X_2\right)^2. \end{equation*}\] Koska oletamme, että satunnaismuutujilla $X_1$ ja $X_2$ on sama varianssi, niin voimme laskea yhdistetyn (eng: pooled) estimaatin varianssille kaavalla \[\begin{equation*} S_p^2 = \frac{(n_1 - 1)S_1^2 + (n_2 - 1)S_2^2}{n_1 + n_2 - 2}. \end{equation*}\] Tällöin kun antamamme oletukset pätevät, niin nollahypoteesin alla testisuure $T$ noudattaa $t$-jakaumaa $n_1 + n_2 - 2$ vapausasteella, \[\begin{equation} \tag{4.4} T = \frac{\bar X_1 - \bar X_2}{S_p\sqrt{\frac{1}{n_1} + \frac{1}{n_2}}} \sim t(n_1 + n_2 - 2). \end{equation}\]

Alla on käytännön esimerkki siitä, miten kahden otoksen t-testi suoritetaan R:llä. Suoritamme testin manuaalisesti ja käyttäen R:n funktiota t.test().

Esimerkki 4.6 (Kahden otoksen t-testi) Esimerkissä 4.1 totesimme jo, että ensimmäisessä tehtaassa laatukriteeri $\mu_1 = 100\mathrm{g}$ liittyen suklaalevyjen keskipainoon täyttyy tarvittavalla tarkkuudella. Nyt johto haluaa vastauksen seuraavaan kysymykseen.

Painolaatukriteerin suhteen, onko toinen tehdas yhtä luotettava kuin ensimmäinen?

Analysoidaan kysymystä kahden otoksen t-testin avulla. Käytetään taas merkitsevyystasoa $\alpha = 0.05$.

Alla näkyy ensimmäistä ja toista tehdasta vastaavat otokset suklaalevyistä.

choco[, 1:2]

## # A tibble: 20 × 2
##    factory1 factory2
##       <dbl>    <dbl>
##  1    103.     102. 
##  2     98.3    111. 
##  3    111.     102. 
##  4    101.      93.7
##  5    107.     104. 
##  6    102.     101. 
##  7    107.      91.2
##  8    101.      87.8
##  9     97.8     93.7
## 10    107.      95.6
## 11     93.4    100.0
## 12    103.     105. 
## 13     91.8     97.7
## 14     88.9     96.8
## 15     97.4     98.4
## 16     98.8     96.4
## 17    108.      98.0
## 18    101.      98.6
## 19     99.3     94.8
## 20    104.      94.9

Ensin laskemme keskiarvot $\bar X_i$ ja otosvarianssit $S_i^2$ tehtaille $j\in\{1,2\}$ sekä yhdistetyn estimaatin yhteiselle varianssille $S^2_p$.

x_bar1 <- mean(choco$factory1)
x_bar2 <- mean(choco$factory2)

var1 <- var(choco$factory1)
var2 <- var(choco$factory2)

var_pooled <- (n - 1) * (var1 + var2) / (2 * n - 2)

Saamme testisuureeksi seuraavanlaisen arvon.

t <- (x_bar1 - x_bar2) / (sqrt(var_pooled * 2/n))
t

## [1] 1.740809

Nyt voimme laskea vastaavan p-arvon.

2 * pt(-t,  df = 2 * n - 2)

## [1] 0.08980959

Koska $p$-arvo ($\approx 0.090$) on suurempi kuin asetettu merkitsevyystaso $\alpha = 0.05$, niin emme hylkää nollahypoteesiä $\mu_1 = \mu_2$. Testin tulos ei siis vielä riitä vakuuttamaan laadunvalvontaa siitä, että toisen tehtaan tuotantolinja olisi epäluotettavampi kuin ensimmäisen. Samaan tulokseen päädytään käyttämällä funktiota t.test().

t.test(choco$factory1, choco$factory2, var.equal = TRUE, paired = FALSE)

## 
##  Two Sample t-test
## 
## data:  choco$factory1 and choco$factory2
## t = 1.7408, df = 38, p-value = 0.08981
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4851756  6.4417612
## sample estimates:
## mean of x mean of y 
## 101.08086  98.10257

Kuten yhden otoksen t-testi, niin myös kahden otoksen t-testi voidaan muotoilla lineaarisen regression avulla. Määritellään lineaarinen regressiomalli \[\begin{equation} \tag{4.5} X_{ij} = \beta_0 + \beta_1 D_{i2} + \varepsilon_i, \quad \varepsilon_i\sim N(0,\sigma^2), \end{equation}\] jossa $D_{i2}$ on osoitinmuuttuja, \[\begin{equation*} D_{ij} = \begin{cases} 0,& \text{jos havainto $X_{ij}$ kuuluu luokkaan $j = 1$}, \\ 1,& \text{jos havainto $X_{ij}$ kuuluu luokkaan $j = 2$}. \end{cases} \end{equation*}\] Nollahypoteesin $\mu_1 = \mu_2$ alla kerroin $\beta_1$ ei ole tilastollisesti merkitsevä. Eli mallissa (4.5)

vakiotermi $\beta_0$ vastaa ensimmäisen ryhmän odotusarvoa $\mu_1$ ja
$\beta_1$ vastaa ryhmien odotusarvojen välistä erotusta $\mu_2 - \mu_1$. Nollahypoteesin alla $\beta_1 = \mu_2 - \mu_1 = 0$.

Esimerkki 4.7 (Kahden otoksen t-testi ja lineaarinen regressio) Suoritetaan esimerkin 4.6 analyysi lineaarisen regression kautta. Ensin dataan täytyy sisällyttää luokkamuuttuja, joka kertoo sen, kummasta tehtaasta suklaalevy tulee. Tätä tarkoitusta varten on useampikin mahdollinen R funktio. Tässä käytämme tidyr paketin funktiota gather(). Fuktion gather() syntaksin opettelu ei ole oleellista kurssin kannalta. On kuitenkin hyvä tietää funktion olemassaolosta, jos tulevaisuudessa joudut siivoamaan ja muokkaamaan dataa.

library(tidyr)
library(dplyr)
choco12 <- choco %>%
  select(factory1, factory2) %>%
  gather("factory", "weight") %>%
  mutate(factory = factor(factory))

choco12

## # A tibble: 40 × 2
##    factory  weight
##    <fct>     <dbl>
##  1 factory1  103. 
##  2 factory1   98.3
##  3 factory1  111. 
##  4 factory1  101. 
##  5 factory1  107. 
##  6 factory1  102. 
##  7 factory1  107. 
##  8 factory1  101. 
##  9 factory1   97.8
## 10 factory1  107. 
## # ℹ 30 more rows

Nyt sovitamme kaavan (4.5) mukaisen lineaarisen regressiomallin. Huomaa, että yllä asetimme luokkamuuttujan tietotyypiksi factor jo valmiiksi.

factory12_lm <- lm(weight ~ factory, data = choco12)

Komennon summary(factory12_lm) tulosteesta voimme lukea kulmakerrointa $\beta_1$ vastaavan merkitsevyystestin p-arvon. Tulosteesta nähdään myös, että tässä tapauksessa R muodosti osoitinmuuttujan tyyliin \[\begin{equation*} D_{ij} = \begin{cases} 0,& \text{jos suklaalevy $X_{ij}$ tulee tehtaasta $1$}, \\ 1,& \text{jos suklaalevy $X_{ij}$ tulee tehtaasta $2$}. \end{cases} \end{equation*}\]

summary(factory12_lm)

## 
## Call:
## lm(formula = weight ~ factory, data = choco12)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.2170  -3.2066  -0.0284   3.2175  12.9144 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      101.081      1.210  83.554   <2e-16 ***
## factoryfactory2   -2.978      1.711  -1.741   0.0898 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.41 on 38 degrees of freedom
## Multiple R-squared:  0.07386,    Adjusted R-squared:  0.04949 
## F-statistic:  3.03 on 1 and 38 DF,  p-value: 0.08981

Huomataan, että p-arvo ($\approx 0.090$) on sama, joka saatiin esimerkissä 4.1 (lue p-arvo riviltä factoryfactory2)!

Huomautus 4.1 (Kahden ryhmän odotusarvojen vertailusta) Tämän luvun versiossa kahden otoksen t-testistä oletimme, että muuttujien $X_1$ ja $X_2$ varianssit ovat yhtä suuret. Näin ei tarvitse kuitenkaan olla ja kahden otoksen t-testin voi myös suorittaa olettaen, että varianssit ovat erisuuret. Tällöin tosin testisuure näyttää erilaiselta kuin kaavassa (4.4).

Kahden otoksen t-testi ei ole soveltuva, jos havaintoparin sisällä $(X_{i1}, X_{i2})^T$ olevat havainnot ovat riippuvaisia toisistaan. Tapaukseen, jossa otoskoot ovat yhtäsuuret $n_1 = n_2$, havaintoparien sisällä on riippuvuutta mutta parit $(X_{11}, X_{12})^T, \ldots, (X_{n_11}, X_{n_12})^T$ ovat riippumattomia, sopii parittainen t-testi (eng: paired t-test). Parittaista t-testiä ei käydä läpi tässä luvussa.

4.3 Yksisuuntainen varianssianalyysi

Yksisuuntainen varianssianalyysi (eng: one-way analysis of variance, ANOVA) on yleistys kahden otoksen t-testistä useamman ryhmän odotusarvon vertailuun. Olkoon $X_1, X_2, \ldots, X_d$ satunnaismuuttujia, jotka noudattavat normaalijakaumaa yhteisellä varianssilla $\sigma^2$ mutta mahdollisesti eri odotusarvoilla $\mu_1, \mu_2, \ldots, \mu_d$. Toisin sanoen $X_j\sim N(\mu_j, \sigma)$ ja haluamme testata seuraavaa nollahypoteesia $H_0$ vaihtoehtoista hypoteesia $H_1$ vastaan \[\begin{align*} &H_0: \mu_1 = \mu_2 = \cdots = \mu_d, \\ &H_1: \mu_j \neq \mu_k \quad\text{jollekin}\quad j,k\in\{1, 2, \ldots, d\}. \end{align*}\] Olkoon $\overrightarrow{X_j} = \{X_{1j}, \ldots, X_{n_jj}\}$ satunnaisotos satunnaismuuttujasta $X_j$. Oletamme satunnaisotoksista, että

otoskoot $n_j$, $j\in\{1, 2, \ldots, d\}$ voivat olla erisuuria ja
yhdistetyn otoksen $\overrightarrow{X_1} \cup \overrightarrow{X_2} \cup \cdots \cup \overrightarrow{X_d}$ kaikki havainnot ovat riippumattomia.

Annetussa asetelmassa voimme soveltaa ANOVA:a. Yhdistetyn otoksen otoskoko on $n = \sum_{j = 1}^d n_j$ ja yhdistetyn otoksen keskiarvo on \[\begin{equation*} \bar X = \frac{1}{n} \sum_{j = 1}^d \sum_{i = 1}^{n_j} X_{ij}. \end{equation*}\] Vastaavasti merkitään ryhmäkeskiarvoja notaatiolla \[\begin{equation*} \bar X_j = \frac{1}{n_j}\sum_{i = 1}^{n_j} X_{ij}. \end{equation*}\] Merkitään myös ryhmien välistä neliösummaa (RVN) ja ryhmien sisäistä neliösummaa (RSN) notaatioilla \[\begin{equation*} \mathrm{RVN} = \sum_{j = 1}^{d} n_j\left(\bar X_j - \bar X\right)^2 \quad\text{ja}\quad \mathrm{RSN} = \sum_{j = 1}^{d}\sum_{i = 1}^{n_j} \left(X_{ij} - \bar X_j\right)^2. \end{equation*}\] Tällöin kun antamamme oletukset pätevät, niin nollahypoteesin alla testisuure $F$ noudattaa $F$-jakaumaa vapausasteilla $d - 1$ ja $n - d$, \[\begin{equation*} F = \frac{\mathrm{RVN}/(d - 1)}{\mathrm{RSN}/(n - d)}\sim F(d - 1, n - d). \end{equation*}\] ANOVA:n ideana on siis vertailla ryhmien välistä vaihtelua ryhmien sisäiseen vaihteluun, \[\begin{equation*} \frac{\text{Ryhmien välinen vaihtelu}}{\text{Ryhmien sisäinen vaihtelu}}. \end{equation*}\] Jos ryhmien välinen vaihtelu on suuri suhteessa ryhmien sisäiseen vaihteluun, niin nollahypoteesi $\mu_1 = \mu_2 = \cdots = \mu_d$ hylätään. Itseasiassa myös testin nimi “varianssianalyysi” liittyy siihen, että odotusarvojen yhtäsuuruuksia testataan vertaamalla erilaisia vaihtelua kuvaavia suureita keskenään (katso huomautus 4.3).

Alla on käytännön esimerkki siitä, miten ANOVA suoritetaan R:llä. Suoritamme testin manuaalisesti ja käyttäen R:n funktiota aov().

Esimerkki 4.8 (ANOVA) Esimerkeissä 4.1 ja 4.6 tehtyjen tilastollisten testien perusteella suklaalevyjen keskimääräiseen painoon $100\mathrm{g}$ liittyvä laatukriteeri täyttyy tehtaissa 1 ja 2 tarvittavalla tarkkuudella. Nyt johto antaa seuraavan tehtävänannon.

Vilin tehtaissa on tärkeää, että laatu pysyy tasaisena. Onko joku tehdas, jossa suklaalevyjen keskimääräinen paino on poikkeava verrattuna muihin tehtaisiin?

Analysoidaan johdon asettamaa kysymystä ANOVA:n avulla. Käytetään merkitsevyystasoa $\alpha = 0.05$.

Alla näkyy kaikkia neljää tehtaita vastaavat otokset suklaalevyistä.

choco

## # A tibble: 20 × 4
##    factory1 factory2 factory3 factory4
##       <dbl>    <dbl>    <dbl>    <dbl>
##  1    103.     102.     101.      89.5
##  2     98.3    111.     110.      91.8
##  3    111.     102.     103.      82.8
##  4    101.      93.7    110.      78.0
##  5    107.     104.      97.2     98.3
##  6    102.     101.      95.6     89.5
##  7    107.      91.2     93.9     92.4
##  8    101.      87.8     88.4     83.5
##  9     97.8     93.7    102.      99.6
## 10    107.      95.6    108.      96.0
## 11     93.4    100.0    102.      92.3
## 12    103.     105.     105.      86.6
## 13     91.8     97.7     97.9     82.7
## 14     88.9     96.8     96.7     81.8
## 15     97.4     98.4    101.      81.5
## 16     98.8     96.4     98.2     88.5
## 17    108.      98.0     95.0     82.6
## 18    101.      98.6    105.      91.3
## 19     99.3     94.8    109.      90.4
## 20    104.      94.9    104.      95.1

Ensin laskemme tehdaskohtaiset keskiarvot $\bar X_j$, yhteisen keskiarvon $\bar X$ sekä tunnusluvut $\mathrm{RVN}$ ja $\mathrm{RSN}$.

n <- nrow(choco) * ncol(choco)
n_j <- nrow(choco)

mean_group <- colMeans(choco)
mean_total <- mean(mean_group)

rvn <- sum(n_j * (mean_group - mean_total)^2)
rsn <- sum(as.matrix(sweep(choco, 2, mean_group, "-"))^2)

rvn

## [1] 2058.938

rsn

## [1] 2454.718

Saamme testisuureeksi arvon $f\approx 21.25$.

d <- ncol(choco)
f <- ((n - d) / (d - 1)) * (rvn / rsn)
f

## [1] 21.24879

Nyt voimme laskea testisuuretta vastaavan p-arvon, joka on tässä tapauksessa lähellä nollaa.

1 - pf(f, d - 1, n - d)

## [1] 4.271928e-10

Koska p-arvo on pienempi kuin asetettu merkitsevyystaso $\alpha = 0.05$, niin hylkäämme nollahypoteesin. Testin tulos antaa tarpeeksi todisteita sen puolesta, että ainakin yhdessä tehtaassa tuotettujen suklaalevyjen keskipaino eroaa merkittävästi muissa tehtaissa tuotettujen suklaalevyjen keskipainosta, jonka pitäisi olla $100\mathrm{g}$.

Sama tulos saadaan käyttämällä valmista R funktiota aov(). Jotta funktiota voidaan käyttää, niin dataa pitää muokata esimerkin 4.7 tyyliin. Muutamme datan siis sellaiseen muotoon, jossa ensimmäinen sarake kertoo suklaalevyn tuotantosijainnin ja toinen sarake levyä vastaavan painon.

choco_mod <- choco %>%
  gather("factory", "weight") %>%
  mutate(factory = factor(factory))

choco_mod

## # A tibble: 80 × 2
##    factory  weight
##    <fct>     <dbl>
##  1 factory1  103. 
##  2 factory1   98.3
##  3 factory1  111. 
##  4 factory1  101. 
##  5 factory1  107. 
##  6 factory1  102. 
##  7 factory1  107. 
##  8 factory1  101. 
##  9 factory1   97.8
## 10 factory1  107. 
## # ℹ 70 more rows

Nyt muokatulle datalla choco_mod voimme suorittaa ANOVA:n. Funktion aov() syntaksi on samantyyppinen kuin funktion lm(). Eli havainnot weight ja luokkamuuttuja factory erotetaan tilden ~ avulla. Testin tuloksen näemme käyttämällä summary funktiota.

choco_anova <- aov(weight ~ factory, data = choco_mod)
summary(choco_anova)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## factory      3   2059   686.3   21.25 4.27e-10 ***
## Residuals   76   2455    32.3                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sarake Sum Sq antaa tunnuslukujen $\mathrm{RVN}$ ja $\mathrm{RSN}$ arvot. Keskineliösummien $\mathrm{RVN}/(d - 1)$ ja $\mathrm{RSN}/(n - d)$ arvot näkyvät sarakkeessa Mean Sq. Viimeiset sarakkeet antavat testisuureen ja $p$-arvon.

Kuten kahden otoksen t-testi, niin myös ANOVA voidaan muotoilla lineaarisen regression avulla. Määritellään lineaarinen regressiomalli \[\begin{equation} \tag{4.6} X_{ij} = \beta_0 + \beta_1 D_{i2} + \beta_2 D_{i3} + \cdots + \beta_{d-1} D_{id} + \varepsilon_i, \quad \varepsilon_i\sim N(0,\sigma^2), \end{equation}\] jossa $D_{ij}$, $j\in\{2, 3, \ldots, d\}$, on osoitinmuuttuja \[\begin{equation*} D_{ij} = \begin{cases} 0,& \text{jos havainto $X_{ij}$ ei kuuluu luokkaan $j$}, \\ 1,& \text{jos havainto $X_{ij}$ kuuluu luokkaan $j$}. \end{cases} \end{equation*}\] Eli jos ryhmiä on $d$ kappaletta, niin osoitinmuuttujia on $d - 1$ kappaletta. Nollahypoteesin $\mu_1 = \mu_2 = \cdots = \mu_d$ alla regressiomalli (4.6) ei ole tilastollisesti merkitsevä. Eli mallissa (4.6)

vakiotermi $\beta_0$ vastaa ensimmäisen ryhmän odotusarvoa $\mu_1$ ja
$\beta_{j-1}$ vastaa ryhmän $j$ ja ensimmäisen ryhmän odotusarvojen välistä erotusta $\mu_{j} - \mu_1$. Nollahypoteesin alla kaikki erotukset $\mu_j - \mu_1$ ovat nollia eli $\beta_1 = \beta_2 = \cdots = \beta_d = 0$.

Tiivistäen: ANOVA on siis vain F-testi lineaariselle regressiomallille (4.6)!

Esimerkki 4.9 (ANOVA ja lineaarinen regressio) Suoritetaan esimerkin 4.8 analyysi lineaarisen regression kautta. Olemme jo valmiiksi luoneet sopivasti muokatun taulukon choco_mod havaintoaineistosta. Voimme siis suoraan sovittaa kaavan (4.6) mukaisen lineaarisen regressiomallin.

choco_lm <- lm(weight ~ factory, data = choco_mod)

Komennon summary(choco_lm) tulosteesta voimme lukea F-testiä vastaavan-arvon. F-testiä vastaava p-arvo on sama kuin esimerkissä 4.8 käytetyn funktion aov() antama p-arvo.

summary(choco_lm)

## 
## Call:
## lm(formula = weight ~ factory, data = choco_mod)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.7190  -3.7660   0.2861   3.4947  12.9144 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     101.08086    1.27081  79.541  < 2e-16 ***
## factoryfactory2  -2.97829    1.79719  -1.657    0.102    
## factoryfactory3   0.01123    1.79719   0.006    0.995    
## factoryfactory4 -12.36215    1.79719  -6.879 1.49e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.683 on 76 degrees of freedom
## Multiple R-squared:  0.4562, Adjusted R-squared:  0.4347 
## F-statistic: 21.25 on 3 and 76 DF,  p-value: 4.272e-10

Huomautus 4.2 (Monivertailuongelma) Esimerkki 4.8 (tai vastaavasti esimerkki 4.9) näyttää, että ainakin yhdessä tehtaassa laatukriteeri liittyen suklaalevyjen painoihin ei täyty. Seuraava luonnollinen askel olisi selvittää, missä kaikissa tehtaissa laatukriteeri ei täyty. Tässä tapauksessa meidän täytyisi suorittaa useampi tilastollinen testi samalle havaintoaineistolle. Esimerkiksi voisimme suorittaa yhden otoksen t-testin nollahypoteesilla $\mu_j = 100\mathrm{g}$ jokaiselle tehtaalle $j\in\{1, \ldots, 4\}$ erikseen. Useamman tilastollisen testin suorittaminen samalle havaintoaineistolle aiheuttaa kuitenkin omat ongelmansa. Tätä niin kutsuttua monivertailuongelmaa (eng: multiple comparisons problem) käsitellään viimeisellä luennolla.

Huomautus 4.3 (Mistä nimitys "yksisuuntainen varianssianalyysi tulee"?) Varianssianalyysi nimitys tulee siitä, että odotusarvojen yhtäsuuruuksia vertaillaan käyttäen tunnuslukuja, jotka mittaavat tietynlaisia vaihteluja. Tarkalleen ottaen vertaamme ryhmien välistä vaihtelua ryhmien sisäisiin vaihteluihin.

Toisaalta etuliite “yksisuuntainen” littyy siihen, että ryhmien odotusarvoja verrataan yhden luokkamuuttujan suhteen – esimerkissä 4.8 luokkamuuttuja on tehdas (tehdas $1/2/3/4$). Jos odotusarvojen erilaisuutta tarkastellaan kahden luokkamuuttujan suhteen niin kyse on kaksisuuntaisesta varianssianalyysistä. Ylipäätään odusarvojen erilaisuutta voidaan tarkastella myös kolmen, neljän, … luokkamuuttujan suhteen.

Huomaa: jos $X\sim N(\mu_1, \sigma^2)$ ja $Y = X + \mu_2$, niin $Y\sim N(\mu_1 + \mu_2, \sigma^2)$.↩︎