6.6 Applikasjoner

opprettelsen av datatabellen kan gjøres som angitt heretter, fra et bord opprettet med kommandoen matrix() :

> tab < – as.tabell(matrise(c(815,115,208,327), now = 2, byrow = SANT,

dimnames = liste ( EST = c(«+»,»−»),

CAD = c(«+»,»-«))))

> tab

Det skal bemerkes at tabellen ble omorganisert for å korrespondere med de vanlige notasjonene, med «positive» hendelser på første rad (vanligvis eksponeringen) og den første kolonnen (vanligvis sykdommen) i tabellen.

Det ville da være mulig å beregne alle nødvendige mengder ved å bruke svært få aritmetiske operasjoner Med R. Pakken epiR inneholder imidlertid alle verktøyene som trengs for å svare på epidemiologiske spørsmål i de 2 × 2 tabellene. Dermed kommandoen epi.tester () gir prevalens, sensitivitet / spesifisitet verdier, samt positive og negative prediktive verdier:

> bibliotek (epiR)

> epi.tester (tab)

Når det gjelder området under kurven, kan en ekstern kommando også brukes, roc.fra.tabell(), i pakken epicalc:

> bibliotek(epicalc)

> roc.fra.tabell (tab, graf = FALSE)

Importere dataene hdis.dat oppnås som heretter:

> bp < – les.tabell («hdis.dat», header = TRUE)

> str(bp)

Som det kan bekreftes, er det ingen etikett knyttet til nivåene av variablene av interesse(bpress for blodtrykk, chol for kolesterolraten). For å generere og knytte etiketter, kan følgende kommandoer settes inn:

> 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 < – faktor(bp$chol, labels = clab)

den siste kommandoen konverterer de opprinnelige variablene til kvalitative variabler og tilknyttede til deres modaliteter etikettene definert av blab og clab. For å bekrefte at databasen nå er i ønsket form, er følgende kommandoer nyttige:

> str(bp)

> sammendrag(bp)

tabellen over de relative frekvensene som er eksponert i setningen, kan nå gjengis:

> round(xtabs(hdis / total bpress + chol, data = bp), 2)

Siden vi skal fokusere på forholdet mellom blodtrykk og hjerteinfarkt, er det nødvendig å samle data om kolesterolnivå. Med andre ord er det nødvendig å samle verdiene for hvert blodtrykksnivå, uavhengig av kolesterolnivået. Muligheten bør også tas for å omdøpe nivåene av variabelen bpress ved hjelp av sentrene i klasseintervallene:

> blab2 < – c (111.5, 121.5, 131.5,141.5,151.5,161.5,176.5,191.5)

> bp$bpress< – rep(blab2, hver = 7)

> dfrm< – aggregat(bp,

list(bpress = bp), sum)

det er den siste kommandoen, aggregat(), som gjør det mulig å aggregere dataene: det summerer alle verdiene (summen) av variabelen chol som ikke er inkludert i listen over variabler som skal beholdes for analyse. Samtidig lagres resultatene i en ny database som heter dfrm. Merk at i dette tilfellet bruker vi ikke formelnotasjonen, men den alternative syntaksen for aggregat () (responsvariabel etterfulgt av en liste bestående av en eller flere forklaringsvariabler). En oversikt over de aggregerte dataene kan oppnås ved å gi hode () som inngang:

> hode(dfrm, 5)

når en andel p er tilgjengelig, er verdien på en skala hvis enheter er logits gitt av relasjonsloggen(p / (1 − p)). Derfra kan følgende kommandoer utledes for å konvertere proporsjonene, beregnet som hdis / total i logit-enheter:

> logit < – funksjon(x) logg(x/(1-x))

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

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

x/(1 − x)). På samme måte kan man skrive:

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

de beregnede verdiene er lagt direkte til datarammen under navnet prop og logit.

den logistiske regresjonsmodellen er skrevet på følgende måte:

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

formuleringen som brukes, cbind (hdis, total-hdis) bpress, tar hensyn til det faktum at vi har tilgang til grupperte data og ikke individuelle svar. Kommandoen glm () med alternativet family = binomial tilsvarer en logistisk regresjon, som, uten å gå mye i detalj, bruker som standard logit-funksjonen som en kanonisk lenke.

regresjonskoeffisientene kan vises ved å bruke sammendrag ():

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

familie = binomial))

det forrige resultatet inneholder viktig informasjon for å svare på spørsmålet om den statistiske signifikansen av sammenhengen mellom blodtrykk og hjerteinfarkt sannsynlighet: regresjonskoeffisienten (på log-odds-skalaen) er lik 0,024 og er signifikant ved den vanlige terskelen på 5% (se kolonne Pr (>|z/)).

