R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Logistinen regressio

Sisällys

Mikä on logistinen regressio?

Hieman teoriaa
  Vedonlyöntisuhde
  Vetosuhde
  Vetosuhteen tilastollinen testaaminen
  Logistinen malli

Logistinen regressio R:ssä
  Aineisto
  Mallin sovittaminen
  Mallin sopivuuden mittarit
  Mallin tulosten tulkinta
  Mallin tulosten esittäminen
  Mallin ennusteiden muuntaminen todennäköisyyksiksi
  Ennustemallin toimivuuden arviointi

Mallin valinta
  Mallien sovittaminen
  Mallien vertailu
  Automaattiset menetelmät

Erikoistapauksia
  Ehdollinen logistinen regressio
  Firthin harhattomampi menetelmä
  Harrellin lrm-funktio
  Tulosten visuaalinen esittäminen

Mikä on logistinen regressio?

Logistista regressiota käytetään tapauksissa, joissa vastemuuttuja on binäärinen eli se voi saada vain kahdenlaisia arvoja. Esimerkiksi ihmisten sairausstatus (sairas / terve) on binäärinen muuttuja. Useimmiten binääristen muuttujien arvot merkitään nollalla (terve) ja ykkösellä (sairas). Logistisella regressiolla ennustetaan vastemuuttujan arvoa yksi.

Hieman teoriaa

Vedonlyöntisuhde

Logistisella regressiolla ennustetaan tapahtuman (sairas) vetoa. Vetoa kutsutaan myös vedonlyöntisuhteeksi. Vedonlyöntisuhde on kahden todennäköisyyden suhde. Esimerkiksi jos sairastumisen todennäköisyys (tai sairaiden suhteellinen osuus aineistossa) on p = 0.6, on sen vedonlyöntisuhde 1 / (1 – p) = 0.6 / (1 – 0.6) = 0.6 / 0.4 = 1.5. Toisinsanoen vedonlyöntisuhde on todennäköisyyden, että tapahtuma sattuu, ja todennäköisyyden että tapahtuma ei satu, suhde. Jos vedonlyöntisuhde (v) tunnetaan, voidaan sen perusteella laskea myös todennäköisyys kaavalla p = v / (1 + v). Jos vedonlyöntisuhde on yllä laskettu 1.5, on sitä vastaava todennäköisyys siis 1.5 / (1 + 1.5) = 1.5 / 2.5 = 0.6.

Vetosuhde

Vetosuhteen avulla ilmaistaan kahden vedonlyöntisuhteen suhdetta. Vetosuhdetta kutsutaan myös odds ratioksi (OR). OR vertaa siis suhteellisia osuuksia tai todennäköisyyksiä niiden vedonlyöntisuhteiden suhdeluvulla. Esimerkiksi tapausverrokkitutkimuksissa verrataan kahta otosta: sellaista, joka on altistunut jollekin sairastumisen kannalta oletetusti merkittävälle asialle (tapaukset) ja sellaista, joka ei ole altistunut (verrokit). Jos verrokkien joukossa sairastumisen todennäköisyys on 0.05 ja tapausten joukossa vastaava todennäköisyys on 0.10, on näiden vetosuhde ( 0.10 / (1 – 0.10) ) / ( 0.05 / (1 – 0.05) ) = 2.11. Vetosuhde ei siis suoraan kerro todennäköisyyksien tai suhteellisten osuuksien suhdetta vaan nimenomaan vedonlyöntisuhteiden suhteen.

Vetosuhde voi olla helpoin ymmärtää erikoistapauksessa, jossa sekä vastemuuttuja, että sen yksi selittäjä ovat molemmat kaksiarvoisia muuttujia.

Esimerkiksi Titanic-onnettomuuden aineistosta voidaan tuottaa menehtyneiden ja sukupuolen välinen ristiintaulukointi, jossa esitetyt numerot ovat matkastajien lukumääriä:

        sukupuoli
menehtyi nainen mies
       0    339  161
       1    127  682

Naisten vedonlyöntisuhde menehtyä on 127 / 339 = 0.375 ja miesten vastaavasti 682 / 161 = 4.24. Vetosuhde miesten menehtymiselle on siten 4.24 / 0.375 = 11.31. Samaan tulokseen päästään ristitulolla: OR = (339*682) / (127*161) = 11.31.

Vetosuhteen tilastollinen testaaminen

Täydellisyyden vuoksi esitetään tässä samassa yhteydessä myös taulukosta lasketun OR:n keskivirheen estimointi. OR:n keskivirhe logit-asteikolla on sqrt( 1/339 + 1/161 + 1/127 + 1/682 ) = 0.136.

Keskivirhettä käyttäen voidaan OR:lle laskea myös 95 % luottamusväli. Sitä varten yllä laskettu OR pitää ensin muuntaa logit-asteikolle. Tällöin log(odds) = log(11.31) = 2.425. Luottamusväli lasketaan lisäämällä tai vähentämällä log(odds) – arvosta 1.96 keskivirhe. Tässä tapauksessa luottamusväliksi logit-asteikolla tulee siis 2.425 – 1.96 * 0.136 … 2.425 + 1.96 * 0.136 = 2.158 – 2.692. Luottamusväli voidaan lopuksi muuntaa samalle asteikolle kuin odds ratio. Tällöin luottamusväliksi tulee exp(2.158) … exp(2.692) = 8.65 – 14.76. Usein tulos esitetään muodossa or = 11.31 (95 % CI = 8.65 – 14.76). Koska OR:n luottamusväli eroaa ykkösestä (logit-asteikolla nollasta), on tulos myös tilastollisesti merkitsevä. Jos luottamusväliksi haluttaisiin vaikkapa 99 % luottamusväli, tulisi kertoimena käyttää 1.96:n sijaan 2.57:ää. Kertoimien arvot tulevat helposti normaalijakauman kvantiilifunktiosta arvolla (1 – ((1 – cipit) / 2)), jossa cipit on tavoitellun luottamusvälin pituus (tässä ci = 0.95) eli R:ssä qnorm(0.975).

Kun keskivirhe on tiedossa, voidaan OR:lle laskea myösp-arvo. Tyypillisesti p-arvon laskenta perustetaan z-testiin, jollloin testisuureen arvo saadaan jakamalla logit-asteikolla olevan OR:n arvo sen keskivirheellä. Tässä tapauksessa testisuureen arvoksi saadaan siis z = 2.425 / 0.136 = 17.83. Testisuureen arvon perusteella voi selvittää p-arvon esimerkiksi kriittisten arvojen taulukoita käyttäen, mutta R:ssä p-arvo voidaan laskea valmiita funktioita käyttäen: 2*(1-pnorm(17.83)). Tässä tapauksessa testisuureen arvoon niin suuri, että sitä vastaava p-arvo on käytännössä nolla. Tulkintahan olisi sitten jotakuinkin, että miehillä oli naisia paljon suurempi todennäköisyys menehtyä Titanicin upotessa.

Logistinen malli

Miten vedonlyöntisuhde sitten liittyy logistiseen regressioon? Logistisessa regressiossa mallinnetaan vedonlyöntisuhdetta, tarkkaan ottaen vedonlyöntisuhteen logaritmia. Vedonlyöntisuhteen logaritmiä käytetään siksi, että logistinen funktio voi ottaa syötteenä minkä tahansa reaaliluvun, ja sen tulos on aina väliltä 0-1. Tämä on näppärää siksi, että todennäköisyydet, joita vedonlyöntisuhteessakin käsitellään, saavat vain arvoja väliltä 0-1.

Logistinen regressiomalli kahdelle selittäjälle voidaan kirjoittaa muotoon:

log(p / 1 – p) = b0 + b1*x1 + b2*x2

Tai vastaavasti muotoon

(p / 1 – p) = exp(b0 + b1*x1 + b2*x2)

Tällöin mallin sanotaan antavan ennusteet log odd tai logit -skaalalla.

