Mallin rakentaminen

Perinteisessä regressioanalyysissä mallin rakentaminen tapahtui aivan toisin kuin ANOVAssa. ANOVA oli suunnitellun koeasetelman analyysi, eli kaikki tekijät voitiin ja haluttiin analysoida. Regressionalyysiä sen sijaan käytettiin ennen kaikkea monimuuttujaisen havaintoaineiston analyysiin. Regressioanalyysin selittävät muuttujat olivat ANOVA-termein epätasapainoisia ja useimmiten satunnaistekijöitä. ANOVAssa oli myös tapana näyttää jokaisen koeasetelman tekijän merkitsevyys, myös silloin kun tekijällä ei ollut merkitsevää vaikutusta. Regressioanalyysiin sen sijaan haluttiin sisällyttää vain merkitseviä tekijöitä. Näin ollen regressioanalyysissä rakennettiin mallit siten, että kaikki mukana olevat tekijät olivat merkitseviä.

R:ssä on mallin rakentamista helpottamaan jo luvussa 7.11 nähty funktio drop1 (``pudota yksi'') sekä vastakkaisesti toimiva add1 (``lisää yksi''). Funktio drop1 kokeilee, kuinka mallin hyvyys muuttuisi, jos tekijä poistettaisiin. Tämän jälkeen tekijä voidaan poistaa ja tutkia pelkistettyä mallia. Ilman lisämääreitä kumpikaan funktio ei tulosta muita vertailukriteerejä kuin AIC:n (yhtälö 7 sivulla [*]), mutta pyydämme myös F-tunnusluvun. Katsokaamme aluksi kuinka aineiston plutakko malli yksinkertaistuu:

> attach(plutakko)
> plut.cur <- lm(Fotosynt ~ Pikopl + Plankt + NH4 + PO4)
> drop1(plut.cur, test="F")
Single term deletions

Model:
Fotosynt ~ Pikopl + Plankt + NH4 + PO4
       Df Sum of Sq     RSS     AIC F value     Pr(F)
<none>               8.0973 -8.0841                  
Pikopl  1    3.4485 11.5458 -2.9882  6.3883 0.0232089
Plankt  1   11.8755 19.9728  7.9728 21.9991 0.0002902
NH4     1    0.0606  8.1579 -9.9350  0.1122 0.7422573
PO4     1    0.3043  8.4016 -9.3462  0.5638 0.4643633
> plut.cur <- update(plut.cur, . ~ . -NH4)  ## Poista NH4
> drop1(plut.cur, test="F")
Single term deletions

Model:
Fotosynt ~ Pikopl + Plankt + PO4
       Df Sum of Sq      RSS      AIC F value     Pr(F)
<none>                8.1579  -9.9350                  
Pikopl  1    3.9250  12.0828  -4.0789  7.6980 0.0135335
Plankt  1   12.0106  20.1684   6.1677 23.5563 0.0001761
PO4     1    0.2696   8.4275 -11.2847  0.5288 0.4776216
> plut.cur <- update(plut.cur, . ~ . -PO4)  ## Poista PO4
> drop1(plut.cur, test="F")
Single term deletions

Model:
Fotosynt ~ Pikopl + Plankt
       Df Sum of Sq      RSS      AIC F value    Pr(F)
<none>                8.4275 -11.2847                 
Pikopl  1    4.0185  12.4460  -5.4867  8.1062  0.01114
Plankt  1   13.6820  22.1095   6.0055 27.5994 6.47e-05
Komennolla drop1 näimme aina huonoimman mallitermin, poistimme sen (update) ja jatkoimme kunnes kaikki jäljellä olevat termit olivat merkitseviä.

Komento add1 toimii samoin, mutta tutkii termien lisäyksen vaikutusta, minkä jälkeen voimme aina lisätä kullakin askelella parhaan termin. Lisäyskomento add1 tarvitsee myös tiedon, mitä termejä malliin voi harkita lisättäväksi (``scope''), eli suurimman harkittavan mallin yksipuolinen mallilause:

> plut.cur <- lm(Fotosynt ~ 1)  
> add1(plut.cur, ~ Plankt+Pikopl+PO4+NH4, test="F")
Single term additions

Model:
Fotosynt ~ 1
       Df Sum of Sq     RSS     AIC F value     Pr(F)
<none>              27.0129  8.0116                  
Plankt  1   14.5669 12.4460 -5.4867 21.0674 0.0002272
Pikopl  1    4.9034 22.1095  6.0055  3.9920 0.0610550
PO4     1    1.8888 25.1241  8.5619  1.3532 0.2599121
NH4     1    0.1573 26.8556  9.8949  0.1054 0.7491703
> plut.cur <- update(plut.cur, . ~ . +Plankt) ## Lisää Plankt
> add1(plut.cur, ~ Plankt+Pikopl+PO4+NH4, test="F")
Single term additions

Model:
Fotosynt ~ Plankt
       Df Sum of Sq      RSS      AIC F value   Pr(F)
<none>               12.4460  -5.4867                
Pikopl  1    4.0185   8.4275 -11.2847  8.1062 0.01114
PO4     1    0.3632  12.0828  -4.0789  0.5110 0.48443
NH4     1    0.6712  11.7748  -4.5954  0.9690 0.33873
> plut.cur <- update(plut.cur, . ~ . +Pikopl) ## Lisää Pikopl
> add1(plut.cur, ~ Plankt+Pikopl+PO4+NH4, test="F")
Single term additions

Model:
Fotosynt ~ Plankt + Pikopl
       Df Sum of Sq      RSS      AIC F value  Pr(F)
<none>                8.4275 -11.2847               
PO4     1    0.2696   8.1579  -9.9350  0.5288 0.4776
NH4     1    0.0259   8.4016  -9.3462  0.0493 0.8272
Lisääminen lopetetaan, kun yksikään jäljellä oleva termi ei ole merkitsevä.

Esimerkkimme edustivat kahta perinteistä mallinrakennustapaa: komennolla add1 eteenpäin valitsevaa (``forward selection'') ja komennolla drop1 taaksepäin eliminoivaa (``backward elimination''). Tällä kertaa kumpikin menetelmä johti samaan lopputulokseen, mutta näin ei ole yleensä laita, vaan eri menetelmin päädytään erilaiseen lopulliseen malliin. Yleensä eliminointi (drop1) johtaa suurempiin (useampitekijäisiin) malleihin kuin valinta. Tämä johtuu selittävien tekijöiden välisistä riippuvuuksista (epätasapainoisuudesta: §7.11): jokin tekijä on merkitsevä vain toisen tekijän kanssa ja tulee valituksi malliin vain jos toinen tekijä on jo valittu malliin. Kumpikaan menetelmä ei välttämättä johda parhaan mahdollisen mallin löytämiseen, koska siihen ei välttämättä ole tietä drop1 tai add1 -menetelmillä, vaan meidän pitäisi tutkia useaa mahdollista tekijää yhtaikaa. Havaintoaineistoissa on sangen usein myös samantapaisia tärkeitä tekijöitä, joiden erot voivat olla hyvin pieniä. On enemmänkin sattumaa, kumpi tekijöistä tulee valituksi malliin. Siihen saattaa vaikuttaa jopa valitsemamme tekijänvalintatapa, mutta ainakin otantavaihtelu. Mallinvalinnan perusteella emme siis voi sanoa, että malliin valitsematon tekijä ei ollut tärkeä: malli jossa se olisi ollut mukana saattaa olla melkein yhtä hyvä kuin nykyinen -- tai jopa parempi, jos vaihtoehtoisen tekijän valinta olisi johtanut toiseen tienhaaraan ja toiseen lopulliseen malliin.11

Eteenpäin valitsevassa menetelmässä voi myös käydä, että jokin aiemmalla askeleella malliin valittu tekijä muuttukiin epämerkitseväksi uusien tekijöiden myötä. Askeltavassa valinnassa (``stepwise regression'') kokeillaan jokaisen add1-askelen jälkeen drop1-askelta mallissa oleville termeille. Erilaisilla säännöillä pyritään estämään ikuiset silmukat (eli termi näyttää merkitsevältä ollessaan ulkona mallista, mutta muuttuu epämerkitseväksi mallissa). Askeltavat regressiot ovat monien tilastopakettien suosituimpia ohjelmia.

Paketissa MASS (Venables & Ripley, 1999) on komento stepAIC, joka tekee askeltavan regression AIC:n perusteella. Käyttö on suoraviivaista, mutta tulostus on pitkä, joten näytän vain komennot:

> library(MASS)
> plut.1 <- lm(Fotosynt ~ 1)
> stepAIC(plut.1, ~ Plankt+Pikopl+NH4+PO4)
Mikä johtaa jälleen samaan lopputulokseen. Komento stepAIC toimii sekä lisäävästi että eliminoivasti. Annoimme esimerkissämme lähtökohdaksi mallin, jossa oli vain vakiotermi ja suurimmaksi malliksi kaikki tekijät, joten malli eteni valikoiden eteenpäin. Jos olisimme antaneet komennon
> stepAIC(plut.full)
eli suurimman mallimme (ja jättäneet scopen pois), stepAIC olisi ryhtynyt yksinkertaistamaan mallia.

Olemme kuitenkin tarkastelleet vain yksinkertaista mallia ilman yhdysvaikutuksia. Mikäli määrittelemme suurimman mallin siten, että myös ensimmäisen asteen interaktiot ovat mukana, tuloksemme on toisenlainen:

> plut.iso <- lm(Fotosynt ~ (Plankt+Pikopl+NH4+PO4)^2)
> stepAIC(plut.iso)
Start:  AIC= -3.46
 Fotosynt ~ Plankt + Pikopl + NH4 + PO4 + Plankt:Pikopl + Plankt:NH4 +  
    Plankt:PO4 + Pikopl:NH4 + Pikopl:PO4 + NH4:PO4

                Df Sum of Sq     RSS     AIC
