Tämän osan verkkomateriaali on vähän aiempia lyhyempi ja filosofisempi, mikä myös jättää aikaa tarvittaessa ottaa kiinni aiempien viikkojen sisältöä. Ylipäätään verkkomateriaalin on tarkoitus toimia linkkinä yleisen R-opetuksen ja kurssikirjan välillä, joka taas ei nojaa R:n. Lisäksi pyrin korostamaan verkkomateriaalissa Howellin tilastotekniikoiden yleistyvämpiä elementtejä, eli tilastomallinnus näkökulmaa. On olennaista, että olet lukenut usean vertailun korjauksista ja kontrasteista (esim. Howell, kpl 12) ennen tätä verkkomateriaalin OSAa 3. Myös aiemmista R-materiaaleista on iloa, mutta tämän viikon R-videoluennon asioita emme suoraan hyödynnä alla vielä.
3.1 Malleja, ei koneita
Tieteellisen kehityksen kiihdyttyä eksponentiaalisesti, on yhä tärkeämpää ymmärtää yleisiä periaatteita nopeasti vanhenevien yksittäisten “temppujen” sijaan – siis ymmärtää yleisperiaatteet temppujen taustalla. Howellinkin toisinaan suosimista vuokaavioista voi syntyä osin virheellinen ajatus, että tilastolliset testit olisivat tiettyyn käyttötarkoitukseen suunniteltuja koneita tai “temppuja”. Ne eivät ole koneita, jotka muuttavat rahan tai aineiston halutuksi vastaukseksi. Ne ovat haluajan ajatusten operationalisoituja määritelmiä. Jos ostat tilastollisen analyysin tilastotieteellisen konsultaation sijaan, ostat jonkun toisen ajatuksen ja käsityksen todellisuudesta. Näin ei käy, jos ostat porakoneen puutöitäsi varten (ainakaan verrattavissa määrin). Tässä OSAssa katsahdamme tähän Howellin esittelemien tekniikoiden kääntöpuoleen OSAn 2 tietojen valossa. Howellkin kyllä asiaa sivuaa (ainakin 7 painoksessa), mutta R ja viimeaikaisemmat tiedekeskustelut yhdessä sallivat meidän tehdä ymmärryksestämme aavistuksen Howellia konkreettisempaa. Lisäksi katsomme kuinka lineaarisia kontrasteja voi R:ssä määrittää.
3.2 Korjausten kääntöpuolia
Edellisessä osassa testasimme ryhmäkeskiarvojen varianssin olemassa oloa. Sellaisen löytyessä, tutkija usein mielellään sanoisi jotain itse ryhmien suhteesta toisiinsa. Toisinaan koko ryhmien välisen varianssin kysymys on toissijainen ja tutkimuksessa vain on useita ryhmiä ja useita näitä koskevia kysymyksen asetteluita, mistä voi seurata tarve tehdä useita tilastollisia testejä. Valitettavasti perinteisimmät tilastollisen epävarmuuden kontrolloinnin menetelmät on luotu yhtä isoloitua testiä silmällä pitäen, eivätkä niiden tulokset sellaisenaan yleisty usean vertailun tilanteeseen. Siksi menetelmiä usein “korjataan” usean vertailun käyttötarvetta silmällä pitäen. Korjauksetkaan eivät kuitenkaan ole ongelmattomia. Lähdetään tutustumaan tilanteeseen tarkemmin.
Ajatellaan, että yhden kahden ryhmän (\(n = 30\) per ryhmä) vertailun sijaan haluamme tehdä viisi kappaletta riippumattomia vertailuja. Vertailut ymmärretään kuitenkin samaksi tutkimuskokonaisuudeksi jossain mielessä. Pyrimme ymmärtämään tilannetta simuloimalla sitä toistuvasti R:ssä (10000–50000 kertaa). Voit palata alla olevaan koodinpätkään tarkemmin luettuasi ensin sen alta asian merkityksestä, mutta kiinnitä alkuun huomiosi vain koodin lopputulokseen, eli siihen millaista aineistoa se tuottaa ja mitä siitä koodiosuuden jälkeen kirjoitan.
#set.seed(22012026)Nr <-10^4# Monte Carlo -toistojen lukumäärän <-30# Otoskoko per ryhmäNt <-5# Testien lukumäärä per tutkimustestitulokset <- xm <- ym <-matrix(0,Nr, Nt) # p-arvo-tulosvektorin alustusfor (i_testi in1:Nt) {for (i in1:Nr) { x <-rnorm(n) # simuloitu n-kokoinen otos ryhmälle 1 y <-rnorm(n) # simuloitu n-kokoinen otos ryhmälle 2 a <-t.test(x, y) # t-testin tiedot apumuuttujaan a testitulokset[i, i_testi] <- a$p.value # p-arvo xm[i, i_testi] <- a$estimate[1] # ensimmäinen keskiarvo ym[i, i_testi] <- a$estimate[2] # toinen keskiarvo }}colnames(testitulokset) <-paste0("Testi", 1:Nt)rownames(testitulokset) <-paste0("Iteraatio", 1:Nr)head(testitulokset)
Olemme nyt luoneet matriisin (aineiston) p-arvoja, jossa jokainen rivi vastaa erillistä “tutkimusta”, joka perustuu muista tutkimuksista riippumattomalle otokselle. Matriisin sarakkaille taas on tutkimuksen tekemät viisi testiä. Esimerkiksi rivin kolme (Iteraatio3) toinen arvo (sarake 2) vastaa kyseisessä tutkimuksessa (simulaatioiteraatiossa) tehtyä toista tilastollista testiä (Testi2), edustaen siitä saatua p-arvoa (0.2802267). Olennaista jatkon kannalta on, että kaikki aineisto on simuloitu nollamallista – jokainen p-arvo siis heijastelee pelkkää otantavaihtelua. Jokainen arvo alittaa luvun 0.05 todennäköisyydellä \(\alpha = 0.05\) ja tällainen tapahtuma tässä aina edustaa tyypin 1 päättelyvirhettä (nollahypoteesin hylkäämistä, kun se on totta).
Testiteoria siis sanoo, että käyttäessämme päätöksentekokriteerinä nollahypoteesin hylkäämistä silloin, kun \(p < 0.05\), tyypin 1 virheitä, eli vääriä hylkäyksi pitäisi tulla keskimäärin vain 5 %. Tarkistetaan nyt tämä teoreettinen ennuste laskemalla keskiarvo kaikista kerroista, joilla simuloitu p-arvo alitti arvon 0.05:
mean(testitulokset <0.05)
[1] 0.05108
Toden totta, noin osuus 0.05, eli 5 %, kaikista testikerroista toteutti ehdon \(p < 0.05\) nollahypoteesin pätiessä. Jos kuitenkin ajattelemme, että yllä nähdyt kunkin rivin 5 arvoa vastaavat saman tutkimuksen eri testejä, todennäköisyys tehdä tyypin virhe jossain testeistä ylittää 0.05, eli asetetun \(\alpha\)-tason. Kuten kirja keskustelee, teoreettinen ennuste on:
\[
\underbrace{1 - \underbrace{(1 - \overbrace{0.05}^{P(\text{virhe 1:ssä testissä})})^5}_{P(\text{ei virheitä 1:ssäkään)}}}_{P(\text{1 tai useampi virhe})} \approx 0.226.
\] Yhden tai useamman virheen todennäköisyys on siis vastatodennäköisyys (komplementtitapahtuman tn.) sille, että viidestä peräkkäisestä riippumattomasta testistä ei mistään synny virhettä, ja riippumattomien tapahtumien todennäköisyys saatiin yksittäisten todennäköisyyksien tulona. Tarkistetaan, että simulaatiomme todella on arvio tästä luvusta tai, kääntäen, että simulaatio varmentaa yo. logiikkamme. Se tehdään laskemalla keskiarvo riveistä, joilla jokin p-arvoista ylittää arvon 0.05:
mean(apply(testitulokset <0.05, 1, any))
[1] 0.2319
Oikealta näyttää.
Jos moni pitääkin normaalia 0.05 virhetasoa (1 per 20 testiä) nykyisin löyhänä, vielä harvempi olisi valmis rutiininomaisesti hyväksymään virhettä lähes joka neljännessä tutkimuksessaan. Kirjasta opimme, että Holmin korjaus on tehokas tapa suojautua tältä perhekohtaisen virhetodennäköisyyden (family-wise error rate, FWER tai FW) inflaatiolta. Todetaan asia R:ssä:
Funktio p.adjust inflatoi p-arvoja siten, että perhekohtainen (koko rivin) tyypin 1 päättelyvirheen todennäköisyys pysyy 5 %:ssa. Homma vaikuttaisi olevan hanskassa siis, mutta Howellkin varoittelee, että prosessissa menetetään tilastollista voimaa. Sittemmin, paljon on tapahtunut alalla (mm. replikaatiokriisi-keskustelua ymv.) ja yleinen ymmärrys on parantunut koskien voiman menettämisen merkityksiä. Katsotaan seuraavaksi mitä tapahtuu niille tilastollisesti merkitseville, mutta väärille, löydöksille, jotka seulasta pääsivät läpi – tarkoituksella, koska ~5% virhetiheys itse hyväksyttiin). Piirretään simulaatiojakaumatietoa (eli nollajakaumatietoa) kaikista absoluuttisista efektikoista, tilastollisesti merkitsevistä sellaisista ja viiden testin perheessä merkitsevistä efektikoista:
Vaikka nollahypoteesin pätiessä keskiarvo ero \(\bar{Y} - \bar{X}\) onkin keskimäärin nolla, itseisarvo \(|\bar{Y} - \bar{X}|\) ei ole keskimäärin nolla. Itseasiassa mediaanin tulisi olla puolessa välissä normaalijakauman \(N(0, 2/n)\)positiivista puoliskoa, koska negatiivinen puolisko peilautuu positiiviselle itseisarvon myötä (voit verrata lukuja qnorm(0.75, sd = sqrt(2/n)) ja median(efektikoot_kaikki) ja miettiä asiaa kuvan avulla). Jos kuitenkin ehdollistamme löydöksen sille, että se oli tilastollisesti merkitsevä, tämän joukon mediaani on puolessa välissä viimeisintä 2,5 prosenttia jakauman hännästä, eli arvoltaan qnorm(0.975 + 0.025/2, sd = sqrt(2/n)) = 0.5787277. Teoreettinen arvio ei ole tarkka, koska t-jakauma ei täysin vastaa normaalijakaumaa, mutta jälleen hämmästyttävän lähellä empiiristä arviota median(efektikoot_merk) = 0.5769533 (ks. myös yo. kuva). Koska Holmin proseduuri korjaa p-arvoja ylöspäin, myös ryhmäerohavainnon on oltava suurempi saavuttaakseen merkitsevyyden. Kääntäen, jos ehdollistamme “löydökselle”, että tulos oli tilastollisesti merkitsevä, FWER-korjatut virheelliset tulokset ovat efektikooltaan suurempia virheitä kuin korjaamattomat tulokset. Tässä on olennaista ymmärtää, että yllä teimme 10 000 replikaatiota, kun tyypillisesti kirjallisuutta lukiessa tai tutkimusta tehdessämme näemme vain yhden. Jos se yksi sattui olemaan tilastollisesti merkitsevä, harva osaa spontaanisti ajatella mitä olisi tapahtunut 9 999 muulle mahdolliselle eri otokseen perustuvalle todellisuuskululle. Olisiko tutkimus kaikissa tapauksissa julkaistu ja silmiemme edessä vai vaikuttaako kenties tapahtuma “vertailu oli merkitsevä” asiaan, jne?
Saavumme keskeisen opetuksen äärelle. Olemme maksaneet hinnan perhekohtaisen virheen kontrolloinnista. Hinta on, että virheen sattuessa, se on tavallista virhettä suurempi. Toinen maksettu hinta koskee efektin koon sijaan sen totuusarvoa. Paradoksaalisesti, yritys olla konservatiinen yhden testiperheen kohdalla maksaa hintaa sen jäseniä koskevien päätelmien luotettavuudessa. OSAn 2 opein voimme esimerkiksi laskea PPV-arvon, eli todennäköisyyden, jolla merkitsevä löydös tuli vaihtoehtoisesta hypoteesista (olisi “todellinen”). Oletetaan esimerkin vuoksi, että priori-vetosuhde on 80-20 nollahypoteeslle ja tilastollinen voima kohdeilmiölle on 0.80 ennen korjausta. Voimme laskea PPV:n ennen ja jälkeen korjauksen seuraavasti:
voima <-0.80# tutkimuksen oletettu tilastollinen voimaR <-20/80# pre-tutkimus vetosuhde 20-80merkitsevyystaso <-0.05# tyypillinen tilastollisen testin merkitsevyystaso# post-tutkimus uskomus H1:n todennäköisyydestävoima*R/(voima*R + merkitsevyystaso)
[1] 0.8
# esimerkkitilanteessa oletusvoiman antava efektikoko olisi:efekti <- pwr::pwr.t2n.test(n1 = n, n2 = n, power =0.8)$defekti
[1] 0.7356292
# uusi voima, jos Bonferroni-korjataan merkitsevyystasoa usealle testille olisi:voima <- pwr::pwr.t2n.test(n1 = n, n2 = n, d = efekti, sig.level = merkitsevyystaso/Nt)$powervoima
[1] 0.5760157
# post-tutkimus uskomus (juuri tämän vertailun) H1:n todennäköisyydestävoima*R/(voima*R + merkitsevyystaso)
[1] 0.7422733
Posterioriuskomme todelliseen ilmiöön tilastollisesti merkitsevän keskiarvoeron taustalla siis laski yli viidellä prosenttiyksiköllä, koska “uhrasimme” tilastollista voimaa neljän kyseiseen keskiarvoeroon liittymättömän testin virheen kontrollointiin. Tilannetta voisi kuvata sanonnalla, että tilastollisessa päättelyssä ei ole “ilmaisia lounaita” (toisinaan tilastotieteellisissä artikkeleissa esitetäänkin “no-free-lunch teoreemia”). Soveltajan kannalta tämä tarkoittaa, että (1) on tarkkaan harkittava minkä epävarmuuden kontrollointi on kriittistä milloinkin, (2) varmuutta ei aina voi saavuttaa ja (3) sen houkutukseen lankeaminen yhtäällä voi viedä turmioon toisaalla.
Yllä olevan problematiikan yksi hallintakeino on väärien löydösten osuuden (false discovery rate, FDR) kontrollointi. Kun \(V\) on virheellisten hylkäyksien lukumäärä useassa testissä ja \(R\) on kaikkien hylkäysten lukumäärä, FDR pyrkii kontrolloimaan odotusarvoa:
\[
FDR = \text{E}[V/R|R>0] \cdot P(R>0).
\]
Nähdään, että FDR-proseduurit, kuten Benjamini-Hochberg-proseduuri, pyrkivät kontrolloimaan virheen todennäköisyyttä (suhteellista osuutta), kun löydös on sattunut. Näin ne huomioivat yllä mainittuja huolia ja ovat tulleet yhä tärkeämmiksi saatavilla olevien aineistomäärien (ja siis mahdollisten kysymysten/testien) kasvaessa digitaalisten teknologioiden vauhdittamana.
Kaikkein tärkein kontrolli ei kuitenkaan ole tilastollinen, ainakaan täysin. Monet tilastotieteilijät ja tilastotieteen soveltajat ajattelevat, että tärkein kontrolli on omistautua merkityksellisen kohdesuureen estimoinnille: siis useiden testien “roiskimisen” sijaan, miettiä tarkasti mikä suure vastaa kaikkein keskeisimpään tarpeeseen, ja estimoida se mahdollisimman tarkasti ja mahdollisimman pienellä virheen mahdollisuudella. Tämän asian arvioiminen on suurelta osin myös alakohtaisen substanssin kysymys. Alla esitellään joitain työkaluja luovemman kysymisen toteuttamiseen.
3.3 Lineaariset kontrastit
3.3.1 Taustoitusta
Lineaariset kontrastit auttavat yllä esiteltyjen ongelmien kanssa. Usein nimittäin tutkija haluaa käsitellä monimutkaisempaa todellisuutta taustalla, vaikka kohdekysymys koskisikin verrattain suoraviivaista vertailua. Useiden välikysymyksenasetteluiden sijaan, hän voi sisällyttää taustatodellisuuden malliin, jonka puitteissa kohdevertailu tapahtuu. Otetaan tästä ensin lyhyt esimerkki.
set.seed(54321)n <-40# havaintojen lukumääräx <-c(rep(1,n/2), rep(0,n/2)) # selittäjä x saa arvoja 1, 1, 0, 0,...,z <-rep(c(rep(0,n/4), rep(1,n/4)),2) # selittäjän z saadessa 1, 0, 1, 0,...y <-0.4*x + z +rnorm(n) # todellinen vastemuuttuja on x + z + e, e ~ N(0,1)fit <-lm(y ~ x) # sovitetaan lineaarinen mallisummary(fit)$coefficients # tarkistetaan estimoidut kulmakertoimet ja p ja t
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.2895803 0.2907420 0.9960042 0.32554844
x 0.7793824 0.4111713 1.8955175 0.06564893
# Huomataan, että x-kertoimen p- ja t-arvot ovat samat kuin ao. t-testissäkin:ttul <-t.test(y[x==1], y[x==0])c(ero =-diff(ttul$estimate), arvo = ttul$statistic, arvo.p = ttul$p.value)
ero.mean of y arvo.t arvo.p
0.77938245 1.89551746 0.06753659
# Mutta, eivät kuitenkaan z:n ollessa mallissa:fit2 <-lm(y ~ x + z)summary(fit2)$coefficients
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.3045369 0.3187534 -0.9553998 0.345576449
x 0.7793824 0.3680647 2.1175148 0.040999091
z 1.1882344 0.3680647 3.2283304 0.002610267
Yllä totesimme, että yksinkertaisen regressiomallin kulmakerroin todella vastaa t-testiä, jos selittäjä on kaksiluokkainen. Näin ei kuitenkaan ole, jos mallissa on muita selittäjiä (z) mukana. Muita selittäjiä usein lisätään kohdemuuttujan (x) vaikutuksen “vakioimiseksi” muille tekijöille (z), ja tämä todella vastaa lineaarista approksimaatioita tämän tekstin OSAn 1 vakiointimenetelmälle. Nyt tekijä z ei kuitenkaan vaikuttanut selittäjään x, eikä siten voinut olla sen sekoittava tekijä. Se kuitenkin vaikutti malliin, jota ajattelimme. Nimittäin ensimmäisessä mallissa ja t-testissä x:stä riippuva vastemuuttuja oli
\[
y = x + z + \varepsilon =: x + \varepsilon_1,
\]
kun taas jälkimmäisessä testissä se (tässä lineaarisessa tapauksessa) oli
\[
y|z = x + \varepsilon =: x + \varepsilon_2.
\]
Yllä olevista esitystavoista huomaamme, että varianssitermi \(Var(\varepsilon_1) = Var(z + \varepsilon) = Var(z) + Var(\varepsilon)\) on suurempi kuin \(Var(\varepsilon_2) = Var(\varepsilon)\). Jälkimmäisen tapauksen pienempi “kohina” vahvistaa x:n “signaalia” tilastollisessa testissä, koska testi vertaa signaalia kohinaan – tämän vuoksi jälkimmäisessä vertailussa saadaan tilastollisesti merkitsevä tulos, vaikkei ensimmäisessä saatukaan. Tässä on olennaista, ettei tutkijaa välttämättä kiinnostunut testata tekijää z. Se edusti t-testiä monimutkaisempaa todellisuutta, jonka tutkija halusi huomioida, kuitenkaan testaamatta sitä koskevaa hypoteesia. Vakiointi siis olennaisesti poisti termin \(Var(z)\) virhevaihtelu-tulkinnan piiristä.
Kirjassa Howell näyttää kuinka lineaariset konstrastit sallivat tutkijan kysyä “alakysymyksiä” mallista, altistamatta itseään lukuisista vertailuista koituville epävarmuuden hallinnan ongelmille. R-ohjelmaa hyödyntävät tilastotutkijat kuitenkin melko harvoin nojaavat lineaarisiin kontrasteihin, ja sille on syynsä. Lineaarinen kontrasti usein implementoi ajatuksen kahden mallin vertaamisesta toisiinsa ja mallien vertailuun on yleistyvämpiäkin työkaluja. Tieteen uudistuminen vahvempiin ja yleistyvämpiin abstraktioihin nojaamalla on varsin tyypillistä ja tässä tapauksessa näkyy siten, että R:ssä historiallisesta terminologiasta juontava anova-niminen funktio ei tee varianssianalyysia (joka oli mm. aov-funktion takana), vaan varsin yleisesti (eri paketeissa) vertaa kahta mallia toisiinsa (usein ns. uskottavuusosamäärä-testillä). Tyypillisesti kontrasteilla ja mallivertailuilla saavutetaan samoja asioita, mutta kontrasteilla voi olla helpompi implementoida monimutkaisia ANOVA-vertailuita, kun taas mallivertailut ovat huomattavasti yleisempi ajattelutapa, joka ei rajoitu parametrien lineaarisiin kuvauksiin.
3.3.2 Mallien vertailuista
Yllä nähty x:n kulmakerroin, jota voidaan kutsua nimellä \(\beta_x\), suoraan edustaa ns. käsittelykontrastia tai hoitokontrastia (treatment contrast). Kun muistetaan, että tässä tapauksessa yllä nähty ensimmäinen regressiomalli edusti ryhmään X (eli ryhmään x == 1) kuuluvia kertoimilla \(\beta_0 + \beta_x\) ja siihen kuulumattomia vakiokertoimella \(\beta_0\), nähdään, että \(\beta_x = (\beta_0 + \beta_x) - \beta_0\) on x-ryhmäläisten ja ei-x-ryhmäläisten käsittelyjen keskiarvoero. ANOVA-malli taas edustaa samoja ryhmäkeskiarvoja parametrein \(\mu + \tau_x\) ja \(\mu + \tau_{\text{ei } x}\), missä \(\mu\) on koko aineiston yleiskeskiarvo. Howellin notaatiolla hoitokontrasti olisi lineaarinen kontrasti \(a = (1,-1)\), joka antaisi siis kulmakerrointa \(\beta_x\) vastaavan vertailun:
Koska parametri \(\beta_x\) tässä jo on kontrasti, eli vertailu, itsessään, ko. vertailun testaaminen on sama asia kuin hypoteesin \(H_0: \beta_x = 0\) testaaminen. Kyseisen parametrin ja kontrastin testaaminen taas on sama asia kuin kahden mallin vertailu uskottavuusosamäärän testillä, joista toinen on sovitettu pakottamalla ehto \(\beta_x = 0\) ja toinen \(\beta_x\):n vapaa estimoituminen sallien. Ehto \(\beta_x = 0\) on käytännössä sama asia kuin selittäjän x jättäminen pois mallista. Näin ollen R:ssä voi toteuttaa tämän kontrastivertailun mm. seuraavasti:
anova(lm(y ~1), lm(y ~ x))
Analysis of Variance Table
Model 1: y ~ 1
Model 2: y ~ x
Res.Df RSS Df Sum of Sq F Pr(>F)
1 39 70.318
2 38 64.244 1 6.0744 3.593 0.06565 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Yo. p-arvon nähdään vastaavan aimpaa kulmakertoimen testiä.# F-arvo taas vastaa kulmakertoimen t-arvon neliötä:summary(lm(y ~ x))$coefficients["x", "t value"]^2
[1] 3.592986
Yksinkertaisen regressiomallin kulmakertoimen testi on siis sama asia kuin ko. regressiomallin selityskyvyn vertaaminen pelkällä vakiolla \(\beta_0\) vastemuuttujaa ennustavaan malliin. Vertailupistemalli voi kuitenkin olla monimutkaisempikin, kunhan kohdemalli vain on vielä monimutkaisempi. Katsotaan mallivertailu yo. tilanteessa, jossa taustatekijä z huomioitiin:
fit_x_and_z <-lm(y ~ x + z)fit_z_only <-lm(y ~ z)anova(fit_z_only, fit_x_and_z)
Analysis of Variance Table
Model 1: y ~ z
Model 2: y ~ x + z
Res.Df RSS Df Sum of Sq F Pr(>F)
1 38 56.199
2 37 50.124 1 6.0744 4.4839 0.041 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Nämä esimerkit olivat tarkoituksellisen yksinkertaisia pedagogisista syistä. Tutkimusmaailmassa yllä nähty mallien vertailu osoittautuu kuitenkin yllättävän joustavaksi tavaksi muotoilla haluamiaan hypoteeseja, mutta anova-funktio toimii (oikein!) tavallisille regressiomalleille vain, kun mallit ovat toisiinsa nähden sisäkkäisiä (nested). Toisin sanoen, anova funktion ensimmäisen argumentin malli on voitava muodostaa sitomalla toisen argumentin mallin parametreja. Jos nyt haluaisimmekin verrata ryhmään X kuulumisen vaikutusta ryhmään Z kuulumisen vaikutukseen, X:n poissulkeva malli ei ole Z:n poissulkevan mallin sisällä, eikä toisin päin (huomaa myös, että neljännes havainnoista kuuluu kumpaankin ryhmään; totea se ajamalla esim. sum((x==1)&(z==1))). Nyt olisi hyödyllistä käyttää lineaarista kontrastia \(a = (1,-1)\) suoraan regressiokertoimiin (ANOVAn ryhmäkeskiarvojen sijaan), jolloin \(\sum_{i \in {x,z}} a_i \beta_i = \beta_x - \beta_z\). [Itse asiassa pärjäisimme kyllä mallivertailuin, jos vain muuttaisimme regressiomallin selittäjäaineistoa (design-matriisia) jo ennen mallin sovitusta, mutta oletetaan tässä kurssisisältöjen hengessä haluavamme “jälkiasentaa” hypoteesin lineaarisin kontrastein].
3.3.3 Kontrastit
Kun kontrasteja käytetään regressiokertoimiin, puhutaan yleisestä lineaarisesta hypoteesista (general linear hypothesis) regressiomallille. Yleinen lineaarinen hypoteesi on mikä hyvänsä oletus \(\sum_{i=0}^k c_i \beta_i = d\), missä \(k\) on selittäjien lukumäärä, \(c_i\) ovat vakioita ja \(d\) on vakio (yleinen lineaarinen hypoteesi voi olla myös matriisimuotoa, \(C\beta = d\)). Tässä tapauksessa haluamme valita \(c\)-vektoriksi \([0, 1, -1]\) ja \(d=0\), jolloin \(\sum_{i=0}^2 c_i \beta_i = \beta_1 - \beta_2\), eli \(\beta_x - \beta_z\), jos x on mallissa ensimmäinen regressioennustaja ja z toinen. Varsinaisen testin voi R:ssä suorittaa paketilla gmodels, seuraavasti:
Test of General Linear Hypothesis
Call:
gmodels::glh.test(reg = fit_x_and_z, cm = matrix(c(0, 1, -1),
nrow = 1))
F = 0.617, df1 = 1, df2 = 37, p-value = 0.4372
On olemassa myös yleisempi menetelmä nimeltään Delta-metodi. Siinä hypoteesin ei tarvitse olla edes lineaarinen, mutta tarkistetaan esimerkin vuoksi, että tulos on suurin piirtein sama lineaariselle hypoteesille:
b <-coefficients(fit_x_and_z) # poimitaan kulmakertoimetb <- b[2:3] # jätetään vakio poisv <-vcov(fit_x_and_z)[2:3,2:3] # kertoimien varianssi-kovarianssi-matriisise <- msm::deltamethod(~(x1-x2), b, v) # x1 = x ja x2 = z tässäp_val <-as.numeric( 2*(1-pnorm(abs(b[1]-b[2])/se)) ) # Waldin p-arvoc(b_xz = b[1] - b[2], se = se, p = p_val)
b_xz.x se p
-0.4088520 0.5205221 0.4321809
Yleisellä lineaarisella hypoteesilla voi helposti tehdä monenlaisia vertailuja, mutta ANOVA-kirjallisuudessa kontrasteja usein asetetaan faktorimuuttujan tasoille/ryhmäkeskiarvoille. Kuten todettua, kyse on pitkälti samasta asiasta, mutta halutessaan tähänkin voi tutustua. Itse luultavasti mieluummin vain koodaisin R:ssä aineiston sellaiseen muotoon, että regressiomalli suorempaan antaa halutun vastauksen, mikä olisi lukijaystävällisempi esitystapa ja myös vähemmän hämmentävää itselleni ja siten vähemmän virhealtista. Koska en voi kuitenkaan varmaksi sanoa, etteikö lukijoillani joskus syntyisi tarvetta (tai tarpeen kokemusta) faktorikontrasteille, katsotaan vielä lopuksi miten perinteisiä ANOVA-kontrasteja R:ssä voi halutessaan muodostaa:
# Ilmaistaan x- ja z-muuttujat yhtenä faktorimuuttujana ja taulukoidaan tasot:xf <-rep("", n)xf[(x==0)&(z==0)] <-"ref"xf[(x==1)&(z==0)] <-"x"xf[(x==0)&(z==1)] <-"z"xf[(x==1)&(z==1)] <-"xz"xf <-factor(xf, levels =c("ref", "x", "z", "xz"))table(xf)
xf
ref x z xz
10 10 10 10
# Verrataan x- ja z-muuttujia keskiarvokontrastein:( testaa <- gmodels::make.contrasts(c(0,1,-1,0)) ) # from human to R readable
summary(lm( y ~ xf, contrasts=list(xf=testaa ))) # alla C1 on haluttu kontrasti
Call:
lm(formula = y ~ xf, contrasts = list(xf = testaa))
Residuals:
Min 1Q Median 3Q Max
-2.85037 -0.66682 0.01306 0.59681 2.78065
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6793 0.1778 3.821 0.000507 ***
xfC1 -0.4089 0.5029 -0.813 0.421528
xfC2 -0.1763 0.3556 -0.496 0.623033
xfC3 1.5380 0.3556 4.326 0.000115 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.124 on 36 degrees of freedom
Multiple R-squared: 0.3527, Adjusted R-squared: 0.2988
F-statistic: 6.539 on 3 and 36 DF, p-value: 0.001208
# Voidaan kerralla testata useita custom-kontrasteja myös:cmat <-rbind( "1 vs 4"=c(-1, 0, 0, 1),"1+2 vs 3+4"=c(-1/2,-1/2, 1/2, 1/2),"1 vs 2+3+4"=c(-3/3, 1/3, 1/3, 1/3))(testaa <- gmodels::make.contrasts(cmat))
1 vs 4 1+2 vs 3+4 1 vs 2+3+4
V1 -1.110223e-16 4.718448e-16 -0.75
V2 2.081668e-16 -1.000000e+00 0.75
V3 -1.000000e+00 1.000000e+00 0.75
V4 1.000000e+00 5.551115e-16 -0.75
# Kokeile halutessasi ajaa ao. kommentoitu rivi (ilman risuaitaa):# summary(lm( y ~ xf, contrasts=list(xf=testaa )))
3.4 Yhteenveto
Tilastollisessa testaamisessa (1) on tarkkaan harkittava minkä epävarmuuden kontrollointi on kriittistä milloinkin, (2) varmuutta ei voi saavuttaa ja (3) sen houkutukseen lankeaminen yhtäällä voi viedä turmioon toisaalla.
Lineaariset kontrastit ovat tapa kysyä monimutkaisia hypoteeseja pilkkomatta niitä useaan testiin, eli tapa hallita epävarmuutta, kun tietää mitä haluaa.
Tilastollisten mallien vertailut ovat yleisempi tapa hallita epävarmuutta monimutkaisemmissa tilanteissa ja yksittäisenkin mallin parametreille voi muodostaa myös epälineaarisia hypoteeseja.
