Varianssianalyysi: tulosten graafinen esitys

Komento plot tuottaa samanlaisen diagnostisen grafiikan kuin regressioanalyysi (kuva 6) riippumatta siitä, onko käytetty komento ollut lm vai aov. Männyntaimiesimerkissä diagnostiikka näyttää mainiolta ja voimme olla malliimme tyytyväisiä.

R käyttää luokitellun aineiston graafisen esityksen oletuskuvana rinnakkaisia boxplot-kuvioita (kuva 8A), mutta varsinkin pienissä aineistoissa niiden tunnusluvut, mediaani ja etenkin kvartiilit ovat epäluotettavia. Alkuperäiset havaintoarvot esittävä pistekuvio (stripchart) on tässäkin esimerkissä parempi vaihtoehto (kuva 8B).

> par(mfrow=c(1,3))                      # 3 osakuvaa
> plot(istutus, pituus, ylab="Pituus (cm)", main="A: Boxplot")
> stripchart(pituus ~ istutus, vertical=T, method="jitter",
+  ylab="Pituus (cm)", main="B: stripchart")

Ekologiassa (kuten myös monilla biologian, lääketieteen ja käyttäytymistieteiden aloilla) on tapana esittää ANOVAn tulokset pylväsdiagrammina, jossa pylvään korkeus näyttää luokkakeskiarvon ja kuhunkin pylvääseen liitetään ns. virhejana. Ekologiassa on kaksi vahvaa koulukuntaa: toiset vaativat käytettäväksi virhejanana keskihajontaa, toiset taas keskiarvon keskivirhettä (pieni vähemmistö käyttää keskiarvon 95 % luottamusväliä). Kumpikaan käytäntö ei välttämättä ole hyvä. Kerrattakoon ensin eräitä perusasioita:

Koska varianssianalyysissä vertaillaan nimenomaan keskiarvoja, keskiarvon keskivirhe näyttää pinnalta katsoen luontevalta tunnusluvulta. Ongelmana on kuitenkin se, että yksittäisten keskiarvojen keskivirheiden avulla ei voi kovinkaan hyvin hahmottaa luokkakeskiarvojen erotusten (tai yleisempien kontrastien) luottamusvälejä tai ``merkitsevyyttä''. Ongelma ei poistu vaikka käytettäisiinkin luokkakeskiarvojen 95% luottamusvälejä. Toisaalta keskihajonnan sisältävä pylväskuvio on sekin data-mustesuhteeltaan huono ja epäinformatiivisempi kuin boxplot ja stripchart.

Vaikka seuraavassa annammekin ohjeet tällaisten pylväskuvioiden piirtämiseksi, emme suosittele niiden käyttämistä. Todettakoon, että yhdessäkään varsinaisessa tilastotieteen oppikirjassa, jossa tavanomaisia tilastokuvioita käsitellään, ei tätä kuviotyyppiä esiinny. Päin vastoin, silloin kun tilastotieteilijät ovat niitä kommentoineet, on arvio ollut kielteinen.

Kuva 8: A: Boxplot on R:n oletuskuva luokitteluille. B: Pienelle aineistolle niitäkin parempi on stripchart. C: Ekologit käyttävät yleensä pylväsdiagrammia ja virhejanaa osoittamaan keskiarvoa ja virhettä (tässä 95 % luottamusväli) sekä merkitsevät samalla kirjaimella pylväät, jotka ``eivät eroa merkitsevästi toisistaan''. Aineisto taimet.

\includegraphics[width=\textwidth]{anovabar.eps}

Pylväsdiagrammien piirtäminen ja varsinkin virhejanojen lisääminen on hieman mutkallista (kuva 8B):