Entä miten vetosuhde sitten tulee kuvioihin? Jos vastemuuttujan vedonlyöntisuhde (odd) selittävien muuttujien suhteen on esimerkiksi exp(b0 + b1*x1), saadaan (ehdollinen) vedonlyöntisuhde yhdelle selittävän muuttujan asteikon muutokselle siuhteesta exp(b0 + b1*(x+1)) / exp(b0 + b1*x1). Regressiomallin parametrin estimaatti saa näin tulkinnakseen, että kutakin selittävän muuttujan x1 yhden pykälän muutosta kohden odds kasvaa exp(b1):llä. Logistinen regressio voidaan siis käsittää OR:n yleistyksenä useampien selittävien muuttujien tapaukseen.

Logistinen malli kuuluu yleistettyjen lineaaristen mallien perheeseen. Logistisen regression tapauksessa linkkifunktio, joka kuvaa selittävien muuttujien lineaarikombinaation ja vastemuuttujan odotusarvon välistä suhdetta, on logit- tai logistinen funktio.

Englannikielinen Wikipedia antaa runsaasti lisätietoa sekä vetosuhteesta että logistisesta regressiosta.

Seuraavaksi tutustutaan logistisen regression sovittamiseen liittyviin käytänteisiin R:ssä.

Logistinen regressio R:ssä

Aineisto

Seuraavassa käytännön työvaiheita kuvaavassa osuudessa käytetään esimerkkiaineistona Titanicin matkustajien selvityvyyttää kuvaavaa aineistoa. Aineisto on saatavilla Valderbiltin yliopiston sivuilta. Aineiston muuttujat on kuvattu erillisellä sivulla.

Aineiston ladataan R:ään esimerkiksi sueraavasti:

library(rio)
d <- import("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/titanic3.csv")

Aineistosta löytyy muuttuja survived, joka kertoo, selvisikö matkustaja matkasta vai ei. Haluamme kuitenkin ennustaa matkaustajien riskiä menehtyä, joten luodaan aineistoon uusi muuttuja seuraavasti:

d$died <- ifelse(d$survived==1, 0, 1)

Näin luotua muuttujaa käytetään vastemuuttujana seuraavissa esimerkeissä.

Mallin sovittaminen

Logistinen regressiomalli sovitetaan R:ssä funktiolla glm(). Esimerkiksi edellä esitettyä ristiintaulukointia vastaava regressiomalli voidaan sovittaa seuraavasti:

fit1 <- glm(died ~ sex, data=d, family=binomial)

Funktiolle glm() määritellään sovitettavan mallin kaava (died ~ sex), jossa ~-merkin vasemmalla puolella on mallin vastemuuttuja, ja sen oikealla puolella mallin selittävät muuttujat. Lisäksi on tarpeen kertoa, mistä data.frame -tyyppisestä objektista R:n muistissa kaavassa mainitut muuttujat löytyvät (data=d). Viimeisenä on kerrottu, millainen malli sovitetaan (family=binomial). Jos valittu malliperhe on binomial, on oletusarvoinen linkkifunktio logit, joten sitä ei tarvitse erikseen määritellä.

Mallin sovituksen jälkeen tulokset voidaan tulostaa ruudulle seuraavasti:

summary(fit1)
 
 
Call:
glm(formula = died ~ sex, family = binomial, data = d)
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.8196  -0.7977   0.6511   0.6511   1.6124  
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9818     0.1040  -9.437   <2e-16 ***
sexmale       2.4254     0.1360  17.832   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 1741.0  on 1308  degrees of freedom
Residual deviance: 1368.1  on 1307  degrees of freedom
  (1 observation deleted due to missingness)
AIC: 1372.1
 
Number of Fisher Scoring iterations: 4

Tuloste koostuu neljästä osasta. Ensimmäinen osa (Call) kertoo sovitetun mallin. Toinen osa (Residuals) antaa mallin sovitteen residuaalien tunnusluvut. Kolmas osa (Coefficients) on mallin sovituksen varsinainen tulos. Neljäs osa puolestaan luettelee mallin sopivuuteen liittyviä tunnuslukuja, kuten deviancen ja siitä johdetun AIC-arvon. AIC-arvon perusteella voidaan tehdä muun muassa mallin valintaa.

Mallin sopivuuden mittarit

Yleensä on tavoitteena sovittaa aineistoon malli, joka kuvaa aineiston piirteitä huvin. Mallin sopivuutta voidaan mitata sen logLikelihood-arvolla ja siitä johdetuilla tunnusluvuilla, kuten deviancella, AIC-arvolla ja selitysasteella. Deviancen perusteella voidaan johtaa myös uskottavuusosamäärätesti mallin yleiselle sopivuudelle aineistoon. Malli sopii aineistoon sitä paremmin mitä lähempänä nollaa mallin residual deviance tai sen AIC-arvo on. Selitysasteen tapauksessa suurempi arvo kielii mallin paremmasta sopivuudesta.

Mallille voidaan laskea log likelihood, logL (uskottavuus). R:ssä mallin uskottavuus voidaan selvittää seuraavasti: logLik(fit1). Null deviance puolestaan johdetaan mallin uskottavuudesta seuraavasti: Null deviance = 2 * (logL(saturated model) - logL(null model)). Residual deviance on vastaavasti 2 * (logL(saturated model) - logL(fitted model)). Saturoidussa mallissa on yksi parametri kutakin havaintoa kohden. Null-mallissa puolestaan on vain yksi parametri (selittäjänä vain keskiarvo eli Intercept). Residual deviance kuvaa mallin selittämättä jättänyttä hajonnan määrää.

Jos siis ensin sovitetaan R:ssä seuraavat kaksi mallia, ja tulostetaan niiden sekä mielenkiintoisen mallin log likelihood:

# Saturoitunut malli
fit1.satur <- glm(died ~ factor(1:nrow(d)), data=d, family=binomial)
 
# Null-malli
fit1.null <- glm(died ~ 1, data=d, family=binomial)
 
logLik(fit1.satur)
'log Lik.' -3.796913e-09 (df=1309) # Käytännössä nolla
 
logLik(fit1.null)
'log Lik.' -870.5122 (df=1)
 
logLik(fit1)
'log Lik.' -684.0515 (df=2)

Voidaan seuraavaksi laskea tulosteestakin saatava Null deviance, joka on siis 2 * (-3.796913e-09 - (-870.5122)) = 1741.024. Residual deviance on vastaavasti 2 * (-3.796913e-09 - (-684.0515)) = 1368.103.

AIC-arvo puolestaan saadaan laskettua mallin parametrien määrästä ja mallin uskottavuudesta kaavalla 2 * k - 2 * logL. Mallin fit1 tuloksista laskettuna AIC-arvo on siis 2 * 2 - 2 * (-684.0515) = 1372.103.

Mallin yleiselle sopivuudelle voidaan laskea myös p-arvo. R:n funktiolla anova voidaan tehdä mallin sopivuudelle likelihood ratio test eli uskottavuusosamäärätesti:

anova(fit1, test="LRT")
 
 
Analysis of Deviance Table
 
Model: binomial, link: logit
 
Response: died
 
Terms added sequentially (first to last)
 
 
     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                  1308     1741.0              
sex   1   372.92      1307     1368.1 < 2.2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Koska testin antama p-arvo on hyvin pieni, on mallin yleinen sopivuus tilastollisesti merkitsevä. Sama testi voidaan laskea myös käsin. Testisuure on mallien devianssien erotus, ja testisuureen vapausasteiden määrä on null-mallin ja sovitetun mallin vapausateiden erotus. Testisuureen arvoksi saadaan siten 1741.024 - 1368.103 = 372.9213. Vapausasteiden määrä on tässä tapauksessa 1. R:ssä sama voidaan tehdä seuraavasti:

# testisuure
lrt <- deviance(fit1.null) - deviance(fit1) # 372.9213
 
# vapausasteet
lrt.df <- df.residual(fit1.null) - df.residual(fit1) # 1
 
# p-arvo
pchisq(lrt, lrt.df, lower.tail=FALSE)
 
[1] 4.326495e-83

Lineaarisen regressiomallin ollessa kyseessä mallin sopivuuden mittarina käytetään usein selitysastetta (R Squared, R2). Yleistettyjen lineaaristen mallien tapauksessa käytetään pseudoselitysastetta. Selitysaste voidaan laskea useilla eri menetelmillä. Tyypillisesti raportoidaan Nagelkerken (all olevassa tulosteessa r2CU) tai McFaddenin selitysaste. Nämä voidaan R:ssä laskea seuraavasti:

library(pscl)
pR2(fit1)
 
         llh      llhNull           G2     McFadden         r2ML         r2CU 
-684.0515253 -870.5121915  372.9213323    0.2141965    0.2479032    0.3370385

Mallin tulosten tulkinta

Yleensä mallin kiinnostavin tulos on regressiotaulukko:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9818     0.1040  -9.437   <2e-16 ***
sexmale       2.4254     0.1360  17.832   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Taulukossa kullekin mallin parametrille on omistettu yksi rivi. Tulosteesta nähdään, että mallissa on kaksi parametriä: Intercept ja sukupuoli. Sukupuoli on listattu muodossa 'sexmale', josta voidaan päätellä, että sukupuoli 'female' on tämä luokittelumuuttujan verrokkitaso, ja taulukossa ilmoitettu tulos on siis miessukupuolelle verrattuna naissukupuoleen.

Kullekin parametrille on neljä saraketta. Estimate kertoo parametrin estimaatin arvon logit-asteikolla. Std. Error on parametrin estimaatin keskivirhe. Sarakkeessa z value on z-testisuureen arvo, joka on saatu jakamalla estimaatin arvo sen keskivirheellä. Sarakkeesa Pr(>|z|) on z-testisuureen arvoa vastaava p-arvo.

Mallin tulosten esittäminen

R:ssä on valmiita funktioita, joilla parametrien estimaateille voidaan estimoida luottamusvälit. Käyttäjän tehtäväksi jää tulosten muuntaminen logit-asteikolta odds ratioiksi. Yllä esitetyssä taulukossa olevat tulokset perustuvat normaalisuusoletukseen. R:n funktiolla confint() voidaan laskea parametrien estimaateille luottamusvälit. Ongelmia tuottaa se pikkuseikka, että funktion oletusarvoisesti laskemta luottamisvälit perustuvat profile likelihoodiin, joka ei menetelmällisesti vastaa summary()-funktion antamia tuloksia.

Normaalijakaumaoletukseen perustuvat luottamusvälit saadaan laskettua seuraavasti:

confint.default(fit1)
 
                2.5 %     97.5 %
(Intercept) -1.185723 -0.7779028
sexmale      2.158845  2.6920316

Tulokset saadaan OR-asteikolle ja kauniimpaan taulukkoon vaikkapa seuraavasti:

fit1coef <- cbind(exp(coef(fit1)), exp(confint.default(fit1)), round(coef(summary(fit1))[,4], 3))
colnames(fit1coef) <- c("OR", "CI 2.5%", "CI 97.5%", "p-value")
round(fit1coef, 3)
 
                OR CI 2.5% CI 97.5% p-value
(Intercept)  0.375   0.306    0.459       0
sexmale     11.307   8.661   14.762       0

Jos halutaan käyttää profile likelihood-arvoja, pitää taulukko koostaa seuraavasti:

fit1coef<-cbind(exp(coef(fit1)), exp(confint(fit1)), round(anova(fit1, test="LRT")$"Pr(>Chi)", 3))
colnames(fit1coef)<-c("OR", "CI 2.5%", "CI 97.5%", "p-value")
round(fit1coef, 3)
 
                OR CI 2.5% CI 97.5% p-value
(Intercept)  0.375   0.305    0.458      NA
sexmale     11.307   8.687   14.809       0

Yllä olevista tuloksista on helppo havaita, että normalisuusoletukseen ja profile likelihoodiin perustuvat tulokset eroavat toisistaan hieman. Profile likelihoodin antamia tuloksia pidetään tarkempina kuin normaalisuusoletukseen perustuvia tuloksia. Normaalisuusoletukseen perustuvat tulokset voidaan laskea mallin tulosteesta käsin, profile likelihoodiin perustuvia tuloksia sen sijaan ei voida.

Mallin ennusteiden muuntaminen todennäköisyyksiksi

Mallin avulla voidaan ennustaa kullekin aineistossa olevalle havainnolle arvo logit-asteikolla tai todennäköisyytenä, että havainto kuuluu luokkaan yksi eli että matkustaja on menehtynyt. Ennusteet saadaan R:ssä funktiolla predict():

# Ennusteet logit-asteikolla
pred1.logit <- predict(fit1)
 
# Ennusteet
pred1.prob <- predict(fit1, type="response")

Tässä tapauksessa sekä logit-ennusteet että todennäköisyysennusteet saavat vain kahdenlaisia arvoja:

table(pred1.logit)
 
pred1.logit
-0.981813020919238   1.44362529285895 
               466                843 
 
table(pred1.prob)
 
pred1.prob
0.272532188841721 0.809015421115058 
              466               843

Sama voidaan laskea myös käsin. Kun mallin tulos näyttää seuraavalta:

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9818     0.1040  -9.437   <2e-16 ***
sexmale       2.4254     0.1360  17.832   <2e-16 ***

voidaan siitä laskea log odds menehtyä, jos sukupuoli on mies, seuraavasti: -0.9818 + 2.4254 = 1.4436. Log odd (logit) voidaan kääntään odds-asteikolle seuraavasti: e^1.4436 eli R:ssä exp(1.4436) = 4.235918. Odds-asteikolta tulos taas voidaan muuntaa todennäköisyydeksi seuraavasti: 4.235918 / (1 + 4.235918) = 0.8090115. Käsin lasketut arvot vastaavat jälleen R:n funktioilla automaattisesti saatavia arvoja.

Ennustemallin toimivuuden arviointi

Edellä on sovitettu ennustemalli, jonka tavoitteena on ennustaa matkustajien kuolleisuutta mahdollisimman hyvin. Ennustemallin ollessa kyseessä pyritään usein minimoimaan väärien ennusteiden määrä. Tätä voidaan myös käyttää mittarina mallin toimivuudelle. Hyvin perinteinen tapa arvioida mallin toimivuutta on Receiver Operating Characteristic (ROC) -käyrä.

Seuraava R-koodi piirtää ROC-käyrän sovitetusta mallista paketin pROC funktioilla. Lisäksi koodi arvioi ROC-käyrän alle jäävän alan. Jotta malli ennustaisi paremmin kuin satunnaisuus, tulee ROC-käyrän erota kuvassa näkyvältä lävistäjältä, ja sen alle jäävän pinta-alan tulee olla suurempi kuin 0.5.

library(pROC)
roc1 <- roc(fit1$model$died, predict(fit1))
plot(roc1)
auc(roc1)
 
Area under the curve: 0.7605

Koodissa on käytetty sovitetun mallin sisältävän R-objektin aineistoa, koska se on varmasti samanmittainen kuin mallin antamien ennusteiden määrä. Pituus voi erota alkuperäisen aineiston pituudesta, koska mallin sovituksen aikana mahdolliset puuttuvat havainnot poistetaan aineistosta automaattisesti. Näin muodostunut, suodatettu aineisto talleteaan malliobjektin osaksi.

Piirretty ROC-käyrä näyttää seuraavalta:

roc1

Käyrän alle jäävä ala on tässä 0.7605. Molemmat tukevat siis oletusta, että malli ennustaa paremmin kuin kolikonheitto.

Mallin ennustavuutta voidaan tutkia myös ristiintaulukointien avulla. Taulukkoa, jossa on esitetty havaittu vastemuuttuja ennustettua vasteen arvoa vasten, kutsutaan usein nimellä confusion matrix. Logistisesta mallista saatavat ennusteet ovat todennäköisyyksiä väliltä 0-1, joten mainittua konfuusiomatriisia varten todennäköisyydet on luokiteltuva samoihin luokkiin kuin vastemuuttujakin, esimerkiksi arvoiksi kyllä/ei tai 1/0. Luokittelua varten on käytettävä jotakin raja-arvoa, jota suuremmat arvot sijoitetaan luokkaan 1 ja sitä pienemmät luokkaan 0. Tässä dikotomisoinnissa käytettävä on raja-arvo on usein 0.5 tai ykkösluokan osuus alkuperäisessä aineistossa. Dikotomisointi ja konfuusiomatriisi voidaan tuottaa seuraavasti:

# Dikotomisointi
pred05 <- ifelse(predict(fit1) > 0.5, 1, 0)
cutoff <- sum(fit1$model$died)/nrow(fit1$model)
predco <- ifelse(predict(fit1) > cutoff, 1, 0)
 
# Konfuusiomatriisit
ctab1 <- table(obs=fit1$model$died, pred=pred05)
ctab1
   pred
obs   0   1
  0 339 161
  1 127 682
 
ctab2 <- table(obs=fit1$model$died, pred=predco)
ctab2
   pred
obs   0   1
  0 339 161
  1 127 682

Konfuusiomatriisista voidaan laskea monia erilaisia tunnuslukuja. Perinteisiä konfuusiomatriisista laskettuja tunnuslukuja ovat sensitiivisyys, spesifisyys ja tarkkuus (accuracy). Tarkkuus on oikeiden postiivisten ja negatiivisten ennusteiden suhteellinen osuus kaikista ennusteita. Spesifisyys on oikeiden negatiivisten ennusteiden suhteellinen osuus kaikkien negatiivisten ennusteiden joukossa. Sensitiivisyys on vastaavalla tavalla oikeiden positiivisten ennusteiden suhteellinen osuus kaikkien positiivisten ennusteiden joukossa. Esimerkiksi R:n laajennuspaketti caret tarjoaa valmiit funktiot näiden laskemiseen, mutta yllä olevasta konfuusiomatriisista nämä voidaan laskea myös seuraavasti:

accuracy <- (ctab1[1,1] + ctab1[2,2]) / sum(ctab1)
sensitivity <- ctab1[2,2] / sum(ctab1[2,])
specificity <- ctab1[1,1] / sum(ctab1[1,])
 
accuracy
[1] 0.7799847
 
sensitivity
[1] 0.8430161
 
specificity
[1] 0.678

Tarkkuuden, sensitiivisyyden ja spesifisyyden lisäksi on olemassa monia muitakin mahdollisia tapoja arvioida eri raja-arvojen toimivuutta. R:ssä muun muassa paketit OptimalCutpoints ja PresenceAbsence tarjoavat tapoja määrittää todennäköisyyksien dikotomisointia varten sopivia raja-arvoja varsin monenlaisilla eri menetelmillä. Paketteja voidaan käyttää seuraavasti:

library(OptimalCutpoints)
methods<-c("MaxSp", "MaxSe")
tmp<-data.frame(observed=fit1$model$died, predicted=predict(fit1, type="response"))
cp <- optimal.cutpoints(X = "predicted", status = "observed", tag.healthy = 0, methods = methods, data = tmp,  pop.prev = NULL, control = control.cutpoints(), ci.fit = TRUE)
print(cp)
 
Call:
optimal.cutpoints.default(X = "predicted", status = "observed", 
    tag.healthy = 0, methods = methods, data = tmp, pop.prev = NULL, 
    control = control.cutpoints(), ci.fit = TRUE)
 
Optimal cutoffs:
   MaxSp  MaxSe
1 0.8090 0.2725
 
Area under the ROC curve (AUC):  0.761 (0.736, 0.785) 
 
 
library(PresenceAbsence)
tmp<-data.frame(plotID=1:nrow(fit1$model), Observed=fit1$model$died, predicted=predict(fit1, type="response"))
optimal.thresholds(tmp)
 
         Method predicted
1       Default  0.500000
2     Sens=Spec  0.540000
3  MaxSens+Spec  0.540000
4      MaxKappa  0.540000
5        MaxPCC  0.540000
6  PredPrev=Obs  0.540000
7       ObsPrev  0.618029
8      MeanProb  0.618029
9    MinROCdist  0.540000
10      ReqSens  0.270000
11      ReqSpec  0.810000
12         Cost  0.540000

Molemmat paketit vaativat erillisen taulukon koostamista raja-arvojen laskemiseksi. Paketti OptimalCutpoints tarjoaa laajemman valikoiman erilaisia menetelmiä, mutta kaikki eivät toimi jokaisessa tapauksessa, joten menetelmät pitää tuntea, jotta tietää milloin niitä voi tai kannattaa käyttää. Paketin PresenceAbsence antamassa tulosteessa on 12 erilaista menetelmää. Esimerkiksi MaxSens+Spec vastaa menetelmää, jossa raja-arvo määritetään siten, että se maksimoi sensitiivisyyden ja spesifisyyden summan. Toinen huomionarvoinen menetelmä on MicROCdist, joka määrittää raja-arvoksi sen arvon, jolla yllä olevan ROC-käyrän ja kuvan vasemman ylänurkan etäisyys tulee kaikkein pienimmäksi. Jos siis haluttaisiin maksimoida sovitetun malli ennustusarvo, voitaisiin raja-arvona käyttää todennäköisyysarvoa 0.54.

pred054 <- ifelse(predict(fit1) > 0.54, 1, 0)
ctab1 <- table(obs=fit1$model$died, pred=pred054)
ctab1
 
   pred
obs   0   1
  0 339 161
  1 127 682

Tällä raja-arvolla laskettu tulos ei eroa muista aiemmin esitetyistä konfuusiomatriiseista, mutta muissa tapauksissa raja-arvon valinnalla voi olla suuri merkitys. Aina todennäköisyyksiä ei tietysti edes kannata dikotomisoida, vaan käyttötarkoituksiin soveltuvat raa'at todennäköisyydetkin.

Kolmas tapa tutkia mallin toimivuutta on ristiinvalidointi (cross-validation). Ristiinvalidoinnissa alkuperäinen aineisto jaetaan osiin, joista tiettyyn osaa sovitetaan malli, ja toisella osalla testataan sen sopivuutta. Kaikkein laajin ristiinvalidointi on leave-one-out (LOO), jossa malli sovitetaan miltei koko aineistoon, ja sitä testaan yhdellä sovitusaineistosta pois jätetyllä havainnolla.

LOO-ristiinvalidointi voidaan tehdä seuraavasti:

# LOO
library(boot)
cv1 <- cv.glm(d, fit1)
 
# Tulos
cv1$delta
[1] 0.1706243 0.1706241

Deltan arvo kertoo ristivalidoinnin virheen. Mitä pienempi virhe on, sitä parempi on mallin sopivuus aineistoon.

Mallin valinta

Edellä on kuvattu yhden logistisen regressiomallin sovittamiseen liittyviä seikkoja. Jos aineistossa on useampia selittäviä muuttujia, joita voidaan käyttää vastemuuttujan selittämiseen erilaisina yhdistelminä, tulee aiheelliseksi myös kysymys siitä, mikä selittävien muuttujien yhdistelmä antaa parhaan tuloksen. Erilaisten vaihtoehtoisten selitysmallien vertailua ja parhaan mallin valitsemista kutsutaan nimellä mallin valinta.

Jos tavoitteena on hyvä ennustemalli, ei mallin parsimonisuus ole kovin keskeinen kysymys. Parsimonialla tarkoitetaan tässä tilannetta, jossa mallista pyritään tekemään mahdollisimman yksinkertainen siten, ettei sen toimivuus kuitenkaan vaarannu. Toisaalta malli ei saa sopia aineistoon liian hyvin, ettei sen yleistettävyys kärsi. Näiden vaatimusten ristipaineissa pyritään siis tekemään paras mahdollinen malli.

Periaatteena voidaan pitää ajatusta, että kutakin mallilla tutkittavaa hypoteesia tai kuvattaa ilmiötä varten tehdään oma malli. Yleensä mallin perustaksi on jo aiempaa tietoa ilmiöstä, joten sitä kannattaa soveltaa mallia valittaessa. Esimerkiksi jos tutkitaan yleistä kuolleisuutta väestössä, kannattaa ihan perustavaa laatua olevaan malliin ottaa mukaan ainakin ikä ja sukupuoli, koska ne vaikuttavat voimakkaasti kuolleisuuteen. Automaattisilla mallin luontiin soveltuvilla menetelmillä voidaan tietysti toimia, jos ilmiö on täysin ennalta tuntematon, mutta silloinkin kannattaa pitää järki kädessä: automaattiset menetelmät eivät välttämättä tuota kunnollista mallia.

