R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Lineaarisen mallin nopea sovittaminen

Tavanomainen lineaarinen regressio sovitetaan R:ssä komennolla lm(). Jos sovitettava malli sisältää enemmän kuin yhden selittävän muuttujan, on malli kätevintä sovittaa kaava-rajapintaa käyttäen, esimerkiksi lm(y ~ x1 + x2, data=dat), jossa y ~ x1 + x2 on siis se kaava, joka spesifioi sovitettavan mallin.

Kaavassa ~ -merkin vasemmalla puolella on ennustettava vastemuuttuja ja sen oikealla puolella luetellaan ennustamiseen käytettävät vastemuuttujaa selittävät muuttujat. Esimerkkimallissa siis muuttujan y arvoja ennustetaan muuttujien x1 ja x2 arvoilla. Kaavassa käytetyt muuttujan nimet ovat samat kuin analyysissä käytettävästä aineistosta (objekti dat) löytyvät sarakkeiden nimet.

Testataanpa muutamien erilaisten regressiomallien sovittamiseen soveltuvien menetelmien nopeutta.

Testiaineiston simulointi

Jos aineisto vain mahtuu R:n muistiin, ei lineaarisen regressiomallin sovittaminen pääsääntöisesti kestä kovin kauaa. Simuloidaanpa ensin kymmenen miljoonan rivi aineisto, jolla asiaa voidaan simuloida:

nr<-10000000  # rivejä
nc<-11        # sarakkeita
set.seed(1)   # siemenluku
 
m<-matrix(ncol=nc, nrow=nr, data=rnorm(nr*nc))
colnames(m)<-c("y", letters[1:10])
 
df<-data.frame(m)

Nyt objekti m sisältää aineiston matriisimuotoisena ja objekti df puolestaan data framena. Objektin tyypin voi aina tarkistaa komennolla class():

class(m)
[1] "matrix"
 
class(df)
[1] "data.frame"

Seuraavia vaiheita varten on tärkeää tietää, missä muodossa aineisto on.

Lineaarinen regressio

Sovitetaan lineaarinen regressiomalli data frame -muotoista objekti mm käyttäen, ja talletetaan tulos objektiin fit1:

fit1<-lm(y~., data=mm)

Kaavassa (Y ~ .) esiintyvä piste tarkoittaa sitä, että kun y-niminen sarake on vastemuuttuja, käytetään kaikkia muita samassa data framessa olevia sarakkeita selittävinä muuttujina.

Mallin tulokset saadaan tulostettua komennolla summary():

summary(fit1)
 
Call:
lm(formula = y ~ ., data = mm)
 
Residuals:
    Min      1Q  Median      3Q     Max 
-5.4282 -0.6743  0.0000  0.6745  5.4845 
 
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  4.038e-04  3.163e-04   1.277    0.202
a           -2.675e-04  3.163e-04  -0.846    0.398
b            1.015e-04  3.162e-04   0.321    0.748
c           -3.021e-04  3.162e-04  -0.955    0.339
d           -6.027e-05  3.164e-04  -0.191    0.849
e           -2.989e-04  3.163e-04  -0.945    0.345
f           -3.846e-04  3.163e-04  -1.216    0.224
g            2.250e-04  3.163e-04   0.711    0.477
h           -2.457e-04  3.163e-04  -0.777    0.437
i            4.233e-05  3.162e-04   0.134    0.894
j            2.830e-04  3.163e-04   0.895    0.371
 
Residual standard error: 1 on 9999989 degrees of freedom
Multiple R-squared:  6.066e-07, Adjusted R-squared:  -3.934e-07 
F-statistic: 0.6066 on 10 and 9999989 DF,  p-value: 0.8097

Mallin sovittamiseen kuluva aika voidaan selvittää komennolla system.time(), jonka sisään mitattava funktio sijoitetaan:

system.time(fit1<-lm(y~., data=mm))
 
   user  system elapsed 
  23.94    1.70   26.25

Tuloksena syntyy tuloste, joka kertoo mallin sovittamiseen kuluneen ajan. Tulosteen luvuista viimeisin eli ”elapsed” kertoo komennon ajamiseen kuluneen seinäkelloajan. Yllä olevan mallin sovittaminen siis vein 26,25 sekuntia.

R:ssä lineaarinen regressiomalli voidaan sovittaa myös komennolla lsfit(). Tällöin datan sisältävän objektin on oltava matriisi:

system.time(fit2<-lsfit(cbind(m[,2:11]), m[,"y"]))
 
   user  system elapsed 
   5.97    1.50    5.80 
 
ls.print(fit2)
 
Residual Standard Error=1.0002
R-Square=0
F-statistic (df=10, 9999989)=0.6066
p-value=0.8097
 
          Estimate Std.Err t-value Pr(>|t|)
