6.6 Applications

datataulukon luominen voidaan tehdä jäljempänä esitetyllä tavalla komentomatriisilla () luodusta taulukosta:

> tab < – as.taulukko (matriisi(c (815,115,208,327), nrow = 2, byrow = TRUE,

dimnames = list (EST = c(”+”,”−”),

CAD = c(”+”,”-”))))

> tab

on huomattava, että taulukko järjestettiin uudelleen vastaamaan tavanomaisia merkintöjä, ja ”positiivisia” tapahtumia oli taulukon ensimmäisellä rivillä (tyypillisesti altistus) ja ensimmäisessä sarakkeessa (tyypillisesti tauti).

tällöin olisi mahdollista laskea kaikki vaaditut suureet käyttämällä hyvin harvoja aritmeettisia operaatioita R: llä. Epir-paketti sisältää kuitenkin kaikki 2 × 2-taulukoiden epidemiologisiin kysymyksiin vastaamiseen tarvittavat välineet. Näin ollen komento epi.testeistä () saadaan esiintyvyys, herkkyys/spesifisyys sekä positiiviset ja negatiiviset ennustearvot:

> kirjasto (epiR)

> epi.testit (tab)

käyrän alle jäävän alueen osalta voidaan käyttää myös ulkoista käskyä, roc.päässä.taulukko (), paketissa epicalc:

> library(epicalc)

> roc.päässä.taulukossa (tab, graph = FALSE)

datan tuominen hdis.dat saavutetaan seuraavasti:

> bp < – read.pöytä (”hdis.dat”, header = TRUE)

> str(bp)

kuten voidaan todentaa, ei ole merkintää, joka liittyisi merkittävien muuttujien tasoihin (bpress verenpaineelle, chol kolesteroliarvolle). Tarrojen luomiseen ja liittämiseen voidaan lisätä seuraavat komennot:

> blab < – c(”<117”,”117-126”,”127-136”,”137-146”,

”147-156”,”157-166”,”167-186”,”>186”)

> clab < – c(”<200”,”200-209”,”210-219”,”220-244”,

”245-259”,”260-284”,”>284”)

> bp$bpress < – factor(bp$bpress, labels = blab)

> bp$chol < – factor(BP$chol, labels = clab)

viimeinen komento muuntaa alkuperäiset muuttujat kvalitatiivisiksi muuttujiksi ja liittää niiden modaliteetteihin blabin ja clabin määrittelemät tarrat. Sen varmistamiseksi, että tietokanta on nyt halutussa muodossa, seuraavat komennot ovat hyödyllisiä:

> str (bp)

> Yhteenveto (bp)

:

> round(xtabs(hdis/total bpress + Chol, data = BP), 2)

koska aiomme keskittyä vain verenpaineen ja sydäninfarktin väliseen yhteyteen, on tarpeen koota kolesteroliarvoja koskevat tiedot. Toisin sanoen on tarpeen kerätä arvot jokaiselle verenpainetasolle riippumatta kolesterolitasoista. Myös muuttujan bpress tasot tulisi nimetä uudelleen luokkavälien keskusten avulla:

> blab2 < – c(111,5,121,5,131.5,141.5,151.5,161.5,176.5,191.5)

> BP$bpress < – rep(blab2, each = 7)

> dfrm < – aggregaatti(bp

lista(bpress = BP), summa)

se on viimeinen komento, aggregaatti (), joka mahdollistaa aineiston yhdistämisen: se summaa kaikki muuttujan chol arvot (summa), jotka eivät sisälly analyysiä varten säilytettävien muuttujien luetteloon. Samalla tulokset tallennetaan uuteen dfrm-nimiseen tietokantaan. Huomaa, että tässä tapauksessa emme käytä kaava-merkintää vaan vaihtoehtoista syntaksia aggregaatille() (vastemuuttuja, jota seuraa luettelo, joka koostuu yhdestä tai useammasta selittävästä muuttujasta). Yleiskuva aggregoiduista tiedoista saadaan antamalla syötteeksi head ():

> head(dfrm, 5)

kun suhde p on saatavilla, sen arvo asteikolla, jonka yksiköt ovat logitteja, saadaan relaatiolokista(p/(1 − p)). Siitä voidaan johtaa seuraavat komennot muuntamaan mittasuhteet, jotka lasketaan kuten hdis / total logit-yksiköinä:

> logit < – function(x) log(x/(1-x))

> dfrm$prop < – dfrm$hdis/dfrm$Total

> dfrm$logit < – logit(dfrm$hdis/dfrm$Total)

on huomattava, että muunnosarvoille X on määritelty pieni funktio, jonka tässä oletetaan olevan mittasuhteet, sen ekvivalenttisessa logissa(x/(1 − x)). Yhtä lailla voisi kirjoittaa:

> log((dfrm$hdis/dfrm$total)/(1-dfrm$hdis/dfrm$total))

lasketut arvot on lisätty suoraan tietokehykseen nimillä prop ja logit.

logistinen regressiomalli on kirjoitettu seuraavasti:

> glm(cbind(hdis,total-hdis) bpress,data = dfrm,family = binominen)

käytetty formulaatio cbind(hdis, total-hdis) bpress ottaa huomioon sen, että käytämme ryhmitettyä dataa emmekä yksittäisiä vastauksia. Komento glm (), jossa valitsinperhe = binomiaali, vastaa logistista regressiota, jossa ei mennä sen kummemmin yksityiskohtiin, käytetään oletusarvoisesti logit-funktiota kanonisena linkkinä.

regressiokertoimet voidaan esittää käyttämällä yhteenvetoa ():

> summary (glm (cbind (hdis, total-hdis) bpress, data = dfrm,

family = binominen))

edellinen tulos sisältää olennaisen tiedon, jotta voidaan vastata kysymykseen verenpaineen ja sydänkohtauksen todennäköisyyden välisen yhteyden tilastollisesta merkitsevyydestä: regressiokerroin (log-odds-asteikolla) on 0, 024 ja se on merkittävä tavanomaisella 5%: n raja-arvolla (KS.sarake Pr(>|z|)).

sydänkohtauksen todennäköisyys tarkasteltavien verenpaineen eritasojen mukaan saadaan seuraavasti:

> glm.res < – glm(cbind(hdis, total-hdis) bpress, data = dfrm,

family = binominen)

> predict(glm.res)

on huomattava, että R: n tuottamat välitulokset on tallennettu nimellä glm.res ennen komennon predict () käyttöä. R: n tuottamat ennusteet ilmaistaan logit-muodossa, ja havaittuja ja ennustettuja logitteja on mahdollista verrata keskenään:

> cbind(dfrm, logit.predit = ennustavat (glm.res))

lähes kaikki tekijät ovat käytettävissä kuvaamaan todettujen ja ennustettujen sydäninfarktien osuuksia graafisesti verenpaineen mukaan. Mallin ennustukset, jotka on ilmaistu mittasuhteina eikä logodds-asteikolla, puuttuvat. Siksi voi olla suotavaa piirtää ennustekäyrä eli sydänkohtauksen todennäköisyys verenpaineen mukaan rajoittamatta jälkimmäistä muuttujan bpress osalta havaittuihin kahdeksaan arvoon. Tässä seuraa mahdollinen ratkaisu:

> dfrm$prop.predit < – predict(glm.res, type = ”response”)

> F < – function(x)

1/(1 + exp(-(coef(glm.res) + coef(glm.res)*x)))

> xyplot(hdis/total ˜ bpress, data = dfrm, aspect = 1.2, cex = .8,

xlab = ”Blood pressure”,ylab = ”Infarction likelihood”,

panel = function(x, y, …) {

panel.xyplot(x, y, col = ”gray30”, pch = 19, …)

panel.curve(f, lty = 3, col = ”gray70”)

panel.points(x, dfrm$prop.predit, col = ”gray70”, …)

})

koska (yksittäisiä) raakatietoja ei ole saatavilla, on tarpeen käyttää suoraan ongelmalausekkeessa annettua frekvenssitaulukkoa:

> alcohol < – matrix(c(666,104,109,96), nr = 2,

dimnames = list(c(”Control”,”Case”),

c(”<> = 80”)))

> alcohol

Regarding the odds ratio, the command oddsratio() from the package vcd will be used:

> library(vcd)

> oddsratio(alcohol, log = FALSE)

vaihtoehto log = FALSE varmistaa, että palautettu tulos vastaa oikein kerroinsuhdetta eikä kerrossuhteen logia. Asymptoottisen luottamusvälin saamiseksi käytetään seuraavaa:

> confint(oddsratio(alcohol, log = FALSE))

vastaavasti tiivistelmää(oddsratio(alcohol)) voidaan käyttää hypoteesitestin tekemiseen log odds ratiolla (H0 : log(θ) = 0).

komentoehdotus.testillä () voidaan testata hypoteesia, jonka mukaan päivittäiskulutuksen saaneiden osuus on ≥ 80 g ja sama sekä tapauksissa että kontrolleissa, mikä osoittaa ongelmalausekkeessa annetusta ristitaulukosta havaitut arvot:

> prop.testi(c(96,109), c (200,775), correct = FALSE)

Tämä testi vastaa täsmälleen Z-testiä, jolla testataan kahden tiedoista estimoidun osuuden välinen ero (jos jatkuvuuskorjausta ei käytetä).

logistisen regressiomallin osalta ennustetaulukko on muunnettava datataulukoksi, jossa kaksi kvalitatiivista muuttujaa (sairaus ja altistus) esitetään oikein, toisin sanoen tietokehyksellä. Tämä voidaan tehdä käyttämällä reshape2-paketin komentoa sula (), joka mahdollistaa datan muuttamisen taulukkomuodossa ”pitkään” muotoon. Tässä vaiheessa vaihtelevan taudin tasot on koodattava uudelleen 0: lla (kontrolli) ja 1: llä (tapaukset), vaikka se ei ole varsinaisesti tarpeen, ja tasoa 0 on pidettävä viiteluokkana (mikä helpottaa tulosten tulkintaa):

>kirjasto(reshape2)

>alkoholi.df < – melt(alcohol)

> nimet(alcohol.df) < – c (”tauti”,” altistuminen”,”n”)

> tasot(alkoholi.df$disease) < – c(1,0)

> alkoholi.df$disease < – relevel(alkoholi.df$disease,”0″)

logistinen regressiomalli kirjoitetaan tämän jälkeen seuraavasti:

> glm.res < – glm(disease exposure, data = alcohol.df,

perhe = binominen, painot = n)

> Yhteenveto(glm.res)

koron tulos on taulukkoon > = 80 liittyvä rivi, koska se antaa tietoa R: llä estimoituun altistusmuuttujaan liittyvän regressiokertoimen arvosta keskivirheineen sekä testin tilastollisen arvon. Tässä regressiokerroin tulkitaan kerroinsuhteen logiksi. Huomaa, että yksi olisi saada täsmälleen samat tulokset vaihtamalla rooli muuttujien edellisen formulaation, altistustauti.

edellä laskettu kerroinsuhde löytyy regressiokertoimesta, joka liittyy merkittävään tekijään (altistuminen), sekä sen 95 prosentin luottamusvälistä, seuraavasti:

> exp(coef(glm.res))

> exp(confint(glm.res))

toisessa ratkaisussa tarkastellaan tapausten määrää ja yksilöiden kokonaismäärää:

> alchol2 < – data.frame(expos = c (”<> = 80”), case = c(104,96),

total = c(104 + 666, 96 + 109))

> summary(glm(cbind(case, total-case) expos, data = alcohol2,

perhe = binominen))

4) nämä tiedot ovat peräisin tutkimuksesta, jossa pyrittiin selvittämään kehon kreatiinikinaasipitoisuuden prognostista paikkansapitävyyttä sydäninfarktin ehkäisyssä .

koska muuttujien nimeä ei ole tiedostossa, ne on annettava heti tietojen tuonnin jälkeen:

> sck < – read.pöytä (”SK.dat”, header = FALSE)

> names(sck) < – c(”ck”, ”pres”, ”abs”)

> summary(sck)

the total koehenkilöiden lukumäärä vastaa kahden muuttujan pres-ja ABS-frekvenssien summaa eli:

> Sum(SCK)

tai vastaavasti: summa (SCK$pres) + summa (SCK$abs) (mutta ei sck$pres + SCK$abs).

näiden kahden muuttujan suhteellisten frekvenssien laskemiseksi on tarpeen tietää kokonaissuureet muuttujittain. Nämä saadaan käyttämällä komentoa apply ja operoimalla sarakkeilla:

> ni < – apply(sck, 2, sum)

tämän perusteella riittää, että jokainen muuttujien Pres ja abs arvo jaetaan aiemmin lasketuilla arvoilla. Saadut arvot tallennetaan kahteen uuteen muuttujaan samaan datataulukkoon:

> sck$pres.prop < – sck$pres/ni