Kun on sovitettu useampia malleja, niitä voidaan vertailla keskenään erilaisilla menetelmillä. Mallin yleisen aineistoon istuvuuden arviointiin voidaan käyttää vaikkapa Akaike Information Criterion (AIC) -arvoa, ja jos tavoitteena on ennustaa hyvin sekä negatiivisiä että positiivisia tapauksia, voidaan valinnassa käyttää esimerkiksi tarkkuutta (accuracy). Jos taas halutaan korostaa positiivisten tai negatiivisten arvojen oikeaa ennustamista, voidaan käyttää myös sensitiivisyyttä ja spesifisyyttä.

Katsotaanpa seuraavaksi miten mallien vertailu R:ssä onnistuu.

Mallien sovittaminen

Jos on malleja sovitettaessa on käytettävissä valmiita hypoteeseja siitä, mikä selittää aineiston kuvaamaa ilmiötä, kannattaa tätä tietämystä käyttää hyväksi malleja sovietttaessa. Esimerkiksi Titanicin osalta muistelen, että ainakin matkustajan sukupuoli ja matkustusluokka selittivät kuolleisuutta. Mutta vaikuttavatko sukupuoli ja matkustusluokka yhdessä kuolleisuuteen siten, että esimerkiksi kolmosluokan miesmatkustajilla on suurempi kuolleisuus kuin ykkösluokan mies- tai naismatkustajilla?

Mallien vertailu ei onnistu kaikilla mittareilla, ellei kaikkia malleja varten ole käytössä tismalleen sama aineisto. Muodostetaan siis vain tarvittavista muuttujista uusi aineisto, josta puuttuvat arvot on poistettu (voitaisiin käyttää myös imputointia):

d2<-na.omit(d[,colnames(d) %in% c("died", "sex", "pclass", "age", "sibsp", "parch", "embarked")])

Nämä mallit voidaan sovittaa ainakin kolmeksi eri malliksi seuraavasti:

fit1 <- glm(died ~ sex, data=d2, family=binomial)
fit2 <- glm(died ~ pclass, data=d2, family=binomial)
fit3 <- glm(died ~ sex + pclass, data=d2, family=binomial)
fit4 <- glm(died ~ sex * pclass, data=d2, family=binomial)

Mutta aineistossa on muitakin muuttujia, kuten ikä, sisarusten ja vanhempien lukumäärä, matkan lähtöpaikka, matkan hinta ja hytin numero. Miten nämä vaikuttavat tuloksiin, jos ne lisätään malleihin? Osan näistä muuttujista sisältäviä malleja voidaan tehdä, esimerkiksi vaikkapa seuraavat:

fit5 <- glm(died ~ sex + pclass + age, data=d2, family=binomial)
fit6 <- glm(died ~ sex + pclass + age + sibsp, data=d2, family=binomial)
fit7 <- glm(died ~ sex + pclass + age + sibsp + parch, data=d2, family=binomial)
fit8 <- glm(died ~ sex + pclass + age + sibsp + parch + embarked, data=d2, family=binomial)
fit9 <- glm(died ~ (sex + pclass + age + sibsp + parch + embarked)^2, data=d2, family=binomial)

Viimeisin malli sisältää kaikki mallissa mukana olevien muuttujien väliset yhteisvaikutukset, ja se antaa varoituksen, että jotkut havainnot saivat tismalleen nollan tai ykkösen todennäköisyysennusteen. Tämä ei ole oikein hyvä juttu, koska se kielii siitä, että on liian monimutkainen.

Mallien vertailu

Edellä sovitettuja malleja voidaan vertailla esimerkiksi AIC-arvojen perusteella:

AIC(fit1, fit2, fit3, fit4, fit5, fit6, fit7, fit8, fit9)
 
     df       AIC
fit1  2 1104.8679
fit2  2 1308.0756
fit3  3 1019.4026
fit4  4  984.4183
fit5  4  990.3451
fit6  5  981.1984
fit7  6  982.5930
fit8  8  970.8747
fit9 28  940.3910

Malli 9 näyttää AIC-arvojen perusteella parhaalta. Entä miten mallit suhtautuvat toisiinsa, jos mittarina käytetään vaikkapa tarkkuutta? Tämähän voidaan tehdä vaikkapa konfuusiomatriiseista, joten keitetäänpä kokoon pieni funktio, jolla asia voidaan selvittää, ja sovelletaan funktiota sitten jokaiseen malliin:

# Havaittu kuolleiden suhteellinen osuus
obsfreq <- sum(d2$died)/nrow(d2)
 
acc <- function(fit, freq, col) {
   dat <- fit$model
   obs <- dat[,colnames(dat) %in% col]
   pred <- ifelse(predict(fit) > freq, 1, 0)
   tab <- table(obs=obs, pred=pred)
   acc <- sum(diag(tab)) / sum(tab)
   print(tab)
   acc
}
 
acc(fit1, obs, "died")
   pred
obs   0   1
  0 290 135
  1  96 523
[1] 0.7787356
 
acc(fit2, obs, "died")
   pred
obs   0   1
  0 179 246
  1 249 370
[1] 0.5258621
 
acc(fit3, obs, "died")
   pred
obs   0   1
  0 290 135
  1 194 425
[1] 0.6848659
 
acc(fit4, obs, "died")
   pred
obs   0   1
  0 218 207
  1 194 425
[1] 0.6159004
 
acc(fit5, obs, "died")
   pred
obs   0   1
  0 298 127
  1 212 407
[1] 0.6752874
 
acc(fit6, obs, "died")
   pred
obs   0   1
  0 298 127
  1 219 400
[1] 0.6685824
 
acc(fit7, obs, "died")
   pred
obs   0   1
  0 299 126
  1 218 401
[1] 0.6704981
 
acc(fit8, obs, "died")
   pred
obs   0   1
  0 300 125
  1 204 415
[1] 0.6848659
 
acc(fit9, obs, "died")
   pred
obs   0   1
  0 294 131
  1 144 475
[1] 0.73659

Mallin tarkkuuden mittareiden perusteella näyttää siltä, että fit1, joka sisältää vain sukupuolen vaikutuksen, on näistä kaikkein paras malli.

Riippuen siitä, onko tarkoituksena tuottaa ennustava malli vai selittävä malli, kohdistuisi valinta joko malliin fit1 tai huomattavasti monimutkaisempaan malliin fit9. Mallin fit1 ennustevirhe on alhainen, mutta sen kokonaissopivuus aineistoon näyttää heikommalta kuin mallin fit9.

Automaattiset mallinvalintamenetelmät

Edellä sovitettiin varsin monimutkainen malli, jossa olivat mukana kaikkien muuttujien väliset kahdenväliset interaktiot muuttujien päävaikutusten lisäksi. Tällaisesta mallista voidaan aloittaa ns. backward elimination, jossa automaattisesti tutkitaan, mikä malli antaa parhaan AIC-arvon, kun parametrejä (muuttujia ja niiden yhteisvaikutuksia) poistetaan yksi kerrallaan mallista. Usein tällaista menetelmää suositellaan käytettäväksi korkeintaan silloin, kun halutaan mahdollisesti yksinkertaistaa jo teoriapohjalta muodostettua suhteellisen monimutkaista mallia.

Käytännössä taaksepäin askeltava mallinvalinta on helppo toteuttaa:

fit10 <- step(fit9)
summary(fit10)
 
Call:
glm(formula = died ~ sex + pclass + age + sibsp + parch + embarked + 
    sex:pclass + sex:age + pclass:age + pclass:sibsp + age:parch, 
    family = binomial, data = d2)
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7487  -0.5593   0.4240   0.6638   3.0869  
 
Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)    -4.745759   1.166560  -4.068 4.74e-05 ***
sexmale         4.210583   0.997176   4.223 2.42e-05 ***
pclass          1.186486   0.393073   3.018 0.002540 ** 
age            -0.028315   0.026338  -1.075 0.282355    
sibsp          -0.892949   0.427709  -2.088 0.036820 *  
parch          -0.839499   0.276899  -3.032 0.002431 ** 
embarkedQ       1.294013   0.446535   2.898 0.003757 ** 
embarkedS       0.791268   0.217267   3.642 0.000271 ***
sexmale:pclass -1.123812   0.297775  -3.774 0.000161 ***
sexmale:age     0.039952   0.016544   2.415 0.015742 *  
pclass:age      0.012902   0.009018   1.431 0.152515    
pclass:sibsp    0.520382   0.162231   3.208 0.001338 ** 
age:parch       0.023634   0.008040   2.940 0.003286 ** 
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 1411.03  on 1043  degrees of freedom
Residual deviance:  894.06  on 1031  degrees of freedom
AIC: 920.06
 