sannsynligheten for å ha et hjerteinfarkt i henhold til de forskjellige nivåene av blodtrykk som vurderes, oppnås som følger:

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

familie = binomial)

> forutsi(glm.res)

det skal bemerkes at mellomresultatene generert Av R har blitt lagret med navnet glm.res før du bruker kommandoen predict (). Forutsigelsene generert Av R er uttrykt i logit-form, og det er mulig å sammenligne de observerte og forutsagte logittene:

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

Nesten alle elementene er tilgjengelige for å representere proporsjonene av myokardinfarkt observert og spådd grafisk i henhold til nivået av blodtrykk. Forutsigelsene til modellen uttrykt i form av proporsjoner og ikke på logodds-skalaen mangler. Det kan derfor være ønskelig å tegne prediksjonskurven, det vil si sannsynligheten for å ha et hjerteinfarkt i henhold til blodtrykket, uten å begrense sistnevnte til de åtte verdiene som er observert for variabelen bpress. Her følger en mulig løsning:

> dfrm$prop.predit< – forutsi(glm.res, type = «respons»)

> f < – funksjon(x)

1/(1 + utløpsdato (- (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»,…)

})

a)

hva er verdien av oddsraten og 95% konfidensintervall (woolf-metoden)? Er det et godt estimat av den relative risikoen?

b)

er andelen forbrukere i risiko den samme blant saker og blant kontroller (vurder α = 0,05)?

C)

Bygge logistisk regresjonsmodell som gjør det mulig å teste sammenhengen mellom alkoholforbruk og individers status. Er regresjonskoeffisienten signifikant?

D)

Finn verdien av den observerte odds ratio, beregnet i (b) og dens konfidensintervall basert på resultatene av regresjonsanalysen.

Siden (individuelle) rådata ikke er tilgjengelige, er det nødvendig å arbeide direkte med frekvenstabellen som er gitt i problemstillingen:

> 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:

>bibliotek(vcd)

>oddsratio(alkohol, log = FALSE)

alternativloggen = FALSE sikrer at resultatet som returneres tilsvarer riktig til et oddsforhold og ikke til loggen for oddsforholdet. For å oppnå et asymptotisk konfidensintervall vil følgende bli brukt:

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

på Samme måte kan sammendrag(oddsratio(alkohol)) brukes til å utføre en hypotesetest på log odds ratio (H0 : log(θ) = 0).

kommandoen prop.test () kan brukes til å teste hypotesen om at andelen personer med et daglig forbruk er ≥ 80 g og er identisk i tilfellene så vel som i kontrollene, noe som indikerer verdiene observert fra krysstabuleringen gitt i problemstillingen:

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

denne testen er nøyaktig ekvivalent Med En Z-test for å teste forskjellen mellom to proporsjoner estimert fra dataene (hvis ingen kontinuitetskorreksjon blir ansatt).

med hensyn til logistisk regresjonsmodell må beredskapstabellen omdannes til en datatabell der de to kvalitative variablene (sykdom og eksponering) er representert riktig, det vil si med en dataramme. Dette kan gjøres ved å bruke kommandoen melt () fra pakken reshape2 som gjør det mulig å transformere dataene i tabellform til» langt » format. På dette stadiet bør nivåene av den variable sykdommen omkodes med 0 (kontroll) og 1 (tilfeller), selv om dette ikke er nødvendig, og nivået 0 bør betraktes som referansekategori (som letter tolkningen av resultatene):

>bibliotek(reshape2)

>alkohol.df < – smelte(alkohol)

> navn(alkohol.df) < – c («sykdom»,» eksponering»,»n»)

> nivåer(alkohol.df$sykdom) < – c(1,0)

> alkohol.df $ sykdom < – relevel (alkohol.df$sykdom, «0»)

den logistiske regresjonsmodellen skrives da på følgende måte:

> glm.res < – glm (sykdomseksponering, data = alkohol.df,

familie = binomial, vekter = n)

> sammendrag(glm.res)

resultatet av interesse er raden forbundet med utstilling > = 80 da den gir informasjon om verdien av regresjonskoeffisienten knyttet til eksponeringsvariabelen som estimert Av R, med sin standardfeil, samt verdien av teststatistikken. Her tolkes regresjonskoeffisienten som loggen for oddsforholdet. Merk at man vil oppnå nøyaktig de samme resultatene ved å bytte rollen til variablene i forrige formulering, eksponeringssykdom.

odds ratio beregnet ovenfor kan bli funnet fra regresjonskoeffisienten forbundet med faktor av interesse(eksponering), så vel som dens 95% konfidensintervall, på følgende måte:

> exp(coef (glm.res))

> exp (confint (glm.res))

