Yksinkertainen regressio

Ranta et al. (1989, Esimerkki 10.1) käyttävät yksinkertaisen regressioanalyysin esimerkkinä aineistoa, jossa tutkitaan ravinnon määrän vaikutusta vesiliskon lisääntymiseen. Kohdemuuttujana on vesiliskon tuottamien munien määrä (muuttuja munat), selittävänä tekijänä vesiliskojen ruokintakerralla saama surviaistoukkien lukumäärä (muuttuja ruoka). Näin pieni aineisto on helppo tallentaa itsekin suoraan:⁸

> munat <- c(4, 1, 6, 2, 5, 11, 7, 11, 9, 15)
> ruoka <- seq(along=munat)

Sovitamme regressiomallin ja pyydämme yhteenvedon tuloksista (olen poistanut osan summaryn riveistä kaikissa esimerkeissä):

> lisko.reg <- lm(munat ~ ruoka)
> summary(lisko.reg)

Call:
lm(formula = munat ~ ruoka)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    0.400      1.750   0.229  0.82489
ruoka          1.218      0.282   4.320  0.00254

Residual standard error: 2.561 on 8 degrees of freedom
Multiple R-Squared:   0.7,      Adjusted R-squared: 0.6625
F-statistic: 18.67 on 1 and 8 degrees of freedom,       p-value: 0.002545

Regressioanalyysin tulos lisko.reg on luokan lm olio. Komento summary tuntee luokan, ja osaa esittää tärkeimmät tiedot. Komento print (jota kutsuu epäsuorasti nimen lisko.reg kirjoittaminen) antaisi niukemmat tiedot. Komennolla attributes(lisko.reg) saamme luettelon olion attribuuteista. On myös joitain erityismenetelmiä eristää joitain tiettyjä luokan attribuutteja. Esim. coef(lisko.reg) antaa pelkät regressiokertoimet. Ainoan selittävän tekijän ruoka regressiokerroin on 1.22: näin paljon munamäärä kasvaa kun ruokintaa lisätään yhdellä surviaisella. R tulostaa myös kertoimien keskivirheen (Std. Error) ja t-arvon, joka on estimaatti jaettuna keskivirheellään. Kertoimen merkitsevyyden arviointi perustuu t-testiin (§6.1). Ruokinta lisää merkitsevästi jälkeläistuotantoa ( p = 0.0025), mikä ei ole järisyttävä tieto. R saattaa tulostaa myös tähdet (*) ``merkitsevien'' kertoimien jälkeen, mutta olen tarkoituksella poistanut tämän toiminnon (katso options), koska tiukat tähtirajat ovat arveluttavia -- onko p = 0.049 todellakin yhden tähden arvoinen, mutta p = 0.051 arvoton?

Vakiotermin (Intercept) merkitys ja merkitsevyys on usein epäselvä ja niin myös tässä tapauksessa. Vakiotermi ilmoittaa suoran arvon pisteessä x = 0 ja niinpä vakio b₀ on mielekäs vain jos piste x = 0 on jotenkin erikoinen ja mielekäs. Origon siirto muuttaa sekä vakiotermiä että sen merkitsevyyttä. Tällä kertaa vakiotermi on ekstrapoloitu sillä yksi toukka oli pienin ruokintamäärä. Vakiotermi ei ole erityisen kiinnostava, sillä se ilmoittaisi jälkeläistuoton silloin kun eläimet eivät saisi lainkaan ruokaa, mikä ei varmaankaan ole kestävä tilanne. Toisaalta on mahdollista, että alimmat ruokintamäärät eivät olleet riittäviä vaan eläimet näkivät jo nyt nälkää ja lisääntyivät vain varantojensa kustannuksella. Mikäli tietäisimme riittävän ravintomäärän kriittisen rajan, voisimme siirtää x-muuttujan nollapisteen tuohon rajaan ja kenties pakottaa vakiotermin nollaksi tuossa pisteessä, mutta meillä ei ole mitään perusteltua syytä pakottaa regressiota nykyisen origon kautta.

Tulostuksen viimeinen rivi kertoo koko mallin merkitsevyyden F-tunnusluvulla eli varianssianalyysillä arvioituna. Koska meillä on vain yhden selittäjän regressio, koko mallin merkitsevyys F-tunnusluvun perusteella on sama kuin ainokaisen selittäjän regressiokertoimen merkitsevyys t-tunnusluvulla arvioituna.

Varianssianalyysi on vielä edessä päin, mutta tässä oppaassa oletetaan, että tilastotieteen perusteet ovat takana päin. Niinpä ei pitäisi olla uutinen, että mallin sopivuutta aineistoon kuvaa jäännösneliösumma (``Residual Sum of Squares'') ja mallin vapausasteiden määrä on havaintojen määrä miinus estimoitujen parametrien määrä. Estimoitujen parametrien lisääntyessä vapausasteet kuluvat ja jäännösneliösummat (todennäköisesti) pienenevät. Jos testattavan ja nollamallin jäännösneliösummat ovat SS₁ ja SS₀ ja vapausasteet df₁ ja df₀, niin F-tunnusluku vapausasteilla df₀ - df₁, df₁ on:

$${\frac{{(SS_0 - SS_1)/(df_0 - df_1)}}{{SS_1 / df_1}}}$$ (5)

R-funktio anova tekee automaattisesti nämä laskut ja arvioi F-tunnusluvun merkitsevyyden kun sille annetaan kaksi luokan lm oliota. Sovitamme ensin mallin, jossa on vain vakiotermi ja sitten vertaamme tätä malliimme lisko.reg, jossa on myös kulmakerroin tekijälle ruoka:

> lisko.null <- lm(munat ~ 1)  ## Vain vakiotermi
> anova(lisko.null, lisko.reg)
Analysis of Variance Table

Model 1: munat ~ 1
Model 2: munat ~ ruoka
  Res.Df Res.Sum Sq Df  Sum Sq F value   Pr(>F)
1      9    174.900                            
2      8     52.473  1 122.427  18.665 0.002545

Mallien jäännösneliösummat olivat SS₀ = 174.9 ja SS₁ = 52.5. Sijoitettuna yhtälöön 5 saamme:

F_{1, 8} = $${\frac{{(174.9-52.473)/(9-8)}}{{52.473/8}}}$$ = $${\frac{{122.427/1}}{{52.473/8}}}$$ = 18.67,

joka on sama tulos kuin suoraan regressionalyysin summaryssa.

Ylläolevan komennon anova tulostuksesta saamme myös muut summaryn alarivien tulokset. Jäännöskeskihajonta, jota R kutsuu nimellä ``Residual standard error'', vastaa muuten tavallista keskihajontaa, mutta se on laskettu regressiosuoran eikä y-muuttujan koko aineiston keskiarvon suhteen. Jäännöskeskihajonta on oletetun mallin teoreettisen virhekeskihajonnan $\sigma$ estimaattori. Se saadaan anova-taulukosta laskuilla

s_resid = $$\sqrt{{\frac{SS_1}{df_1}}}$$ = $$\sqrt{{\frac{52.473}{8}}}$$ = 2.56

``Multiple R-Squared'' on suomeksi determinaatiokerroin eli ``yhteiskorrelaatiokertoimen neliö'', ja se saadaan niin ikään anova-taulukosta:

R² = 1 - $${\frac{{SS_1}}{{SS_0}}}$$ = 1 - $${\frac{{52.437}}{{174.9}}}$$ = 0.70

R² tunnetaan myös nimellä ``selitysaste'': se kuvaa kuinka suuren osan vastemuuttujan arvojen varianssista kokonaisneliösummasta (SS<sub>0</sub>) regressiosuora kykenee ao. aineistossa ``selittämään''. ``Adjusted R-squared'' on samantapainen, mutta sitä on rangaistu (``penalisoitu'') estimoitujen parametrien määrällä.

Näistä alarivien tuloksista suosituin ja yliarvostetuin on determinaatiokerroin eli R². ``Selitysaste'' ei esimerkissämme ole erityisen mielekäs, sillä jälkeläistuotannon varianssi on syntynyt keinotekoisessa kokeessa jossa on systemaattisesti varioitu ruoan määrää. Mikäli ruokintamäärän varianssi olisi ollut vähäisempi (vähemmän ruokintatasoja tai enemmän toistoja keskimmäisillä tasoilla), myös munatuotannon varianssi olisi alentunut ja yhtä hyvällä regressiolla ``selitysaste'' olisi alentunut. Myös havaintoaineistoissa tutkittavan muuttujan vaihtelun kasvattaminen on usein helpoin tapa saada hyviä ``selitysasteita''. Regressiomallin antaman ennusteen hyvyyttä kuvaa paremmin absoluuttinen vaihtelu regressiosuoran ympärillä eli residuaalivaihtelu s<sub>resid</sub>: jos vaihtelu on vähäistä, regressio on ennustusvoimainen riippumatta ``selitysasteesta''. Esimerkissämme R² = 0.7, mikä vaikuttaa korkealta, mutta s_resid = 2.56, mikä ei minusta vaikuta kovinkaan tyydyttävältä.