Number of Fisher Scoring iterations: 6

Taaksepäin askeltavalla mallinvalinnalla löydettiin käytännössä vielä aiempaa fit9-mallia pienemmän AIC-arvon antava malli fit10. Mallin ennustekyky ei kuitenkaan ole yhtä hyvä kuin pelkän sukupuolen sisältävän mallin fit1 tai mallinvalinnan alkupiste malli fit9. Mallin valintaa voidaan tehdä myös toiseen suuntaan samalla funktiolla.

Toinen vaihtoehto automaattiseksi mallinvalintatavaksi on sovittaa kaikki mahdolliset regressiomallit, jotka muuttujista voidaan koostaa, ja sitten valita niistä paras. Tällaisen toiminnallisuuden logistiselle mallille tarjoaa paketti glmulti. Tavanomaisille lineaarisille regressiomalleille voidaan käyttää myös esimerkiksi leaps-pakettia. Muitakin vaihtoehtoja on, ja niistä on kohtuullisen hyvä katsaus esimerkiksi osoitteessa https://rstudio-pubs-static.s3.amazonaws.com/2897_9220b21cfc0c43a396ff9abf122bb351.html. Tehdään seuraavassa analyysi kaikki päävaikutukset ja kahdenväliset interktiot sisältäville malleille:

library(glmulti)
 
glmulti.lm.out <-
    glmulti(died ~ sex + pclass + age + sibsp + parch + embarked, 
            data = d2,
            level = 2,               # No interaction considered
            method = "h",            # Exhaustive approach
            crit = "aic",            # AIC as criteria
            confsetsize = 50,        # Keep 5 best models
            plotty = F, report = T,  # No plot or interim reports
            # glm function with logit link
            fitfunction = "glm", family=binomial)

Analyysi kestää jonkin aikaa, joten on mukava saada väliaikatietoja sovituksen etenemisestä. Siksi tässä on annettu argumentti report = T, koska se tulostaa sovitettujen mallien ja siihen mennessä parhaan löydetyn mallin tulokset säännöllisin väliajoin ruudulle:

After 256050 models:
Best model: died~1+sex+pclass:sex+age:sex+sibsp:pclass+sibsp:age+parch:sex+parch:age+embarked:pclass
Best model: died~1+sex+age+pclass:sex+age:sex+sibsp:pclass+sibsp:age+parch:sex+parch:age+embarked:pclass
Best model: died~1+sex+pclass+pclass:sex+age:sex+sibsp:pclass+sibsp:age+parch:sex+parch:age+embarked:pclass
Crit= 921.212297509378
Mean crit= 922.80260099537

Tulosteesta tulkiten analyysi on käynyt läpi yli 250 000 mallia, ja löytänyt toistaiseksi parhaaksi malliksi ruudulle tulostetun mallin. Mallin AIC-arvo on noin 921.21.

Erikoistapauksia

Katsahdetaan vielä lopuksi eräisiin logistisen regression erikoistapauksiin, ja tulosten esittämiseen visuaalisesti sekä taulukkona että kuvana.

Ehdollinen logistinen malli

Ehdollista logistista regressiota käytetään silloin, jos aineiston tapaukset ja verrokit on jollakin tapaa sovitettu toisiinsa. Yleensä menetelmää käytetään esimerkiksi tapauksissa, joissa on ensin poimittu tapaukset, joille sitten etsitään verrokit, jotka ovat ominaisuuksiltaan varsin samanlaisia kuin tapaukset. Jos tapaus on vaikkapa 40-vuotias nainen, hänelle valittaisiin verrokiksi samanikäinen nainen. Tällä pyritään vähentämään tutkimuksen harhaa. Englanninkieliseltä nimeltään tällainen tutkimusasetelma on matched case-control study.

Kuvatunlainen koeasetelma voitaisiin generoida myös Titanic-aineistosta, mutta käytetään seuraavassa esimerkissä valmista aineistoa. Aineisto löytyy nimellä backpain paketista HSAUR2, ja tutkimuksessa on selvitetty autoilun assosioitumista alaselkäkipuun. Varsinainen analysointi tapahtuu Coxin regressiolla (paketissa survival), sillä ehdollisen logistisen mallin uskottavuusfunktio vastaa Coxin mallin uskottavuusfunktiota tietyllä struktuurilla.

Tällaisen mallin sovittaminen onnistuu esimerkiksi seuraavasti apufunktiolla clogit():

library(HSAUR2)
library(survival)
 
data(backpain)
 
painfit <- clogit(I(status == "case") ~ driver + suburban + strata(ID), data = backpain)
 
Call:
coxph(formula = Surv(rep(1, 434L), I(status == "case")) ~ driver + 
    suburban + strata(ID), data = backpain, method = "exact")
 
  n= 434, number of events= 217 
 
              coef exp(coef) se(coef)     z Pr(>|z|)  
driveryes   0.6579    1.9307   0.2940 2.238   0.0252 *
suburbanyes 0.2555    1.2911   0.2258 1.131   0.2580  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
            exp(coef) exp(-coef) lower .95 upper .95
driveryes       1.931     0.5180    1.0851     3.435
suburbanyes     1.291     0.7746    0.8293     2.010
 
Rsquare= 0.022   (max possible= 0.5 )
Likelihood ratio test= 9.55  on 2 df,   p=0.008457
Wald test            = 8.85  on 2 df,   p=0.01195
Score (logrank) test = 9.31  on 2 df,   p=0.0095

Tulosten perusteella autolla ajaminen on alaselkäsäryn riskitekijä, ja esikaupunkialueella asuminenkin nostaa riskiä, joskaa vaikutus ei ole tilastollisesti merkitsevä.

Jos yhteensovittaminen (matching) jätetään analyysissä huomiotta, on tulos hieman erilainen, joskin tässä tapauksessa tulosten tulkinta ei paljoa muuttuisi:

painfit2<-glm(I(as.numeric(status=="case")) ~ driver + suburban, family=binomial, data=backpain)
 
summary(painfit2)
 
Call:
glm(formula = I(as.numeric(status == "case")) ~ driver + suburban, 
    family = binomial, data = backpain)
 
Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.26431  -1.17576   0.07019   1.09292   1.42033  
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5550     0.2257  -2.460   0.0139 *
driveryes     0.5511     0.2671   2.063   0.0391 *
suburbanyes   0.2059     0.2097   0.982   0.3261  
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 601.65  on 433  degrees of freedom
Residual deviance: 593.61  on 431  degrees of freedom
AIC: 599.61
 
Number of Fisher Scoring iterations: 4

Firthin harhattomampi menetelmä

Tavanomainen logistinen regressio ajautuu ongelmiin erityisesti tilanteessa, jossa aineistossa on jonkin selittävän muuttujan suhteen selkeä separaatio. Separaatio tarkoittaa tilannetta, jossa jollakin (kategorisen) selittävän muuttujan arvolla tai selittäjien yhdistelmällä esiintyy vain toista vastemuuttujan arvoa.

Separaation aiheuttaman ongelman voi ratkaista käyttämällä eksaktia logistista regressiota (R:n laajennuspaketti elrm) tai Firthin muokattua uskottavuusfunktiota (Firth, 1993, Kosmidis et al., 2010) mallin parametrien estimointiin (R:n laajennuspaketit logistf, brglm, mbest). Jos ei ole tarpeen sovittaa aineistoon monimuuttujamallia, voi Fisherin testikin tulla kyseeseen.

Seuraavassa simuloidaan pieni aineisto, jossa on separaatio-ongelma. Sovitetaan aineistoon tavanomainen GLM-malli, ja sama malli kolmella eri funktiolla, jotka kaikki käyttävät Firthin uskottavuusfunktiota. Lisäksi verrataan logististen mallien antamia tuloksia Fisherin testin tulokseen. Poikkeuksena tavanomaiseen esitysmuotoon, ajetut koodirivit alkavat tässä ”>” -merkillä.