den andre løsningen består i å vurdere antall tilfeller og totalt antall individer:

> alcohol2 < – data.ramme(expos = c(«<> = 80»), sak = c (104,96),

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

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

familie = binomial))

4) disse dataene kommer fra en studie som søker å fastslå prognostisk gyldighet av kreatinkinasekonsentrasjon i kroppen på forebygging av forekomst av hjerteinfarkt .

siden navnet på variablene ikke er i datafilen, vil det være nødvendig å tilordne dem umiddelbart etter at data er importert:

> sck< – les.tabell («sck.dat», header = FALSE)

> navn(sck) < – c(«ck», «pres», «abs»)

> sammendrag(sck)

summen antall fag tilsvarer SUMMEN av frekvensene for de to variablene pres og abs, det vil si:

> sum(sck)

eller tilsvarende: sum (sck$pres) + sum ( sck$abs) (men ikke sck$pres + sck$abs).

for å beregne de relative frekvensene til disse to variablene, er det nødvendig å kjenne de totale mengdene med variabel. Disse kan oppnås ved å bruke kommandoen bruk og drift av kolonner:

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

Basert på dette er det nok å dele hver verdi av variablene pres og abs med verdiene beregnet tidligere. Verdiene som er oppnådd, lagres i to nye variabler i samme datatabell:

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

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

det er mulig å verifisere at beregningene er riktige: summen av verdiene for hver variabel må nå være lik 1:

> bruk(sck, 2, sum)

da blir de oppnådde proporsjonene bare representert i et enkelt diagram, med tanke på verdiene av variabelen ck som abscissae:

> xyplot(pres.prop + abs.prop ck, data = sck, type = c («b»,»g»),

auto.key = TRUE, ylab = «Frequency»)

instruksjonstypen = c («b»,» g») betyr at vi vil at punktene skal vises forbundet med linjer («b» = » o » + «l’) samt et rutenett («g»).

regresjonsmodellen for grupperte data er skrevet som:

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

> sammendrag(glm.res)

prognosene, uttrykt i form av sannsynligheter og ikke på log-odds-skalaen, oppnås ved hjelp av kommandoen predict ved å angi alternativet type = «response»som følger:

> glm.< – forutsi(glm.res, type = «svar»)

> navn (glm.pred) < – sck$ck

ved å vurdere at sannsynlighetene ≥ 0.5 utpeker personer som er syke, oppnås følgende fordeling:

> glm.Pred

det konkluderes med at folk vil bli ansett som syke, ifølge denne modellen, for verdier ck på 80 eller mer. Dette er lett verifisert med et diagram der sannsynlighetene spådd basert på verdiene av variabelen ck er representert:

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

> xyplot(glm.pred ~ sck$ck, type = «l»,

ylab = «Sannsynlighet», xlab = «ck»,

panel = funksjon(…) {

panel.xyplot ( … )

panel.xyplot (sck$ck, sck$sick, pch = 19, col = «grey»)

})

i et første trinn er det nødvendig å «dekomprimere» de grupperte dataene og lage et bord med to kolonner: den første representerer variabelen ck og den andre representerer tilstedeværelsen eller fraværet av sykdom. Vi vil bruke tellingene etter undergruppe tidligere beregnet og tilgjengelig i variabelen ni:

> sck.utvid < – data.ramme(ck = c(rep(sck$ck, sck$pres),

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

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

> tabell(sck.utvid $ syk)

> med (sck.utvid, tapply (sick, ck, sum))

de to siste kommandoene er ment å sikre at de samme tellingene oppstår, og at fordelingen av syke mennesker i henhold til verdien av ck er riktig. Det kan også verifiseres at de samme resultatene oppnås med hensyn til logistisk regresjon, og samtidig kan verdiene spådd i forrige datatabell legges til (samme prosedyre som tidligere kunne følges og spådommene kunne replikeres, men dette er enklere på denne måten):

> glm.res2 < – glm (syk ck, data = sck.utvid, familie = binomial)

> sck.utvid$prediksjon < – ifelse (forutsi (glm.res2,

type = «respons») > = 0.5,

1, 0)

> med(sck.utvid, tabell (sick, prediction))

den siste kommandoen gjør det mulig å vise en forvirringsmatrise der den virkelige diagnostikken og de som forutses av regresjonsmodellen, krysses. De riktige klassifiseringsratene kan sammenlignes når vi varierer referansegrensen:

> classif.tab < – med (sck.utvid, tabell (syk, prediksjon))

> sum (diag (classif.tab)) / sum (classif.for å vise ROC kurven, vil vi gjøre bruk av pakken ROCR:

> bibliotek(ROCR)

> pred < – prediksjon(forutsi (glm.res2, type = «svar»),

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”