- NH4:PO4        1    0.0668  5.6661 -5.2245
- Plankt:NH4     1    0.0881  5.6874 -5.1496
- Plankt:Pikopl  1    0.0892  5.6886 -5.1455
- Plankt:PO4     1    0.1476  5.7469 -4.9413
- Pikopl:NH4     1    0.4697  6.0690 -3.8507
<none>                        5.5993 -3.4617
- Pikopl:PO4     1    1.2824  6.8817 -1.3372

#### ... Tulostusta poistettu ...
#### Tässä lopullinen malli:

             Df Sum of Sq     RSS     AIC
<none>                     6.3597 -8.9150
- Pikopl:NH4  1    0.6805  7.0402 -8.8818
- Pikopl:PO4  1    1.7242  8.0839 -6.1172
- Plankt      1   13.4084 19.7681 11.7667

Call:
lm(formula = Fotosynt ~ Plankt + Pikopl + NH4 + PO4 + Pikopl:NH4 + Pikopl:PO4)

Coefficients:
(Intercept)       Plankt       Pikopl          NH4          PO4   Pikopl.NH4  
   -0.20783      0.59332      0.75285     -0.03788      0.31123      0.21538  
 Pikopl.PO4  
   -0.82396  
Lopullinen malli on selvästi isompi. Tämä johtuu paljolti siitä, että kriteerinä on AIC eikä F -- AIC:n käyttöhän merkitsee samaa kuin ``epämerkitsevän'' rajan F = 2 käyttö. Mukana ovat myös muuttujan Pikopl ja ravinnemuuttujien interaktiot. Huomattakoon, että sekä stepAIC että drop1 ymmärtävät sen verran mallinrakennuksen päälle, että ne yrittävät poistaa vain tekijöitä, joille poistaminen on sallittu: sellaista tekijää, joka on jonkin muun tekijän osatekijä ei saa ajatellakaan poistettavaksi. Ensimmäisessä askelessa harkittiin vain interaktiotermien poistamista: päävaikutusta ei voi harkita poistettavaksi ennen kuin on poistettu kaikki interaktiot, joihin se osallinen. Viimeisessä vaiheessa harkittiin poistettavaksi vain termejä Pikopl:NH4, Pikopl:PO4 ja Plankt. Muut jätettiin rauhaan, sillä ne sisältyivät mallissa vielä oleviin interaktioihin.

Sama sääntö pätee myös polynomisiin ja muihin vastaaviin regressioihin. Jos meillä on regressiomalli

E(y) = b0 + b1x + b2z + b3x2 + b4z2 + b5xz

Niin ensimmäisellä askeleella saamme harkita vain termien x2, z2 ja xz poistamista. Saamme harkita termin x poistamista vasta kun sekä x2 että xz on poistettu. Loogisesti edeten, saamme harkita vakiotermin poistamista vasta kun kaikki muut termit on poistettu. Tästä yleissäännöstä saa poiketa vain hyvin perustellusta syystä silloin kun syntyneellä mallillamme on selvä tulkinta. Esimerkiksi malli E(y) = b1x määrittelee regression origon kautta ja voi olla mielekäs kun origolla on absoluuttinen ja mielekäs tulkinta. Tällainen malli voi olla mielekäs, vaikka useimmat näkemäni origon kautta pakotetut regressiot ovat mielettömiä ja virheelliseen päättelyyn pohjautuvia. Malli E(y) = b0 + b2x2 taas määrittelee paraabelin jonka kulmakerroin on 0 pisteessä x = 0. Tällainenkin malli voi olla mielekäs, vaikkakin harvoin. Tällainen mallin rakentaminen on hankalampaa, sillä R ei automaattisesti osaa pitää huolta polynomisista riippuvuuksista (ainakaan tätä kirjoitettaessa).

Olen viipynyt mallin rakentamisessa sangen pitkään, sillä aihe on hankala. Vielä tärkeämpi syy on, että sovitettaessa epätasapainoisia, laajoja lm-malleja, normaali käytäntö on regressioanalyysin perinteestä kumpuava pikemmin kuin ANOVA-taulukoista. Samoja periaatteita käytetään sekä jatkuville että luokkamuuttujille selittäjinä. Daphnia-esimerkissä (§7.10) sovelsimme jo itse asiassa näitä periaatteita: interaktiotermin poistaminen perustui mallien vertailuun ja päätekijän (lämpötila) tunnusluvut arvioitiin tästä pelkistetystä mallista.

Mallin rakentaminen on paitsi vaikeaa myös vastuullista. Tutkija ei voi piiloutua automaattisten proseduurien taakse vaan hänen on otettava vastuu mallistaan. Tämä merkitsee tietystikin tässä luvussa esitettyjen muodollisten sääntöjen noudattamista. Ennen kaikkea se merkitsee kuitenkin, että tutkija pyrkii rakentamaan mielekkään mallin. Mitä tahansa tekijöitä ei sovi heittää selittäjäksi vain satunnaisen korrelaation vuoksi, vaan lopullisen yhtälön on oltava myös ekologisesti järkevä. Tämä edellyttää yleensä pitkällistä työskentelyä ja harkintaa. On myös muistettava, että tilastollisilla perusteilla meillä voi olla monta suunnilleen yhtä hyvää vaihtoehtoista mallia.


Jari Oksanen 2003-01-21