R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Tilastolliset mallit R:ssä

Johdanto

R:ään on toteutettu varsin ilmaisuvoimainen tilastollisten mallien määrittelytapa. R:n manuaali ”An Introduction to R” kattaa suuren osan erilaisista määrittelyistä. Alkuperäinen idea lienee kotoisin Wilkinsonin ja Rogersin artikkelista.

Seuraavassa käydään läpi muutamia erilaisia mallivaihtoehtoja eräälle esimerkkiaineistolle. Esimerkkiaineiston löytää tämän postauksen lopusta R-koodimuodossa.

Kaavojen perusteet

R:n kaavoissa on kolme keskeistä osaa, niin sanottu vasen puoli (vastemuuttuja(t)), oikea puoli (selittävät muuttujat) ja näitä erottava tilde (~). Mallissa on aina oletusarvoisesti mukana ”y-akselin leikkauspiste” eli Intercept. Intercept:iä merkitään mallissa numerolle 1. Se voidaan eksplisiittisesti lisätä malliin käyttämällä ilmaisua ”+1” tai poistaa se mallista ilmaisulla ”-1”. Muuttujien yhteisvaikutuksia kuvataan kahdella eri merkinnällä: * ja :. Hierarkkisia ja muita vastaavia rakenteita (nesting) kuvataan merkillä /.

Lineaarinen regressio, ANOVA ja ANCOVA

R:ssä lineaarinen regression, varianssianalyysi (ANOVA) ja covarianssianalyysi (ANCOVA) voidaan kaikki sovitta samalla komennolla lm(). Näiden eronahan on se, että lineaarisessa regressiossa on ainostaan jatkuvia selittäjiä, ANOVA:ssa vain kategorisia selittäjiä ja ANCOVA:ssa sekä jatkuvia että kategorisia selittäjiä.

ANOVA
Esimerkkiaineistoon voidaan sijoittaa yksinkertainen ANOVA esimerkiksi seuraavasti:

fit.anova<-lm(Weight~Sex, data=d)

Komennossa ”Weight~Sex” siis määrittelee mallin, jossa painoa (Weight) ennustetaan sukupuolella (Sex). Argumentti data määrittelee sen data framen nimen, josta kyseiset mallissa esiintyvät muuttujat löytyvät (saman nimiset sarakkeet!).

Komennolla summary() saadaan mallin tulokset luettavassa muodossa. Ennen niiden tulkitsemista on syytä myös tarkastella mallin oletusten toteutumista graafisesti (plot(fit.anova)), mutta hypätään tässä esityksessä ko. vaiheen yli.

# Vaihtoehto 1
summary(lm(Weight~Sex, data=d))
# Vaihtoehto 2
summary(fit.anova)
 
Call:
lm(formula = Weight ~ Sex, data = d)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-2.1385 -0.9463  0.1057  0.7799  1.8856 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   9.0703     0.1984  45.722   <2e-16 ***
Sexmale      -0.8316     0.2806  -2.964   0.0044 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.087 on 58 degrees of freedom
Multiple R-squared:  0.1316,    Adjusted R-squared:  0.1166 
F-statistic: 8.786 on 1 and 58 DF,  p-value: 0.004398

Tuloksen Coefficients-taulukosta nähdään, että Intercept:in arvoksi on arvioitu 9.0703 ja sukupuolen vaikutukseksi -0.8316. Sukupuoli esiintyy tuloksessa muodossa Sexmale, mikä pitää tulkita siten, että uros-sukupuolta on verrattu naaras-sukupuoleen, ja urosten paino on pienempi kuin naarasten, koska vaikutus on negatiivinen.

Miten Intercept:in arvo sitten määräytyy? Tässä tapauksessa se on naarasten keskiarvo, ja urosten vaikutus on naarasten keskiarvon ja urosten keskiarvon erotus. Tämän nähdään vaikkapa laskemalla ryhmäkohtaiset keskiarvot:

aggregate(d$Weight, list(d$Sex), mean)
  Group.1        x
1  female 9.070312
2    male 8.238705

Kaksisuuntaisessa ANOVA Intercept on vasteen estimoitu arvo, kun kaikki selittäjät saavat arvokseen verrokkitason. Seuraava malli antaa Intercept:ksi arvon 8.8358. Tuo on hyvin lähellä painon havaintoarvojen keskiarvoa, kun keskiarvo lasketaan niistä tapauksista, joissa Sex=female ja Genotype=CloneA. Arvo ei kuitenkaan ole tarkalleen sama, sillä ennustettu arvo eroaa tästä hieman. Ennustettu arvo saadaan laskettu komennolla predict():

fit.anova<-lm(Weight~Sex+Genotype, data=d)
coef(summary(fit.anova))
 
                 Estimate Std. Error   t value     Pr(>|t|)
(Intercept)     8.8357555  0.1725066 51.219807 8.032504e-47
Sexmale        -0.8316062  0.1304027 -6.377214 4.520596e-08
GenotypeCloneB  0.9677839  0.2258642  4.284805 7.743319e-05
GenotypeCloneC -1.0436077  0.2258642 -4.620510 2.485289e-05
GenotypeCloneD  0.8239635  0.2258642  3.648049 6.039358e-04
GenotypeCloneE -0.8753990  0.2258642 -3.875776 2.947099e-04
GenotypeCloneF  1.5345955  0.2258642  6.794329 9.664664e-09
 
mean(d$Weight[d$Genotype=="CloneA" & d$Sex=="female"]) # 8.849022
newd<-data.frame(Sex="female", Age=1, Genotype="CloneA", Score=4)
predict(fit.anova, newd) # Sama kuin mallin antama Intercept:in arvo
 
       1 
8.835756

Lineaarinen regressio

Lineaarinen regressio määritetään samaan tyyliin kuin ANOVA:kin, mutta selittäjänä on jatkuva muuttuja, tässä ikä:

fit.lm<-lm(Weight~Age, data=d)
summary(fit.lm)
 
Call:
lm(formula = Weight ~ Age, data = d)
 
Residuals:
     Min       1Q   Median       3Q      Max 
-1.95515 -0.90267  0.05664  0.93925  2.14454 
 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.7558     0.3280  23.644  < 2e-16 ***
Age           0.2996     0.0989   3.029  0.00366 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 1.083 on 58 degrees of freedom
Multiple R-squared:  0.1366,    Adjusted R-squared:  0.1217 
F-statistic: 9.175 on 1 and 58 DF,  p-value: 0.003659

Tuloksesta nähdään, että iän (Age) vaikutus on 0.2996 eli toisin sanoen paino lisääntyy 0.2996 yksikköä kutakin vuotta kohden. Intercept:n arvo on mallilla estimoitu vasteen keskiarvo, kun selittävän muuttujan arvo on 0:

newd<-data.frame(Sex="female", Age=0, Genotype="CloneA", Score=4) # Selittäjä (Age) = 0
predict(fit.lm, newd) # Mallin intercept
 
      1 
7.75576

ANCOVA

ANCOVA-mallissa on sekä jatkuvia että kategorisia muuttujia. Erilaisia mallivaihoehtoja on useammalla muuttujall useita erilaisia:

# Analysis of covariance, just main effects
fit.ancova.noint<-lm(Weight~Age+Sex, data=d)
# Analysis of covariance, interactions
# Different intercepts and slopes for every combination of Age and Genotype
fit.ancova.int<-lm(Weight~Age*Sex, data=d)
fit.ancova.int<-lm(Weight~Age+Sex+Age:Sex, data=d) # Sama kuin edellinen
# Separate regression models of type ~1+Age in each Genotype level
fit.ancova.compslope<-lm(Weight~Genotype/Age-1, data=d)
fit.ancova.compslope<-lm(Weight~Genotype+Genotype:Age-1, data=d) # Sama kuin edellinen

Mallin tulkinta on samatapainen kuin aiemmissa esimerkeissäkin. Esimerkiksi pelkkien päävaikutusten mallin tulos näyttää seuraavalta:

Call:
lm(formula = Weight ~ Age + Genotype, data = d)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-0.7704 -0.4364  0.0027  0.4043  0.9081 
 
Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)     7.52120    0.20816  36.133  < 2e-16 ***
Age             0.29958    0.04542   6.595 2.02e-08 ***
GenotypeCloneB  0.96778    0.22253   4.349 6.25e-05 ***
GenotypeCloneC -1.04361    0.22253  -4.690 1.96e-05 ***
GenotypeCloneD  0.82396    0.22253   3.703 0.000509 ***
GenotypeCloneE -0.87540    0.22253  -3.934 0.000245 ***
GenotypeCloneF  1.53460    0.22253   6.896 6.63e-09 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
Residual standard error: 0.4976 on 53 degrees of freedom
Multiple R-squared:  0.8336,    Adjusted R-squared:  0.8147 
F-statistic: 44.24 on 6 and 53 DF,  p-value: < 2.2e-16

Interceptin tulkinta on samanlainen kuin aiemmissakin tapauksissa eli se on ennustettu keskiarvo, kun AGE=0 ja Genotype=CloneA:

newd<-data.frame(Sex="female", Age=0, Genotype="CloneA", Score=4)
predict(fit.ancova.noint, newd)
 
       1 
7.521204

ANCOVA:n avulla on mahdollista siis sovittaa malli, jossa kullekin kategorisen muuttujan tasolle on sovitettu oma lineaarinen regressiomalli. ANCOVA voi olla helpointa hahmottaa käyttämällä erilaisia graafisia esityksiä. Eräs näppärimmistä plot-komennoista tätä tarkoitusta varten on coplot():

# From Zuur, Ieno: A Beginners Guide to R
panel.lm<-function(x, y, ...) {
   tmp<-lm(y~x, na.action=na.omit)
   abline(tmp, col="#CC0000", lwd=2)
   points(x, y, ...)
}
# Palkki (|) määrittelee kategorisen muuttujan
coplot(formula=Weight~Age|Genotype, data=d, panel=panel.lm, pch=16, cex=0.75)

Tuloksena on seuraavanlainen kuva, johon on sovitettu samaan tapaan lineaarinen malli kuin ylläkin on tehty. Kuhunkin genotyyppikategoriaan on siis sovitettu erillinen lineaarinen regressiomalli. Kuten edellä olevasta mallin tuloksesta voi päätellä, on mallissa vain yksi intercept, ja vaikutukset on laskettu ikäänkuin poikeamina siitä.

Sekamallit

Edellä esitetyissä malleissa on mukana ainoastaan ns. kiinteitä tekijöitä (fixed factors). Tällaisille muuttujille arvioidaan aina vaikutuksen suuruus ja esimerkiksi tilastollinen merkitsevyys. Sekamalleissa on samassa mallissa sekä kiinteitä tekijöitä että satunnaistekijöitä (random factors). Sekamalleissa satunnaistekijöiden eri tasojen ajatellaan usein olevan vain otos suuremmasta määrästä erilaisia mahdollisia tasoja, eikä satunnaistekijöille arvioida esimerkiksi tilastollista merkitsevyyttä, vaan ne vaikuttavat vain malliin sitä kautta, että ne vaikuttavat varianssin arviointiin.

Eräs esimerkki usein käytetystä sekamallien sovelluksesta on surveytutkimus, jossa otos muodostaa hierarkkisen rakennelman. Esimerkiksi survey, jossa ensin valitaan satunnaisesti kouluja, joista valitaan luokkia, ja luokista oppilaita, muodostaa hierarkkisen asetelman, jossa varsinainen havaintoyksikkö on oppilas, joka sitten kuuluu niin yhteen luokkaan kuin kouluunkin.

Toinen yleinen sekamallien sovellus on aikasarja-analyysi, jossa havaintoykisköitä on ajan kuluessa otettu toistomittauksia. Tällöin tavanomaisen lineaarisen regression oletukset havaintojen riippumattomuudesta rikkoutuvat, mutta tilanne voidaan silti analysoida esimerkiksi sekamallilla.

Sovitetaan ensin aineistoon maksimaalinen lineaarinen (ANCOVA) malli (fit1) ja hieman yksinkertaistettu malli (fit2):

fit1<-lm(Weight~Age*Sex*Genotype, data=d)
fit2<-lm(Weight~Age+Sex*Genotype, data=d)

Malli voidaan sovittaa myös sekamallina riippuen siitä, millainen koeasetelma oikeasti on ollut. R:n paketti nlme mahdollistaa sekamallien sovittamisen. Muitakin paketteja on, kuten esimerkiksi lme4, joka on saman tekijän kuin nlme:kin.

Jos oletetaan, että genotyyppiä halutaan käsitellä satunnaistekijänä, on erilaisia mahdollisia kaavamalleja useita. Käydään läpi muutamia yleisimpiä:

Random Intercept
Sovitetaan kullekin genotyypille malli, jossa sovitettavan regressiosuoran kulmakerroin säilyy samana kaikilla ikä- ja sukupuoli-yhdistelmillä, mutta Intercept saa vaihdella eri genotyyppien välillä:

library(nlme)
fit.lme1<-lme(fixed=Weight~Age+Sex, random=~1|Genotype, data=d)

Random-argumentti spesifioi satunnaistermien mallin, tässä siis Intercept (1) saa vaihdella eri genotyypin tasojen välillä (|Genotype).

Random Intercept ja slope

Ensimmäisessä mallissa iän vaikutukselle arvioidaan sekä random Intercept että slope. Jälkimmäisessä mallissa vastaava sovitetaan sukupuolen suhteen.

library(nlme)
fit.lme2<-lme(fixed=Weight~Age+Sex, random=~1+Age|Genotype, data=d)
# tai
fit.lme3<-lme(fixed=Weight~Age+Sex, random=~1+Sex|Genotype, data=d)

Useita satunnaisvaikutuksia

Seuraavassa sekä iän että sukupuolen suhteen sovitetaan satunnaisvaikutus (random Intercept ja slope):

library(nlme)
fit.lme4<-lme(fixed=Weight~Age+Sex, random=list(Age=~Genotype, Sex=~Genotype), data=d)

Repeated measures ANOVA