> # Simuloidaan aineisto
> y<-c(rep(0, 10), rep(1,10))
> x1<-c(rep(0, 5), rep(1,5), rep(0,10))
> set.seed(123)
> x2<-sample(c(0,1), 20, replace=T)
 
> # Taulukointi
> tab1<-table(y, x1)
> tab2<-table(y, x2)
 
> # Y:n ja X1:n nelikentässä on tyhjiä soluja
> print(tab1)
   x1
y    0  1
  0  5  5
  1 10  0
 
> print(tab2)
   x2
y   0 1
  0 4 6
  1 5 5
 
> # Sovitetaan tavallisella GLM-menetelmällä
> fit1<-glm(y~x1, family=binomial)
> fit2<-glm(y~x2, family=binomial)
> fit3<-glm(y~x1+x2, family=binomial)
 
> # Mallin fit1 antama estimaatti muuttujalle X1 on epärealistisen pieni (käytännössä nolla)
> summary(fit1)
 
Call:
glm(formula = y ~ x1, family = binomial)
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4823  -0.3707   0.4502   0.9005   0.9005  
 
Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    0.6931     0.5477   1.266    0.206
x1           -19.2592  2917.0127  -0.007    0.995
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 27.726  on 19  degrees of freedom
Residual deviance: 19.095  on 18  degrees of freedom
AIC: 23.095
 
Number of Fisher Scoring iterations: 17
 
> summary(fit2)
 
Call:
glm(formula = y ~ x2, family = binomial)
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.2735  -1.1010  -0.0084   1.1271   1.2557  
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.2231     0.6708   0.333    0.739
x2           -0.4055     0.9037  -0.449    0.654
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 27.726  on 19  degrees of freedom
Residual deviance: 27.524  on 18  degrees of freedom
AIC: 31.524
 
Number of Fisher Scoring iterations: 3
 
> summary(fit3)
 
Call:
glm(formula = y ~ x1 + x2, family = binomial)
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.5829  -0.3503   0.4101   0.8576   0.9695  
 
Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    0.9163     0.8367   1.095    0.273
x1           -19.2500  2903.5527  -0.007    0.995
x2            -0.4055     1.1106  -0.365    0.715
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 27.726  on 19  degrees of freedom
Residual deviance: 18.961  on 17  degrees of freedom
AIC: 24.961
 
Number of Fisher Scoring iterations: 17
 
> # Fisherin testi
> # Fisherin testi kärsii samasta ongelmasta kuin GLM-malli yllä:
> # estimoitu OR on käytännössä nolla.
> fisher.test(tab1)
 
        Fisher's Exact Test for Count Data
 
data:  tab1
p-value = 0.03251
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.0000000 0.8365428
sample estimates:
odds ratio 
         0 
 
 
> # Firthin uskottavuusfunktio löytyy kolmesta eri paketista
> library(logistf)
> fit4<-logistf(y~x1)
> summary(fit4)
logistf(formula = y ~ x1)
 
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood Profile Likelihood 
 
                  coef  se(coef) lower 0.95 upper 0.95    Chisq          p
(Intercept)  0.6466272 0.5436237 -0.3519213  1.7589162 1.588980 0.20747191
x1          -3.0445232 1.7069594 -7.9915263 -0.6036788 6.464714 0.01100373
 
Likelihood ratio test=6.464714 on 1 df, p=0.01100373, n=20
Wald test = 3.181209 on 1 df, p = 0.07448959
 
Covariance-Matrix:
           [,1]       [,2]
[1,]  0.2955267 -0.2955267
[2,] -0.2955267  2.9137102
> 
> library(brglm)
> fit5<-brglm(y~x1)
> summary(fit5)
 
Call:
brglm(formula = y ~ x1)
 
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.6466     0.5436   1.189   0.2343  
x1           -3.0445     1.7070  -1.784   0.0745 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 22.606  on 19  degrees of freedom
Residual deviance: 19.973  on 18  degrees of freedom
Penalized deviance: 19.71627 
AIC:  23.973 
 
> 
> library(mbest)
> fit6<-glm(y~x1, family=binomial(), method="firthglm.fit")
> summary(fit6)
 
Call:
glm(formula = y ~ x1, family = binomial(), method = "firthglm.fit")
 
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.4616  -0.6790   0.2498   0.9177   0.9177  
 
Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.6471     0.5437   1.190   0.2340  
x1           -3.0404     1.7041  -1.784   0.0744 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
(Dispersion parameter for binomial family taken to be 1)
 
    Null deviance: 27.726  on 19  degrees of freedom
Residual deviance: 19.976  on 18  degrees of freedom
AIC: 23.976
 
Number of Fisher Scoring iterations: 3

Tässä demoaineistossa olevaa separaatio-ongelmaa ei pysty ratkaisemaan yksinkertaisella Fisherin testillä, vaan pitää käyttää esimerkiksi Firthin uskottavuusfunktiota. Heinze et al. (2002) argumentoi, että Firthin uskottavuusfunktio tarjoaa soveltuvan ratkaisun myös tähän ongelmaan.
Firthin uskottavuusfunktio on implementoitu kolmeen eri pakettiin, joista mbest ja brglm antavat keskenään saman tuloksen niin parametrien estimaattien kuin niiden p-arvojen osaltakin (neljänteen desimaaliin asti). Paketti logistf sen sijaan poikkeaa kahdesta muusta p-arvojen osalta, sillä ne ovat kahden muun tulokseen verrattuna selkeästi pienempiä, vaikka parametrien estimaatit ovatkin samanlaisia. Erot näyttävät johtuvan siitä, että sekä mbest- että brglm-pakettien yhteenvedossa p-arvot on laskettu t-testisuureen perusteella (2 * pnorm(-abs(tvalue))) kun taas logistf käyttää Khiin neliötestiä (1 – pchisq(2 * (fit.full$loglik – fit.i$loglik), 1)).

Harrellin lrm-funktio

Frank Harrellin tuottamassa rms-paketissa on funktio lrm(), jolla voidaan sovittaa varsin monipuolisesti logistisia malleja. Yleensä on tapana sovittaa splineihin perustuvia additiivisia malleja (GAM) esimerkiksi mgcv-paketin kautta, mutta lrm-funktiolla niidenkin käyttäminen onnistuu.

Perusmuodossaa käytettynä funktio antaa samantyyppisen tulosteen kuin glm()-funktiokin, mutta mukana on monenlaisia indeksejä valmiiksi laskettuina:

fit.lrm <- lrm(died ~ pclass + sex + age + sibsp + parch, data=d)
summary(fit.lrm)
 
Frequencies of Missing Values Due to Each Variable
  died pclass    sex    age  sibsp  parch 
     1      1      1    264      1      1 
 
Logistic Regression Model
 
 lrm(formula = died ~ pclass + sex + age + sibsp + parch, data = d)
 
 
                       Model Likelihood     Discrimination    Rank Discrim.    
                          Ratio Test           Indexes           Indexes       
 Obs          1046    LR chi2     443.40    R2       0.466    C       0.846    
  0            427    d.f.             5    g        1.900    Dxy     0.691    
  1            619    Pr(> chi2) <0.0001    gr       6.685    gamma   0.694    
 max |deriv| 6e-09                          gp       0.342    tau-a   0.334    
                                            Brier    0.148                     
 
           Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept -4.9568 0.4341 -11.42 <0.0001 
 pclass     1.1582 0.1129  10.26 <0.0001 
 sex=male   2.5533 0.1732  14.74 <0.0001 
 age        0.0386 0.0065   5.89 <0.0001 
 sibsp      0.3449 0.1051   3.28 0.0010  
 parch     -0.0768 0.1000  -0.77 0.4422

Harrellin funktion käyttäminen mahdollistaa myös siloitettujen termien käytön mallissa:

fit.lrm <- lrm(died ~ pclass + sex + rcs(age, 3) + sibsp + parch, data=d)
 
Frequencies of Missing Values Due to Each Variable
  died pclass    sex    age  sibsp  parch 
     1      1      1    264      1      1 
 