> sck$abs.prop < – SCK$abs/ni

on mahdollista tarkistaa laskelmien paikkansapitävyys: kunkin muuttujan arvojen summan on nyt oltava yhtä suuri kuin 1:

> sovelletaan(sck, 2, summa)

silloin saadut osuudet esitetään yksinkertaisesti yhdessä kaaviossa, kun otetaan huomioon muuttuja ck kuten abscissae:

> xyplot(Pres.prop + abs.prop ck, data = sck, type = c (”b”,”g”),

auto.key = TRUE, ylab = ”Frequency”)

ohjetyyppi = c(”b”, ”g”) tarkoittaa, että haluamme pisteiden näyttävän yhdistettyinä viivoilla (”b” = ”o” + ”l’) sekä ruudukolla (”g”).

ryhmitetyn aineiston regressiomalli kirjoitetaan seuraavasti:

> glm.res < – glm(cbind(pres,abs) ck, data = sck, family = binominen)

> summary(glm.res)

todennäköisyyksinä ilmaistut ennusteet, jotka eivät ole log-odds-asteikolla, saadaan komennon predict avulla määrittelemällä optio type = ”response” seuraavasti:

> glm.pred < – predict(glm.res, type = ”response”)

> names(glm.pred) < – SCK$ck

ottamalla huomioon todennäköisyyksien ≥ 0, 5 osoittavan sairastuneita, saadaan näin seuraava jakauma:

> glm.pred

on päätelty, että ihmisiä pidetään tämän mallin mukaan sairaina, jos arvot ck ovat 80 tai enemmän. Tämä on helppo todentaa kaaviolla, jossa esitetään muuttujan ck arvoihin perustuvat todennäköisyydet:

> sck$sick < – SCK$pres/(sck$pres + SCK$abs)

> xyplot(GLM.pred ~ sck$ck, type = ”l”,

ylab = ”Probability”, xlab = ”CK”,

panel = function(…) {

panel.xyplot ( … )

paneeli.xyplot (sck$CK, SCK$sick, PCH = 19, col = ”grey”)

})

ensimmäisessä vaiheessa on ”purettava” ryhmitelty tieto ja luotava taulukko, jossa on kaksi saraketta: ensimmäinen kuvaa muuttujaa ck ja toinen kuvaa taudin esiintymistä tai puuttumista. Käytämme muuttujassa ni aiemmin laskettuja ja saatavilla olevia aliryhmittymiä:

> sck.expand < – data.frame(ck = c(rep(SCK$ck, sck$pres),

rep(sck$ck, sck$abs)),

sick = c(rep(1,ni), rep(0,ni)))

>taulukko (sck.expand$sick)

> kanssa(sck.laajenna, tapply (sick, ck, sum))

kahdella viimeisellä käskyllä pyritään varmistamaan, että samat lukemat kohdataan ja että sairaiden jakauma ck: n arvon mukaan on oikea. Voidaan myös todentaa, että logistisesta regressiosta saadaan samat tulokset ja samalla voidaan lisätä edellisessä datataulukossa ennustetut arvot (samaa menettelyä kuin aiemmin voitiin noudattaa ja ennusteita voidaan toistaa, mutta tämä on yksinkertaisempaa näin):

> glm.res2 < – glm(sairas ck, data = sck.expand, family = binominen)

> sck.expand$prediction < – ifelse(predict(glm.res2,

type = ”response”) > = 0.5,

1, 0)

> kanssa(sck.expand, table(sick, prediction))

viimeinen komento mahdollistaa sekaannusmatriisin näyttämisen, jossa reaalidiagnostiikka ja regressiomallin ennustamat risteytetään. Oikeita luokitusasteita voidaan verrata, kun vaihtelemme viitekynnystä:

> classif.tab < – kanssa(sck.expand, table(sick, prediction))

> sum(diag (classif.tab))/summa(classif.tab)

ROC – käyrän näyttämiseksi käytetään pakettia ROCR:

> library(ROCR)

> pred < – prediction(predict(glm.res2, type = ”response”),

sck.expand$sick)

> perf < – performance(pred, ”tpr”,”fpr”)

> plot(perf, ylab = ”Sensitivity”, xlab = ”1-Specificity”)

> grid()

> abline(0, 1, lty = 2)

The value of the area under the curve is obtained as follows:

> performance(pred, ”auc”)@”y.values”