R-ohjelmointi.org

Tilastotieteellistä ohjelmointia R-kielellä

Geneettinen epidemiologia – OR vastaan ennustearvo

Geneettinen epidemiologia kattaa suuren joukon erilaisia tutkimusasetelmia, joiden tarkoituksena on löytää tiettyyn tautitilaan liittyviä geneettisiä piirteitä. Eräs tyypillinen asetelma on tapaus-verrokkitutkimus, jossa tutkitaan yhden emäksen muutoksia genomissa. Varsin tyypillisesti tieteellisissä julkaisuissa esitetään tällaisissa tapauksissa geneettisen ominaisuuden odds ratio (vetokerroin), muttei välttämättä sitä, miten paljon lisäarvoa ominaisuus tuo jo tunnettujen taudin riskitekijöiden lisäksi. Lisäarvo olisi kuitenkin kohtullisella vaivalla selvitettävissä. Katsahdetaanpa seuraavassa R:n suomiin mahdollisuuksiin.

R:n paketti gap tarjoaa sopivan demoaineiston nimellä nep499. Aineistossa on terveitä verrokkeja ja sairaita Alzheimerin tautipotilaita. Iän tiedetään altistavan Alzheimerin taudille, ja sukupuolten välillä on eroa taudin ilmaantuvuudessa. Aineistosta löytyvän APO-E -geenin alleelimuoto 4 on aiemmin kuvattu taudin riskitekijäksi. Gap-paketin lisäksi tarvitaan paketit caret, epicalc ja boot erinäisten mallin validointiin tarvittavien laskutoimitusten tekemiseksi. Sama asia voidaan hoitaa myös rms-paketilla, mistä tämän postauksen lopussa lisää.

Ladataan ensin paketit ja aineisto:

require(caret)
require(epicalc)
library(boot)
require(gap)
 
data(nep499)

Sovitetaan sitten joukko erilaisia regressiomalleja:

fit1<-glm(status~age, data=nep499, family=binomial())
fit2<-glm(status~sex, data=nep499, family=binomial())
fit3<-glm(status~apoe4, data=nep499, family=binomial())
fit4<-glm(status~age+sex, data=nep499, family=binomial())
fit5<-glm(status~age*sex, data=nep499, family=binomial())
fit6<-glm(status~age+sex+apoe4, data=nep499, family=binomial())
fit7<-glm(status~age*sex+apoe4, data=nep499, family=binomial())
fit8<-glm(status~age*apoe4*sex, data=nep499, family=binomial())

Mallien sopivuutta aineistoon voidaan mitata vaikkapa AIC-arvoilla. Niillä voidaan myös tehdä mallin valintaa:

aic<-AIC(fit1,fit2,fit3,fit4,fit5,fit6,fit7,fit8)
aic$delta<-aic$AIC-min(aic$AIC)
aic
#    df      AIC      delta
#fit1  2 675.6579 15.7878400
#fit2  2 693.0749 33.2048456
#fit3  2 689.7413 29.8712294
#fit4  3 672.8274 12.9573665
#fit5  4 660.6442  0.7742166
#fit6  4 671.9815 12.1114944
#fit7  5 659.8700  0.0000000
#fit8  8 663.3539  3.4839207

AIC-arvoilla mitattuna paras malli on siis fit7, mutta fit5 on käytännössä miltei yhtä hyvä, koska AIC-arvot ovat varsin lähellä toisiaan.

Perinteinen tapa mitata mallin antamien ennusteiden osuvuutta on muuttaa logistisen regression antamat todennäköisyysestimaatit 0/1-asteikolle, ja verrata näitä sitten alkuperäisiin havaintoihin. Näin mallille voidaan laskea muun muassa sensitiivisyys ja spesifisyys. Caret-paketin funktiot avustavat tässä:

# fit5
sensitivity(table(pred=as.numeric(predict(fit5, type="response")>(sum(nep499$status)/nrow(nep499))), obs=fit5$model$status))
#[1] 0.4876033
specificity(table(pred=as.numeric(predict(fit5, type="response")>(sum(nep499$status)/nrow(nep499))), obs=fit5$model$status))
#[1] 0.688716
 
# fit7
sensitivity(table(pred=as.numeric(predict(fit7, type="response")>(sum(nep499$status)/nrow(nep499))), obs=fit7$model$status))
#[1] 0.5661157
specificity(table(pred=as.numeric(predict(fit7, type="response")>(sum(nep499$status)/nrow(nep499))), obs=fit7$model$status))
#[1] 0.6692607

Luokittelumatriisit näyttävät seuraavilta:

# fit5
(table(pred=as.numeric(predict(fit5, type="response")>(sum(nep499$status)/nrow(nep499))), obs=fit5$model$status))
#    obs
#pred   0   1
#   0 118  80
#   1 124 177
 
# fit7
(table(pred=as.numeric(predict(fit7, type="response")>(sum(nep499$status)/nrow(nep499))), obs=fit7$model$status))
    obs
#pred   0   1
#   0 137  85
#   1 105 172

Tällä perusteella näyttää siltä, että apo-e4:n lisääminen malliin parantaa mallin sensitiivisyyttä, siis sitä, että terveet henkilöt tulevat ennustettua oikein. Toisaalta, samalla sairaiden potilaiden ennustamistarkkuus heikkenee, mikä ei ole ollenkaan toivottavaa.

Mallien tulokset näyttävät seuraavilta:

round(data.frame(OR=exp(coef(fit5)), CI=exp(confint(fit5))), 3)
#Waiting for profiling to be done...
#                     OR CI.2.5..    CI.97.5..
#(Intercept) 4178287.119 6318.196 4.168488e+09
#age               0.815    0.745 8.870000e-01
#sex               0.001    0.000 5.300000e-02
#age:sex           1.096    1.044 1.153000e+00
 
round(data.frame(OR=exp(coef(fit7)), CI=exp(confint(fit7))), 3)
#Waiting for profiling to be done...
#                     OR CI.2.5..    CI.97.5..
#(Intercept) 3338446.936 4930.103 3.398887e+09
#age               0.817    0.747 8.890000e-01
#sex               0.001    0.000 5.200000e-02
#apoe4             1.479    0.934 2.365000e+00
#age:sex           1.096    1.044 1.153000e+00

Muuttuja apoe4 ei ole mallissa tilastollisesti merkitsevä, mutta sen odds ratio (OR) on n. 1.5, joten se voi olla heikko Alzheimerin taudin riskitekijä. Tätä ei useinkaan myöskään huomioida artikkelien johtopäätöksissä: vaikka jokin vaikutuksen estimaatti ei olekaan tilastollisesti merkitsevä, voi se silti olla merkityksellinen, muttei tietysti välttämättä. Tässä APO-E4:n OR on aika pieni, joten sen merkitys lienee kokonaisuudessaan varsin vähäinen.

Sensitiivisyys ja spesifisyys eivät kuitenkaan välttämättä ole parhaat mahdolliset mittarit mallille, koska ne riippuvat sitä, mistä logistisen mallin antaman todennäköisyyden arvot on katkaistu dikotomisointia varten. Mallin ennustevoimaa voidaan nimittäin arvioida muillakin mittareilla. Eräs tällainen on Receiver Operating Curve:n alle jäävän alueen pinta-ala (auc). Suurempi auc-arvo kielii korkeammasta ennustevoimasta. Tämä on mahdollista arvioida mm. epicalc-paketin funktiolla:

lroc(fit5)$auc
#[1] 0.6529247
lroc(fit7)$auc
#[1] 0.6634965

Mallin auc-arvo on vain marginaalisesti parempi, jos siihen on lisätty APO-E4 muuttujaksi.

Kaiken kaikkiaan APO-E4:n vaikutus ei tässä aineistossa ole kovin merkittävä, eikä se kenties juurikaan tarkenna mallin antamia ennusteita verrattuna malliin, jossa ovat vain ikä ja sukupuoli. Merkittävää on mm. tehdä päätös, onko tärkeämpää ennustaa oikein terveet vai sairaat, ja valita malli sen mukaan.

Edellä esitetyn voi tehdä yksinkertaisesti myös rms-paketilla seuraavasti:

lrm1<-lrm(formula=status~age*sex, data=nep499, model=TRUE, x=TRUE, y=TRUE)
lrm2<-lrm(formula=status~age*sex+apoe4, data=nep499, model=TRUE, x=TRUE, y=TRUE)
lrm1
 
#Logistic Regression Model
#
#lrm(formula = status ~ age * sex, data = nep499, model = TRUE, 
#    x = TRUE, y = TRUE)
#
#                      Model Likelihood     Discrimination    Rank Discrim.    
#                         Ratio Test            Indexes          Indexes       
#Obs           499    LR chi2      38.67    R2       0.099    C       0.653    
# 0            242    d.f.             3    g        0.638    Dxy     0.306    
# 1            257    Pr(> chi2) <0.0001    gr       1.893    gamma   0.311    
#max |deriv| 9e-09                          gp       0.146    tau-a   0.153    
#                                           Brier    0.231                     
#
#          Coef    S.E.   Wald Z Pr(>|Z|)
#Intercept 15.2454 3.4084  4.47  <0.0001 
#age       -0.2042 0.0442 -4.62  <0.0001 
#sex       -6.7134 1.9556 -3.43  0.0006  
#age * sex  0.0917 0.0251  3.65  0.0003  
 
validate(lrm1)
#          index.orig training    test optimism index.corrected  n
#Dxy           0.3058   0.3181  0.2938   0.0243          0.2816 40
#R2            0.0994   0.1098  0.0946   0.0153          0.0842 40
#Intercept     0.0000   0.0000 -0.0103   0.0103         -0.0103 40
#Slope         1.0000   1.0000  0.9409   0.0591          0.9409 40
#Emax          0.0000   0.0000  0.0152   0.0152          0.0152 40
#D             0.0755   0.0842  0.0715   0.0126          0.0629 40
#U            -0.0040  -0.0040  0.0011  -0.0051          0.0011 40
#Q             0.0795   0.0882  0.0704   0.0178          0.0617 40
#B             0.2313   0.2288  0.2332  -0.0045          0.2358 40
#g             0.6380   0.6685  0.6115   0.0570          0.5810 40
#gp            0.1459   0.1503  0.1406   0.0097          0.1362 40

Yllä olevassa taulukossa C vastaa ROC-käyrän auc-arvoa, mutta tuloste antaa myös monia muita arvoja (discriminant index). Mallin validointi on yllä hoidettu saman nimisellä funktiolla, joka arvioi myös indeksi optimisimin eli mahdollisen ylisovituksen määrän indeksien arvoissa.

Paketti rms onkin aiemmin esitetty yksinkertaisempi tapa hoitaa eri analyysin vaiheet, koska se tekee monet seikat automaattisesti.


Vastaa

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