library(nlme)
fit.lme5<-lme(fixed=Weight~Sex, random=~1|Sex/Genotype, data=d)
# Usein myös yllä esitettyä random slope-mallia kutsutaan
# toistomittaus-ANOVA:ksi

Tyypillinen toistomittaus-ANOVA vaatii pienen lisäyksen dataan. Luodaan dataan kutakin yksilöä kuvaava muuttuja id:

d$id<-rep(1:12, each=5)

Tämän jälkeen malli voidaan sovittaa esimerkiksi seuraavasti, jos oletetaan että ajan suhteen tehdyt mittaukset ovat toistoja samasta yksilöstä:

fit.lme5<-lme(fixed=Weight~Sex*Genotype*Age, random=~1|id/Genotype, data=d)
# tai (ei vastaava kuin yllä oleva malli, vaan yksinkertaisempi)
fit.lme5<-lme(fixed=Weight~Age*Sex, random=~1|id/Genotype, data=d)

Tiivistelmä

Edellä on käsitelty useita erilaisia tapoja määrittää tilastollisia malleja. Tiivistelmänä nämä voidaan esittää listanomaisesti:

# Yksisuuntainen ANOVA
fit.anova<-lm(Weight~Sex, data=d)
# Kaksisuuntainen ANOVA
fit.anova<-lm(Weight~Sex+Genotype, data=d)
# Lineaarinen regressio
fit.lm<-lm(Weight~Age, data=d)
# Kovarianssianalyysi (ANCOVA)
fit.ancova.int<-lm(Weight~Age*Sex, data=d)
 
# Sekamalli, random intercept
fit.lme1<-lme(fixed=Weight~Age+Sex, random=~1|Genotype, data=d)
# Sekamalli, random intercept ja slope
fit.lme2<-lme(fixed=Weight~Age+Sex, random=~1+Age|Genotype, data=d)
# Sekamalli, repeated measures ANOVA
fit.lme5<-lme(fixed=Weight~Age*Sex, random=~1|id/Genotype, data=d)

Mallien piirteet ja sovittaminen on tässä kuvattu pääpiirteissään mutta suurpiirteisesti. Jossakin vaiheessa sekamallien sovittaminen sekä nlme- että lme4-pakettia käyttäen on tarkoitus kattaa nykyistä tarkemmalla tasolla.

Esimerkiaineisto

# http://www.cas.miamioh.edu/~schaefrl/data/Crawley%20Data/Gain.txt
# Muuttujat:
#   vaste:      Weight
#   selittäjät: Sex, Age, Genotype
d<-structure(list(Weight = c(7.44563039, 8.000223376, 7.705105123, 
8.348736726, 8.454716107, 8.778212882, 8.340044871, 9.092006302, 
9.105216761, 9.845136836, 6.100193712, 6.513768952, 7.236199692, 
7.228504191, 7.650989006, 7.874350905, 8.511566673, 8.607624762, 
8.994112473, 9.689926545, 6.479922339, 6.828618805, 7.376614806, 
7.245332705, 7.923819602, 9.143637003, 9.170335854, 9.379750168, 
10.03772954, 10.0531328, 8.128602512, 8.516594654, 8.987419296, 
8.956231029, 9.656264805, 9.18797551, 9.618807211, 9.816943243, 
9.71091241, 10.38210711, 7.685253334, 7.384318914, 7.907385712, 
7.946400479, 8.110432545, 9.06934822, 9.580191173, 9.568835855, 
10.04826226, 10.4949405, 7.52193091, 7.394792392, 7.767411841, 
8.263021356, 8.644069167, 9.714612618, 9.721849885, 10.79904741, 
10.95591672, 10.56946689), Sex = structure(c(2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("female", "male"), class = "factor"), 
    Age = c(1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
    4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
    4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
    4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 
    4L, 5L), Genotype = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 
    5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
    2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 
    5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L), .Label = c("CloneA", "CloneB", 
    "CloneC", "CloneD", "CloneE", "CloneF"), class = "factor"), 
    Score = c(4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 3, 3, 3, 3, 3, 4.5, 
    4.5, 4.5, 4.5, 4.5, 3.2, 3.2, 3.2, 3.2, 3.2, 5.5, 5.5, 5.5, 
    5.5, 5.5, 4.8, 4.8, 4.8, 4.8, 4.8, 5.8, 5.8, 5.8, 5.8, 5.8, 
    3.8, 3.8, 3.8, 3.8, 3.8, 5.3, 5.3, 5.3, 5.3, 5.3, 4, 4, 4, 
    4, 4, 6.3, 6.3, 6.3, 6.3, 6.3)), .Names = c("Weight", "Sex", 
"Age", "Genotype", "Score"), class = "data.frame", row.names = c(NA, 
-60L))


Vastaa

Sähköpostiosoitettasi ei julkaista. Pakolliset kentät on merkitty *