Monen selittävän tekijän regressio

Yhden selittävän tekijän regressioyhtälö (4) on helposti laajennattavissa monen selittävän tekijän malliksi:

E(y) = b0 + b1x1 + b2x2 + ...   + bpxp (6)

Tämäkin on lineaarinen malli, vaikka tulokset saattavat joskus olla käyriä, kuten

E(y) = b0 + b1x + b2x2

joka on paraabeli. Lineaarisuus regressiomallissa tarkoittaa, että malliyhtälössä on vakiotermin lisäksi vain kerroin x tekijä -tyyppisiä termejä. Hienommin sanottuna, malli on lineaarinen parametrien suhteen (linear in parameters).

Yhtälöä 6 vastaa R:ssä mallilauseke (vakiotermin ``1'' voi tietysti jättää kirjoittamatta näkyviin):

lm(y ~ 1 + x1 + x2 + ... + xp)

Monen selittävän muuttujan regression sovittaminen R:ssä on aivan yhtä helppoa kuin mallilausekkeen kirjoittaminen. Tässä suhteessa ei paljonkaan uutta. Monen selittäjän malleihin liitty kuitenkin ongelmia tulkinnassa, esitystekniikassa ja mallin rakentamisessa, joten niitä on tarkasteltava erikseen.

Biometrian kirjassa käsitelty esimerkki (Esimerkki 10.8 Ranta et al., 1989, aineisto plutakko paketissa rekola) tutkii kalliolammikoiden perustuotantoa, oletettavasti Tvärminnessä. Vastemuuttujana on hiilen yhteytys $ \mu$g dm-3 h-1 (muuttuja Fotosynt), mahdollisina selittäjinä pikoplanktonin (muuttuja Pikopl) ja suuremman planktonin (muuttuja Plankt) klorofylli-a sekä levien käyttämättä jättämä NH4 (muuttuja NH4) ja PO4 (muuttuja PO4), kaikki selittäjät $ \mu$g dm-3.

Biometrian kirjassa tarkastellaan kahta mallia (Ranta et al., 1989). Ensimmäisessä selittäjinä ovat planktonfraktioiden klorofyllitiheydet:

> data(plutakko)
> plut.lm <- lm(Fotosynt ~ Pikopl + Plankt, data=plutakko)
> summary(plut.lm)

Call:
lm(formula = Fotosynt ~ Pikopl + Plankt, data = plutakko)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.18388    0.26309   0.699   0.4941
Pikopl       0.50003    0.17562   2.847   0.0111
Plankt       0.45827    0.08723   5.254 6.47e-05

Residual standard error: 0.7041 on 17 degrees of freedom
Multiple R-Squared: 0.688,      Adjusted R-squared: 0.6513
F-statistic: 18.75 on 2 and 17 degrees of freedom,      p-value: 5.013e-05

Molemmat tekijät ovat hyvin merkitseviä ja niiden kertoimet keskenään melkein yhtäsuuria.

Toisessa mallissa käytettiin kaikkia selittäviä tekijöitä, joten R:n mallilausekkeen voi mukavasti lyhentää:

> plut.full <- lm(Fotosynt ~ . , data=plutakko)
> summary(plut.full)

Call:
lm(formula = Fotosynt ~ ., data = plutakko)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.15390    0.65204   0.236  0.81660
Pikopl       0.52813    0.20895   2.528  0.02321
Plankt       0.49325    0.10516   4.690  0.00029
NH4          0.03771    0.11255   0.335  0.74226
PO4         -0.15870    0.21136  -0.751  0.46436

Residual standard error: 0.7347 on 15 degrees of freedom
Multiple R-Squared: 0.7002,     Adjusted R-squared: 0.6203
F-statistic:  8.76 on 4 and 15 degrees of freedom,      p-value: 0.0007443

Ravinteet eivät ole merkitseviä, mutta planktonfraktiot ovat. Planktonfraktioiden kertoimet muuttuivat ja kertoimien p-arvot jopa huononivat. Myös anova on sitä mieltä, että ravinnetermien lisääminen ei merkitsevästi parantanut mallia:

> anova(plut.lm,plut.full)
Analysis of Variance Table

Model 1: Fotosynt ~ Pikopl + Plankt
Model 2: Fotosynt ~ Pikopl + Plankt + NH4 + PO4
  Res.Df Res.Sum Sq Df Sum Sq F value Pr(>F)
1     17     8.4275                         
2     15     8.0973  2 0.3302  0.3058  0.741

anova testaa nyt kumpaakin ravinnetekijää yhtaikaa, joten vapaustasteiden muutos on 2. Tilastolliselta kannalta meillä ei ole mitään syytä pitää mallissa mukana ravinnetekijöitä. Meillä ei ole siihen mitään syytä myöskään ekologiselta kannalta: ravinnepitoisuudet on mitattu vedestä eli ne ovat kasvien käyttämättä jättämiä ravinteita, joiden ei voi odottaakaan selittävän fotosynteesiä.

Tarkastelemme perusteellisemmin ensimmäistä, vain planktonfraktiot sisältävää mallia. Koska meillä on kaksi selittävää tekijää, emme voi esittää tuloksia yksinkertaisena regressiokuvaajana, jossa x-akselilla on selittävä tekijä. Sen sijaan voimme tutkailla luokan lm olion oletuskuvaa (kuva 6):

> par(mfrow=c(2,2))       ## 4 osakuvaa
> plot(plut.lm)

Kuva 6: Lineaarisen mallin tuloksen diagnostinen oletuskuva (aineisto plutakko)

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

Komento plot.lm tuottaa neljä diagnostista kuvaa, joita voi käyttää mallin arviointiin:

  1. Residuals vs Fitted: Regression residuaalit vastaan ennustetut arvot. Mikäli malli sopii, residuaalit muodostavat tasaisen vyön suoran Residuals = 0 ympärille. Näin ei käy vaan residuaalit tuntuvat taipuvan alaspäin korkeammissa sovitetuissa arvoissa. Tämä viittaa siihen, että yksinkertainen lineaarinen mallimme ei sovi, vaan meidän pitäisi tutkia muita muotoja. Lisäksi hajonta saattaa hieman kasvaa kohti korkeampia ennustettuja arvoja. Meillä näkyy myös olevan muutama hyvin poikkeava arvo; havainnot 2, 3 ja 5 on merkitty erityisen poikkeaviksi.

  2. Normal Q-Q plot: Vaaka-akselilla on normaalijakauman mukaiset teoreettiset kvantiilit, pystyakselilla standardisoidut residuaalit. Mikäli residuaalit ovat normaalisesti jakautuneita, kuvassa on suora viiva. Mitä ilmeisimmin normaalisuusoletusta vastaan on rikottu. Etenkin havainnot 2, 3 ja 5 käyttäytyvät huonosti.

  3. Scale-Location plot: Kuten osakuva 1, mutta nyt näytetään vain poikkeaman suuruus, ei suuntaa. Tämä kuva osoittaa selkeämmin kuin osakuva 1 varianssin kasvun suurissa ennustetuissa arvoissa. Etenkin havainnot 2, 3 ja 5 käyttäytyvät huonosti.

  4. Cook's distance plot: Suuri Cookin etäisyys tarkoittaa, että yksittäisellä pisteellä on suuri vaikutus regressiosuoraan. Etenkin pisteet 1, 3 ja 5 ovat huolestuttavia.

Näissä kuvissa voisi raakaresiduaalien sijaan käyttää ulkoisesti studentoituja (externally studentized) residuaaleja (katso help(rstudent)), jotka ovat mm. herkempiä paljastamaan oudokkeja (outlier).

Kaikki diagnostiset kuvat osoittivat mallimme jossain määrin huonosti sopivaksi. Havainnot 3 ja 5 olivat poikkeavia kaikissa osakuvissa ja 2 useimmissa. Sangen tavallinen ratkaisu on poistaa poikkeavat havainnot aineistosta ja koettaa uudelleen. Meillä saattoi kuitenkin olla muitakin huolen aiheita kuin kolme outoa havaintoa (mikä on 12 % aineistostamme), joten tutkimme mallia hieman paremmin.

Planktonfraktioiden kertoimien erot olivat hyvin pienet verrattuna keskivirheisiinsä, joten voisimme kokeilla näiden kertoimien pakottamista identtiseksi. Käytännössä tämä tapahtuu laskemalla yhteen muuttujat Pikopl ja Plankt. Tämä on itse asiassa myös biologisesti mielekästä: klorofyllihän yhteyttää, eikä ratkaisevaa välttämättä ole minkäkokoisen levän sisällä yhteyttävä rakenne on (itse asiassa planktonin jako kokofraktioihin on tyypillistä patriarkaatin ekologiaa, missä luonto pakotetaan tutkimuslaitteen asettamiin rajoihin). Ensimmäinen osakuva (kuva 6) viittasi myös käyräviivaiseen suhteeseen. Samaan viittaa muuttujien sirontakuvioiden tarkastelu:

> pairs(plutakko, panel=panel.smooth)

Tämän komennon kuvat (joita en näytä) viittaavat siihen, että voimme mahdollisesti aproksimoida fotosynteesin ja klorofyllitiheyden suhdetta toisen asteen yhtälöllä. On aivan varmaa, että paraabeli ei ole oikea suhde, mutta se on sopiva alustavaksi ratkaisuksi.

R:ssä on funktio poly, jonka avulla saamme nopeasti toisen asteen polynomin. Sen sijaan yhteenlaskun sijoittaminen mallilausekkeeseen vaatii pienen tempun, sillä mallilauseessa ``+'' tarkoittaa yleensä kummankin tekijän sisällyttämistä erikseen malliin. R:ssä on funktio I, joka eristää (``isolates'') operaattorit, joita käytetään tavanomaisessa matemaattisessa merkityksessä eikä mallitermeinä. Lopullinen mallimme on:

> plut.new <- lm(Fotosynt ~ poly(I(Pikopl+Plankt),2), data=plutakko)
> summary(plut.new)

Call:
lm(formula = Fotosynt ~ poly(I(Pikopl + Plankt), 2), data = plutakko)

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)
(Intercept)                    1.3465     0.1407   9.567 2.95e-08
poly(I(Pikopl + Plankt), 2)1   4.3086     0.6295   6.845 2.85e-06
poly(I(Pikopl + Plankt), 2)2  -1.3090     0.6295  -2.080    0.053

Residual standard error: 0.6295 on 17 degrees of freedom
Multiple R-Squared: 0.7507,     Adjusted R-squared: 0.7213
F-statistic: 25.59 on 2 and 17 degrees of freedom,      p-value: 7.461e-06

Myös diagnostinen kuva (plot(plut.new)) on paljon parempi - sitä ei kuitenkaan näytetä, vaan se on tuotettava itse R-istunnossa. Aineistossa on yhä vielä yksi häirikkö (havainto 5), mutta muuten kaikki diagnostiset kuvaajat ovat paljon parempia. Toisen asteen termi ei itse asiassa ole merkitsevä perinteisten rajojen mukaan, mutta sen mukana olo tuntuu tekevän regressiostamme realistisemman mallisen.

Nyt meillä onkin itse asiassa vain yksi selittävä muuttuja, joten voimme piirtää perinteisen regressiokuvan (kuva 7):

> plut.ci <- predict(plut.new,interval="confidence") # 1
> plut.pi <- predict(plut.new,interval="prediction") # 2
> i <- sort.list(Pikopl+Plankt)                      # 3
> par(mfrow=c(1,1))                                  # 4
> plot(Pikopl+Plankt, Fotosynt,xlab="Klorofyllitiheys",ylab="Yhteytys")
> matlines((Pikopl+Plankt)[i], plut.ci[i,], lty=c(1,2,2), col="black")
> matlines((Pikopl+Plankt)[i], plut.pi[i,], lty=c(1,3,3), col="black")
> identify(Pikopl+Plankt, Fotosynt)                  # 8
[1]  1  2  3  5  8 10 11

Kuva 7: Lopullinen malli: regressiosuora, sen luottamusrajat sekä ennusteen luottamusrajat (aineisto plutakko)

\includegraphics[width=0.6\textwidth]{plutpoly.eps}

Komentojono voi näyttää kauhistuttavalta, muttei kuitenkaan paljon poikkea aikaisemmista -- nyt vain on hieman enemmän piirtämistä. Lisäksi kuva syntyi ``iteratiivisesti'', vaikka välivaiheet onkin yltä poistettu. Meitä kiinnostaa nyt sekä regressiokäyrän (rivi 1) että ennusteiden luottamusvälit (rivi 2). Käytämme ennustamiseen vanhaa aineistoa.9Sen takia meidän on järjestettävä kaikki ennusteet x-muuttujan mukaan ja rivillä 3 laskemme tähän tarvitsemamme indeksin i, jota käytämme myöhemmin riveillä 6 ja 7. Piirrämme vain yhden kuvan, joten poistamme edellisen osakuvamäärittelyn (rivi 4) ja piirrämme peruskuvan akseleineen ja datapisteineen (rivi 5). Sen jälkeen lisäämme sovitetut arvot sekä luottamusrajat (rivit 6 ja 7): tällä kertaa matlines helpottaa huomattavasti viivojen piirtoa. Komento matlines vaihtaa viivatyyppiä ja viivan väriä jokaiselle sarakkeelle, joten määritämme värin ja viivatyypin kullekin sarakkeelle komennon määreissä. Lopulta identifioimme joitain pisteitä: rivi 8 käynnistää interaktiivisen komennon identify, joka lisää havaintonumerot niihin pisteisiin, joita näpäytämme hiiren vasemmalla näppäimellä.

Vaikka teimme periaatteessa lineaarisen regression, saimme tulokseksi käyrän (kuva 7). Mallimme on silti lineaarinen, sillä se on lineaarisen regressioyhtälön erikoistapaus (yhtälö 6). Polynomimalli ei varmastikaan ole tässä tapauksessa oikea muoto: ekstrapoloitu käyrä painuisi alas pian tutkitun alueen ulkopuolella. Ilmeisestikin oikea muoto olisi jonkinlainen saturaatiokäyrä, missä fotosynteesi nousee aluksi melkein lineaarisesti, mutta korkeammassa klorofyllitiheydessä tasoittuu melko vakaaksi. Tällaiselle käyrälle olisi helppo keksiä ekologisesti mielekkäitä selityksiä, mutta pidättäydyn tästä, sillä meillä ei ole mitään perusteita tietää, onko selitys oikea. Vaikka regression luottamusvälit ovat melko kapeat, yksittäisen kalliolammikon yhteytyksen ennustaminen näyttää hyvin epävarmalta (kuva 7).


Jari Oksanen 2003-01-21