Kenties vielä paremman kuvan mallimme toimivuudesta saamme tarkastelemalla regressiosuoran ja ennusteiden luottamusrajoja: rajoja, joiden sisällä oikea regressiosuora tai ennustettu yksittäinen arvo on esim. 95 % todennäköisyydellä. Ne voidaan laskea komennolla predict.lm, kuten alla tehdään kuva piirrettäessä.

Luokan lm olioille on oma plot-komento, joka tuottaa lähinnä diagnostiikkaan tarkoitettuja kuvia. Koska meillä on vain yksi selittävä tekijä, tarkastelemme ennemmin perinteisiä regressiokuvia; monimutkaisemmissa malleissa tutkimme oletuskuvia. Peruskuva syntyy komennoilla:

> plot(ruoka,munat,xlab="Ravinto (toukkaa/kerta)",ylab="Munien tuotto")
> abline(lisko.reg)          # Regressiosuora
> abline(h=mean(munat),lty=2)  
> abline(v=mean(ruoka),lty=2)

Viivoja kuvaan lisäävä abline osaa eristää tuloksen regressiokertoimet ja piirtää regressiosuoran. Kaksi viimeistä abline-komentoa lisää viivat muuttujien keskiarvon kohdalle muistuttamaan, että regressiosuora kulkee kummankin muuttujan keskiarvon kautta.

Tarvitsemme vielä luottamusrajat. Ne on laskettava erikseen komennolla predict.lm:

> lisko.ci <- predict(lisko.reg, interval="confidence")

Tuloksena on matriisi, jossa on kolme saraketta: ennustettu arvo (fit), alempi luottamusraja (lwr) ja ylempi raja (upr). Oletuksena on käyttää 95 % rajoja. Selittävät muuttujat ovat valmiiksi oikeassa järjestyksessä, joten voimme lisätä kuvaan rajat komennoilla:

> lines(ruoka,lisko.ci[,"lwr"], lty=2)
> lines(ruoka,lisko.ci[,"upr"], lty=2)

Mikäli ``[,"lwr"]'' näyttää merkilliseltä, parasta kerrata matriisit (§3.6).

Kuva 5: Yksinkertainen regressioanalyysi sekä regression ja ennusteen 95 % luottamusrajat. Pysty- ja vaakasuorat katkoviivat osoittavat kummankin muuttujan keskiarvoja.

Nämä ovat regressiosuoran luottamusrajat: 95 % todennäköisyydellä regressiosuora on rajojen sisällä. Kuten näkyy, luottamusrajat ovat kapeimmillaan selittävän tekijän keskiarvon kohdalla (kuva 5). Mikäli haluamme käyttää regressiosuoraa ennustamaan yksittäisen vesiliskon munatuottoa, on käytettävä vielä laajempia rajoja, eli ennusteen luottamusväliä. Tähän saakka käyttämämme rajat esittävät vain suoran epävarmuutta, mutteivät yksittäisten arvojen poikkeamista regressiosuorasta. Laskemme vielä ennusteen luottamusrajat ja lisäämme ne kuvaan, jolloin saamme lopulliseksi tulokseksi kuvan 5. Paitsi komentoa lines, voimme käyttää myös komentoa matlines, joka lisää molemmat rajat yhtaikaa.

> lisko.pi <- predict(lisko.reg, interval="prediction")
> matlines(ruoka, lisko.pi[,c("lwr","upr")], lty=3, col="black")

Komento matlines vaihtaa väriä ja viivatyyppiä matriisin jokaiselle sarakkeelle, joten jouduimme estämään tämän määritteillä lty ja col. Ennusteen luottamusrajat ovat todella väljät eikä mallimme näytä olevan kovinkaan käyttökelpoinen mikäli tarkoitus on ennustaa yksittäisen vesiliskon tuottama munamäärä. Esim. viidellä toukalla per ruokintakerta, 95 % vesiliskoista tuottaa 0.3...12.7 munaa:

> lisko.pi[5,]
       fit        lwr        upr
 6.4909091  0.2882827 12.6935355

Ranta et al. (1989) esittelevät vain regressiosuoran luottamusrajat, joita he harhaanjohtavasti nimittävät ``ennusteen luottamusrajoiksi''.