Logistic Regression Model
 
 lrm(formula = died ~ pclass + sex + rcs(age, 4) + sibsp + parch, 
     data = d)
 
 
                       Model Likelihood     Discrimination    Rank Discrim.    
                          Ratio Test           Indexes           Indexes       
 Obs          1046    LR chi2     453.61    R2       0.475    C       0.849    
  0            427    d.f.             7    g        1.926    Dxy     0.699    
  1            619    Pr(> chi2) <0.0001    gr       6.862    gamma   0.702    
 max |deriv| 8e-09                          gp       0.345    tau-a   0.338    
                                            Brier    0.145                     
 
           Coef    S.E.   Wald Z Pr(>|Z|)
 Intercept -6.0092 0.5649 -10.64 <0.0001 
 pclass     1.1173 0.1129   9.90 <0.0001 
 sex=male   2.6049 0.1760  14.80 <0.0001 
 age        0.1156 0.0258   4.48 <0.0001 
 age'      -0.2106 0.0759  -2.78 0.0055  
 age''      0.8132 0.3384   2.40 0.0163  
 sibsp      0.4228 0.1097   3.86 0.0001  
 parch      0.0115 0.1044   0.11 0.9126

Lisäksi lrm-objektille voidaan tehdä helposti mm. bootstrapping-menetelmään perustuva validointi. Tällöin tulosobjektiin halutaan mukaan myös aineisto (argumentit x=T, y=T):

fit.lrm <- lrm(died ~ pclass + sex + age + sibsp + parch, data=d, x=T, y=T)
validate(fit.lrm)
          index.orig training   test optimism index.corrected  n
Dxy           0.6915   0.6874 0.6880  -0.0006          0.6921 40
R2            0.4660   0.4650 0.4619   0.0031          0.4629 40
Intercept     0.0000   0.0000 0.0087  -0.0087          0.0087 40
Slope         1.0000   1.0000 0.9965   0.0035          0.9965 40
Emax          0.0000   0.0000 0.0025   0.0025          0.0025 40
D             0.4229   0.4223 0.4183   0.0040          0.4189 40
U            -0.0019  -0.0019 0.0001  -0.0020          0.0001 40
Q             0.4249   0.4242 0.4182   0.0060          0.4188 40
B             0.1481   0.1483 0.1493  -0.0011          0.1491 40
g             1.8998   1.9000 1.8825   0.0174          1.8824 40
gp            0.3417   0.3408 0.3398   0.0010          0.3407 40

Funktion lrm käyttämisessä on useita piirteitä, joita ei ole helppo toteuttaa muilla funktioilla. Yksi on nomogrammin luonti mallin perusteella. Ennen mallin sovittamista pitää ajaa funktio datadist, joka päättelee asioita datasta valmiiksi. Tämän jälkeen sovitetaan malli ja piirretään nomogrammi. Syntyvä kuva on esitetty koodiblokin alapuolella.

ddist <- datadist(d); options(datadist='ddist')
fit.lrm <- lrm(died ~ pclass + sex + age + sibsp + parch, data=d, x=T, y=T)
plot(nomogram(fit.lrm))

nomogrammi

Nomogrammia tulkitaan siten, että eri muuttujien arvot yhdistetään viivaimella, ja pisteasteikolta luetaan yhdistelmän saama pistemäärä, joka tässä tapauksessa kuvaa kuolemanvaaraa Titanicin onnettomuudessa. Nomogrammista näkee helposti, että esimerkiksi kakkosluokan miesmatkustajalla on paljon pienempi riski menehtyä kuin kolmosluokan miesmatkustajalla, sillä näiden yhdistelmien saamat pistemäärät Points-asteikolla ovat varsin eri suuruiset.

Paketin toiminnallisuuksiin kannattaa ehdottomasti tutustua tarkemmin muun muassa Frank Harrelin kirjasta Regression Modeling Strategies.

Tulosten visuaalinen esittäminen

Edellä katsahdettiin jo lyhyesti nomogrammimuotoiseen tulosten esittämiseen. Tavanomainen tapa esittää tulokset on toki taulukko. Yllä on jo katsahdettu taulukoidenkin tuottamiseen, mutta paketti stragazer tarjoaa tähänkin tarkoitukseen erinäisiä apukeinoja. Stargazer mahdollistaa kauniiden HTML- ja LaTeX-taulukoiden tuottamisen, mutta se osaa myös piirtää kauniihkoja tesktimuotoisia taulukoita. Taulukoissa voidaan esittää yksi tai useampia malleja, joiden niiden antamia tuloksia voidaan helposti vertailla:

library(stargazer)
stargazer(fit1, fit2, fit3, fit5, type="text")
=========================================================
                            Dependent variable:          
                  ---------------------------------------
                                   died                  
                     (1)       (2)       (3)       (4)   
---------------------------------------------------------
sexmale           2.460***            2.524***  2.491*** 
                   (0.152)             (0.163)   (0.166) 
                                                         
pclass                      0.796***  0.855***  1.131*** 
                             (0.080)   (0.095)   (0.112) 
                                                         
age                                             0.034*** 
                                                 (0.006) 
                                                         
Constant          -1.106*** -1.350*** -2.995*** -4.583***
                   (0.118)   (0.183)   (0.256)   (0.406) 
                                                         
---------------------------------------------------------
Observations        1,044     1,044     1,044     1,044  
Log Likelihood    -550.434  -652.038  -506.701  -491.173 
Akaike Inf. Crit. 1,104.868 1,308.076 1,019.403  990.345 
=========================================================
Note:                         *p<0.1; **p<0.05; ***p<0.01

Toinen näppärä tapa esittää tulokset on piirtää ne kuvaksi. Kuvassa on esitetty estimoitu mallin parametrin arvo ja sen 95 % luottamusväli. Kuvan tuottaminen onnistuu automaattisesti arm-paketin funktiolla coefplot():

library(arm)
 
par(mar=c(5,7,2,2))
coefplot(fit8, xlim=c(-3,3))

Muodostuva kuva on seuraavanlainen:

coefplot

Lisäksi logistisen malli käyrä voidaan esittää kuvana siten, että siihen on merkitty havaintopisteet tai vähintäänkin mallin ennusteet selittäjein eri arvoilla. Tämän mahdollistaa esimerkiksi paketti effects. Valitettavasti paketti on hieman nirso syötteensä suhteen, joten sovitetaan ensin malli 8 uudestaan, kun muuttujat sex ja embarked on ensin muutettu faktoreiksi. Piirretään tämän jälkeen kuva mallin ennusteista:

library(effects)
 
d2$sex<-factor(d2$sex)
d2$embarked<-factor(d2$embarked)
fit8 <- glm(died ~ sex + pclass + age + sibsp + parch + embarked, data=d2, family=binomial)
plot(allEffects(fit8))

Muodostuvassa kuvassa on yksi kuva kutakin selittävää muuttujaa kohden:

effects

Yhteenveto

Tässä artikkelissa on esitelty logistisen regression sovittamiseen liityviä käytänteitä R-ohjelmistolla. Mallin sovittamiseen liittyy tiettyä erikoispiirteitä, jotka on yllättävän helppo sivuuttaa.

Eräs ärsyttävimmistä piirteistä on R:n tapa ilmoittaa glm()-funktion tulosteessa ns. Waldin menetelmällä muodostettujen testien tulokset, mutta mallin tuloksille automaattisesti luottamusvälit laskeva funktio puolestaan ilmoittaa automaattisesti profile likelihoodiin perustuvat luottamusvälit.

Harrellin tuottamassa rms-paketissa on monia kiinnostvia ominaisuuksia, mutta sen käyttö vaatii rinnalleen saman nimisen kirjan, jotta kaikista nyansseista saa kiinni. Lisäksi se ei ole yhteensopiva R:n perus-glm:n kanssa.

Nykyisin artikkeleissa on hieman yleistynyt käytäntö, jossa mallin tuloksia ei välttämättä enää esitetä taulukkona, vaan kuvana. Käytäntö tukee kyllä tulosten silmäilyä, mutta estää tulosten käyttämisen esimerkiksi meta-analyyseissä. Kuvia kannattaa siis viljellä hienovaraisesti, ja mahdollisuuksien mukaan sisällyttää myös taulukkomuotoiset esitykset vaikkapa supplementary-informaatioksi, jos tieteellinen julkaisu sellaisen käytännön sallii.

Tags: , , ,