> ka <- tapply(pituus,istutus,mean)        # Luokkakeskiarvot
> hajonta <- tapply(pituus,istutus,sd)     # Luokkahajonnat
> n <- tapply(!is.na(pituus),istutus,sum)  # n luokittain
> se <- hajonta/sqrt(n)                    # Keskivirheet
> se <- tapply(pituus,istutus, function(x) sqrt(var(x)/length(x)))
> ci <- se*qt(0.975,n-1)                   # Luottamusrajat
> mp <- barplot(ka, ylim=c(0,max(ka+ci)), ylab="Pituus (cm)",
+ main="C: Pylväät", col=NULL)
> arrows(mp,ka-ci, mp,ka+ci, angle=90, code=3, length=0.1)
> text(mp, 5, c("a","b","a"))

Meidän on ensin laskettava luokkakeskiarvot sekä virhejanan pituus. Tällä kertaa käytämme 95 % luottamusrajoja, lähinnä näyttääksemme kuinka ne lasketaan. Keskihajontaa varten on valmis funktio, mutta keskivirhe on laskettava itse ja sitä varten tarvitsemme lukumäärät luokissa (n). Keskivirheen voi laskea myös ilman välivaiheita, käyttämällä nimetöntä funktiota tapply-komennossa. Hieman yllättävän näköinen tapa laskea havaintojen luvut luokittain (n) takaa, että tulos ottaa huomioon mahdolliset puuttuvat tiedot muuttujassa pituus. Puuttuvia tietoja ei tällä kertaa ole, joten ilmeisempi table(istutus) olisi tuottanut samat tulokset. Komento qt(p,df) antaa kriittisen t-arvon. Käytimme todennäköisyyttä p = 0.975, joka jättää häntään 2.5 %, ja näin kaksisuuntaisena antaa 95 % luottamusrajat.

Komento barplot tuottaa pylväsdiagrammin ja ylim takaa että y-akseli skaalataan niin, että virhejanatkin mahtuvat kuvaan. Yleensä R:n oletusgrafiikka on mustavalkoista, ja värit on lisättävä komennoilla, mutta barplot tuottaa värikuvia, jotka eivät sovi tähän mustavalkoiseen oppaaseen, joten olen poistanut värit (col=NULL). Komento barplot palauttaa pylväiden keskiviivojen x-koordinaatit (mp), joita käytämme komennossa arrows sijoittamaan virhejanat x-suunnassa. Kulma 90o määrittelee nuolenpään kulmaksi suoran, ja code=3 kertoo, että nämä suorat ``nuolenpäät'' piirretään kumpaankin päähän janaa.

Ekologien yleensä käyttämät kirjaimet (a, b,...) tarkoittavat, että luokat ovat parittaisessa vertailussa merkitsevästi erilaiset, mikäli niissä ei ole yhteistä kirjainta. Merkitsevä ANOVA-tuloksemme näyttäisi aiheutuvan yhdestä ainokaisesta merkitsevästä erosta: ruukkutaimet ovat pitempiä kuin muut. Tällainen post hoc vertailu on tavallinen käytäntö ANOVAn jälkeen ja useimmissa tilasto-ohjelmistoissa on laaja valikoima vertailumenetelmiä. Ei kuitenkaan R:ssä, missä voimme käyttää funktioita pairwise.t.test6.1.4, sivu [*]) tai funktiota TukeyHSD. Erillisissä paketeissa on todennäköisesti lisää vaihtoehtoja, mutta R:n kehittäjät suhtautuvat monivertailuihin penseästi, koske niitä yleensä käytetään väärin ja ne ovat hyväksyttäviä vain harvinaisissa tapauksissa. ANOVA olettaa että ainoa ero luokkien välillä on keskiarvo ja kaikkien luokkien sisävaihtelu on sama. Erimittaisten hajontajanojen piirtäminen on vastoin ANOVAn oletusta ja näin ollen testimme ja kuvamme ovat ristiriidassa. Samaten monivertailut tehdään jälkikäteen ja ne eivät välttämättä ole yhteneväisiä ANOVAn kanssa: on mahdollista saada merkitsevä ANOVA ilman ainokaistakaan parittain merkitsevää eroa tai epämerkitsevä ANOVA ja useita parittain merkitseviä eroja.


Jari Oksanen 2003-01-21