Intercept    4e-04   3e-04  1.2766   0.2017
a           -3e-04   3e-04 -0.8458   0.3977
b            1e-04   3e-04  0.3211   0.7481
c           -3e-04   3e-04 -0.9554   0.3394
d           -1e-04   3e-04 -0.1905   0.8489
e           -3e-04   3e-04 -0.9451   0.3446
f           -4e-04   3e-04 -1.2158   0.2241
g            2e-04   3e-04  0.7113   0.4769
h           -2e-04   3e-04 -0.7768   0.4373
i            0e+00   3e-04  0.1339   0.8935
j            3e-04   3e-04  0.8947   0.3709

Mallin saa sovitettua vieläkin pikkuhiukan nopeammin, jos ei tarvitse muuta kuin selittävien muuttujien vaikutusten estimaatit. Komentoon lsfit() verrattuna lm.fit() vaatii selittävät tekijät koeasetelmamatriisina, jossa on tarpeen vaatiessa myös intecept-termi esitettynä yhdellä sarakkeella:

system.time(fit3<-lm.fit(cbind(1, m[,2:11]), m[,"y"]))
 
   user  system elapsed 
   4.39    0.91    5.30 
 
data.frame(fit2$coefficients)
 
  fit2.coefficients
       4.037930e-04
a     -2.675151e-04
b      1.015406e-04
c     -3.020950e-04
d     -6.027140e-05
e     -2.988859e-04
f     -3.845862e-04
g      2.249807e-04
h     -2.456845e-04
i      4.232559e-05
j      2.830163e-04

Logistinen malli

R:ssä logistinen malli sovitetaan komennolla glm(). Lineaarisen mallin sovittamisen verrattuna yleistetyn lineaarisen mallin tapauksessa on myös määritettävä linkkifunktio argumentilla family, jos se ei ole gaussinen.

Muutetaan ensin aiemmin tehtyä aineistoamme siten, että vastemuuttuja on binäärinen;

set.seed(1)
mm$y<-rbinom(10000000, 1, 0.5)

Tämän jälkeen sovitetaan logistinen malli, ja katsotaan sen antamat tulokset:

system.time(fit4<-glm(y~., family="binomial", data=mm))
 
summary(fit4)
Call:
glm(formula = y ~ ., family = "binomial", data = mm)
 
Deviance Residuals: 
   Min      1Q  Median      3Q     Max  
-1.182  -1.177   1.174   1.177   1.182  
 
Coefficients:
              Estimate Std. Error z value Pr(>|z|)  
(Intercept)  3.943e-05  6.325e-04   0.062   0.9503  
a            1.025e-05  6.324e-04   0.016   0.9871  
b            1.044e-03  6.323e-04   1.651   0.0987 .
c           -2.271e-04  6.323e-04  -0.359   0.7195  
d            5.170e-05  6.326e-04   0.082   0.9349  
e            1.070e-03  6.324e-04   1.692   0.0906 .
f            8.973e-04  6.325e-04   1.419   0.1560  
g           -4.689e-04  6.324e-04  -0.741   0.4584  
h            3.352e-04  6.324e-04   0.530   0.5961  
i           -2.573e-04  6.323e-04  -0.407   0.6841  
j           -1.090e-04  6.325e-04  -0.172   0.8632  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 13862944  on 9999999  degrees of freedom
Residual deviance: 13862935  on 9999989  degrees of freedom
AIC: 13862957
 
Number of Fisher Scoring iterations: 3

Logistisen mallin sovittaminen on oleellisesti hitaampaa kuin lineaarisen mallin sovittaminen. Paketti speedglm kuitenkin sisältää toiminnallisuutta, joka nopeuttaa tällaisen mallin sovittamista:

system.time(fit4<-speedglm(y~., family=binomial(), data=mm))
 
#   user  system elapsed 
#  28.53    4.69   33.29

Logistisen mallin sovittamista voidaan vieläkin nopeuttaa, jos siirrytään käyttämään Microsoftin RClientin (9.0.1) tarjoamaa rxGlm() -komentoa. rxGlm() ei osaa tulkita kaavaa, jossa esiintyy piste (.) selittävien muuttujien sijaan, joten koko kaava on kirjoitettava kokonaan auki. Komennossa tulee myös käyttää muutamia lisäargumentteja, jotta tulokset vastaavat R:n perus-glm:n parametrisointia, koska oletusarvoinen parametrisointi vastaa SAS:ia:

system.time(fit5<-rxGlm(y~a+b+c+d+e+f+g+h+i+j, family="binomial", data=mm, dropFirst=TRUE, dropMain=FALSE))
 
   user  system elapsed 
   0.23    0.53   12.42

Yhteenveto

Lineaarisen regressiomallin voi sovittaa kaikkein nopeimmin komennolla lsfit(). Logistisen mallin sovittaminen taas käy nopeimmin Microsoftin RClientin komennolla rxGlm(). Jos sitä ei ole saatavilla, sovitusta voi nopeuttaa myös speedglm-paketin komennolla speedglm().

Mainittakoon, että esitetyt nopeudet ovat suuntaa-antavia, koska tässä on raportoitu vain yhden toiston tuloksia. Toisettaessa eri menetelmien keskinäiset nopeussuhteet kuitenkin säilyvät, vaikkei vakuuttavampia, useampiin toistoihin perustuvia tuloksia olekaan esitetty.