R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Logistisen regression ongelmatilanteita ja niiden ratkaisuja

Logistinen regressio on esimerkiksi epidemiologian tilastollinen perustyökalu. Katsahdetaanpa seuraavassa kahteen usein vastaan tulevaan ongelmaan ja niiden mahdollisiin ratkaisuihin.

Ongelma 1: separaatio

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)).

Ongelma 2: todennäköisyyksistä luokiksi

Logistisella regressiolla sovitteesta voidaan ennustaa kullekin havainnolle todennäköisyys, että vastemuuttuja saa arvon 1. Jos ennustetut todennäköisyydet halutaan muuttaa vastemuuttujan arvoiksi, pitää todennäköisyydet jakaa kahteen luokkaan jonkin katkaisupisteen perusteella. Usein katkaisupisteenä käytetään 0.5:ttä tai vastemuuttujan saamien ykkösarvojen suhteellista osuutta aineistossa.

Muitakin vaihtoehtoja katkaisupisteen päättelemiseksi on. Eräs yleisesti käytetty tapa on päätellä katkaisupiste ROC-käyrästä.

Erilaisten menetelmien käyttämien sopivan katkaisupisteen määrittämiseksi on mahdollista automatisoida R:n laajennuspaketeilla PresenseAbsense ja OptimalCutpoints. Seuraavassa koodiesimerkki pakettien käytöstä:

> library(PresenceAbsence)
> library(OptimalCutpoints)
> 
> # OptimalCutpoints: käytettävät metodit
> meths<-c(
+ "ROC01",
+ "MinPvalue",
+ "ObservedPrev",
+ "MeanPrev")
> 
> # OptimalCutpoints:
> # menetelmä vaatii erikseen muotoillun syötedatan (havaitut vastearvot ja mallin ennusteet)
> tmp<-data.frame(observed=y, predicted=predict(fit6, type="response"))
> cp <- optimal.cutpoints(X = "predicted", status = "observed", tag.healthy = 0, methods = meths, data = tmp,  pop.prev = NULL, control = control.cutpoints(), ci.fit = TRUE)
Warning: Sensitivity CI: "Exact" method may not be valid for some values (see Help Manual).
 
Warning: Specificity CI: "Exact" method may not be valid for some values (see Help Manual).
 
Warning: Positive Predictive Value CI: "Exact" method may not be valid for some values (see Help Manual).
 
Warning messages:
1: In chisq.test(tabl) : Chi-squared approximation may be incorrect
2: In chisq.test(tabl) : Chi-squared approximation may be incorrect
> print(cp)
 
Call:
optimal.cutpoints.default(X = "predicted", status = "observed", 
    tag.healthy = 0, methods = meths, data = tmp, pop.prev = NULL, 
    control = control.cutpoints(), ci.fit = TRUE)
 
Optimal cutoffs:
   ROC01 MinPvalue ObservedPrev MeanPrev
1 0.6564    0.6564       0.6564   0.0837
2      -         -            -   0.6564
 
Area under the ROC curve (AUC):  0.75 (0.587, 0.913) 
> 
> # PresenceAbsence
> # menetelmä vaatii erikseen muotoillun syötedatan (ID, havaitut vastearvot ja mallin ennusteet)
> tmp<-data.frame(plotID=1:length(y), Observed=y, predicted=predict(fit6, type="response"))
> optimal.thresholds(tmp)
         Method predicted
1       Default  0.500000
2     Sens=Spec  0.370000
3  MaxSens+Spec  0.370000
4      MaxKappa  0.370000
5        MaxPCC  0.370000
6  PredPrev=Obs  0.370000
7       ObsPrev  0.500000
8      MeanProb  0.513186
9    MinROCdist  0.370000
10      ReqSens  0.650000
11      ReqSpec  0.660000
12         Cost  0.370000
Warning messages:
1: In optimal.thresholds(tmp) : req.sens defaults to 0.85
2: In optimal.thresholds(tmp) : req.spec defaults to 0.85
3: In optimal.thresholds(tmp) : costs assumed to be equal

Mallin antaa havainnoille vain kahdenlaisia ennusteen arvoja: 0.08368114 ja 0.65635428. Käytännössä kaikilla menetelmillä määritetyt todennäköisyyksien dikotomisointipisteet sijoittuvat noiden arvojen väliin, joten tässä käytettävällä menetelmmällä ei ole väliä. Useimmiten tulokset kuitenkin vaihtelevat hieman menetelmästä riippuen. Esimerkiksei PresenceAbsence-paketin tuloksissa MinROCdist vastaa ROS-käyrämenetelmää, ja ObsPrev tilannetta, jossa dikotomisointikohtana käytetään vastemuuttujan suurempien arvojen suhteellista osuutta.

Mallin tulokset voidaan muutta binäärisiksi ja esittää vaikkapa confusion matriisina seuraavasti:

> table(obs=y, pred=as.numeric(predict(fit6, type="response")>0.37))
 
   pred
obs  0  1
  0  5  5
  1  0 10

Yhteenveto

Edellä esitetyt R:n laajennuspaketit lisäävät logistista regressiota käyttävien työkalupakkiin kaksi tärkeää työkalua: mahdollisuuden sovittaa logistinen malli myös aineistoon, jossa esiintyy separaatiota ja mahdollisuuden määrittää havainnoille ennustettujen todennäköisyyksien katkaisukohdat useilla eri menetelmillä automaattisesti.

Kirjallisuus

Firth D (1993) Bias Reduction of Maximum Likelihood Estimates, Biometrica, 80, 27-38.
Heinze G, Schemper M (2002) A solution to the problem of separation in logistic regression, Statist. Med., 21, 2409-2419.
Kosmidis I, Firth D (2010) A generic algorithm for reducing bias in parametric estimation, Electr. J. of Statistics, 4, 1097-112.

Tags: