6.6 applikationer

oprettelsen af datatabellen kan udføres som angivet nedenfor fra en tabel oprettet med kommandomatricen ():

> tab< – as.tabel (matrice(c (815,115,208,327), nul = 2, BYV = sand,

dimnames = liste (EST = c(“+”,”−”),

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

> tab

det skal bemærkes, at tabellen blev omorganiseret for at svare til de sædvanlige notationer med “positive” begivenheder på den første række (typisk eksponeringen) og den første kolonne (typisk sygdommen) i tabellen.

det ville da være muligt at beregne alle de krævede mængder ved hjælp af meget få aritmetiske operationer med R. Pakken epiR indeholder dog alle de værktøjer, der er nødvendige for at reagere på epidemiologiske spørgsmål i tabellerne 2 til 2. Således kommandoen epi.test () giver prævalens, følsomhed/specificitetsværdier samt positive og negative prædiktive værdier:

> bibliotek (epiR)

> epi.test (tab)

vedrørende området under kurven kan en ekstern kommando også anvendes, roc.fra.tabel (), i pakken epicalc:

> bibliotek (epicalc)

> roc.fra.tabel (tab, graf = falsk)

import af data HDI ‘ er.dat opnås som herefter:

>bp< – Læs.tabel(“HDI’ er.dat”, header = TRUE)

> str(bp)

da det kan verificeres, er der ingen etiket forbundet med niveauerne af variablerne af interesse (bpress for blodtryk, chol for kolesterolhastigheden). For at generere og knytte etiketter kan følgende kommandoer indsættes:

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

den sidste kommando konverterer de originale variabler til kvalitative variabler og associerer til deres modaliteter etiketterne defineret af blab og clab. For at kontrollere, at databasen nu er i den ønskede form, er følgende kommandoer nyttige:

> str(bp)

> oversigt (bp)

tabellen over de relative frekvenser, der er eksponeret i udsagnet, kan nu gengives:

>runde(htabs(hdis/total bpress + chol, data = bp), 2)

da vi kun vil fokusere på forholdet mellem blodtryk og myokardieinfarkt, er det nødvendigt at samle data om kolesterolniveauer. Det er med andre ord nødvendigt at akkumulere værdierne for hvert blodtryksniveau, uanset kolesterolniveauet. Muligheden bør også tages for at omdøbe niveauerne af variablen bpress ved hjælp af centrene i klasseintervallerne:

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

liste(bpress = BP), sum)

det er den sidste kommando, aggregate(), som gør det muligt at samle dataene: det opsummerer alle værdier (sum) af variablen Chol, der ikke er inkluderet i listen over variabler, der skal bevares til analyse. Samtidig gemmes resultaterne i en ny database med navnet dfrm. Bemærk, at vi i dette tilfælde ikke bruger formelnotationen, men den for den alternative syntaks for aggregat() (svarvariabel efterfulgt af en liste bestående af en eller flere forklarende variabler). En oversigt over de aggregerede data kan opnås ved at give head() som input:

> head(dfrm, 5)

Når en andel p er tilgængelig, er dens værdi på en skala, hvis enheder er logits, givet af relationsloggen(p/(1 − p)). Derfra kan følgende kommandoer udledes for at konvertere proportionerne, beregnet som HDI ‘ er / total i logit-enheder:

> logit < – funktion(h) log(H/(1-h))

> dfrm$prop < – dfrm$HDI ‘ er/dfrm$i alt

> dfrm$logit < – logit(dfrm$HDI ‘er/dfrm$i alt)

det skal bemærkes, at der er defineret en lille funktion til konvertering af værdier, som her antages at være en proportioner, i dens ækvivalente log(dfrm $ HDI’ er/dfrm $ i alt)

(1−)). Ligeledes kunne man skrive:

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

de beregnede værdier er blevet tilføjet direkte til datarammen under navnet prop og logit.

den logistiske regressionsmodel er skrevet på følgende måde:

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

formuleringen, der anvendes, cbind(hdis, total-hdis) bpress, tager højde for det faktum, at vi får adgang til grupperede data og ikke individuelle svar. Kommandoen glm () med indstillingen family = binomial svarer til en logistisk regression, som uden at gå meget i detaljer bruger som standard logit-funktionen som et kanonisk link.

regressionskoefficienterne kan vises ved at gøre brug af resume ():

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

family = binomial))

det forrige resultat inkluderer de væsentlige oplysninger til at besvare spørgsmålet om den statistiske betydning af sammenhængen mellem blodtryk og sandsynlighed for hjerteanfald: regressionskoefficienten (på log-odds skalaen) er lig med 0,024 og er signifikant ved den sædvanlige tærskel på 5% (se kolonne Pr (>|å|)).

sandsynligheden for at få et hjerteanfald i henhold til de forskellige niveauer af blodtryk under overvejelse opnås som følger:

> glm.res <-glm(cbind(HDI ‘er, total-HDI’ er) bpress, data = dfrm,

familie = binomial)

> forudsige(glm.res)

det skal bemærkes, at de mellemliggende resultater genereret af R er blevet gemt med navnet glm.res før du bruger kommandoen forudsige (). Forudsigelserne genereret af R udtrykkes i logitform, og det er muligt at sammenligne de observerede og forudsagte logits:

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

næsten alle elementer er tilgængelige for at repræsentere proportionerne af myokardieinfarkt observeret og forudsagt grafisk i henhold til blodtryksniveauet. Forudsigelserne af modellen udtrykt i form af proportioner og ikke på logodds skalaen mangler. Det kan således være ønskeligt at tegne forudsigelseskurven, det vil sige sandsynligheden for at få et hjerteanfald i henhold til blodtrykket uden at begrænse sidstnævnte til de otte værdier, der observeres for variablen bpress. Her følger en mulig løsning:

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

> f < – funktion(H)

1/(1 + eksp(-(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”,…)

})

da de (individuelle) rådata ikke er tilgængelige, er det nødvendigt at arbejde direkte med frekvenstabellen 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)

indstillingen log = FALSE sikrer, at det returnerede resultat svarer korrekt til et oddsforhold og ikke til loggen for oddsforholdet. For at opnå et asymptotisk konfidensinterval vil følgende blive anvendt:

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

tilsvarende, resume(oddsratio(alkohol)) kunne bruges til at udføre en hypotesetest på Log odds ratio (H0 : log(Kurt) = 0).

kommandoen prop.test () kan bruges til at teste hypotesen om, at andelen af personer med et dagligt forbrug er 80 g og er identisk i tilfælde såvel som i kontrollerne, hvilket angiver de værdier, der er observeret fra krydstabuleringen i problemstillingen:

> prop.test(c(96,109), c (200,775), korrekt = falsk)

denne test svarer nøjagtigt til en å-test for at teste forskellen mellem to proportioner estimeret ud fra dataene (hvis der ikke anvendes nogen kontinuitetskorrektion).

med hensyn til den logistiske regressionsmodel skal beredskabstabellen omdannes til en datatabel, hvor de to kvalitative variabler (sygdom og eksponering) er repræsenteret korrekt, det vil sige med en dataramme. Dette kan gøres ved at bruge kommandoen melt() fra pakken reshape2, som gør det muligt at omdanne dataene i tabelform til “langt” format. På dette stadium bør niveauerne af den variable sygdom omkodes med 0 (kontrol) og 1 (tilfælde), selv om dette ikke er virkelig nødvendigt, og niveauet 0 bør betragtes som referencekategori(hvilket letter fortolkningen af resultaterne):

> bibliotek (reshape2)

> alkohol.df < – smelte(alkohol)

> navne(alkohol.df) < – c(“sygdom”, “eksponering”, “n”)

> niveauer(alkohol.DF $ sygdom) < – c(1,0)

> alkohol.DF $ sygdom < – relevel (alkohol.DF $ sygdom, “0”)

den logistiske regressionsmodel skrives derefter på følgende måde:

> glm.res < – glm(sygdomseksponering, data = alkohol.DF,

familie = binomial, vægte = n)

> resume(glm.res)

resultatet af interesse er rækken forbundet med eksponering> = 80, da den giver information om værdien af regressionskoefficienten forbundet med eksponeringsvariablen som estimeret af R med dens standardfejl samt værdien af teststatistikken. Her fortolkes regressionskoefficienten som loggen for oddsforholdet. Bemærk, at man ville opnå nøjagtigt de samme resultater ved at bytte variablernes rolle i den tidligere formulering, eksponeringssygdom.

oddsforholdet beregnet ovenfor kan findes ud fra regressionskoefficienten forbundet med interessefaktoren (eksponering) såvel som dens 95% konfidensinterval på følgende måde:

> eksp(coef(glm.res))

> eksp(confint(glm.res))

den anden løsning består i at overveje antallet af tilfælde og det samlede antal individer:

> alcohol2 < – data.ramme (ekspos = c (“<> = 80”), case = c (104,96),

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

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

family = binomial))

4) disse data kommer fra en undersøgelse, der søger at fastslå den prognostiske gyldighed af kreatinkinasekoncentration i kroppen om forebyggelse af forekomsten af myokardieinfarkt .

da navnet på variablerne ikke er i datafilen, er det nødvendigt at tildele dem umiddelbart efter at data er importeret:

>sck< – Læs.tabel (“sck.dat”, header = FALSE)

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

> sammendrag(sck)

den samlede antal emner svarer til summen af frekvenserne for de to variabler pres og ABS, det vil sige:

> sum(SCK)

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

for at beregne de relative frekvenser af disse to variabler er det nødvendigt at kende de samlede mængder efter variabel. Disse kan opnås ved hjælp af kommandoen Anvend og drift af kolonner:

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

baseret på dette er det tilstrækkeligt at dividere hver værdi af variablerne pres og abs med de tidligere beregnede værdier. De opnåede værdier gemmes i to nye variabler i den samme datatabel:

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

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

det er muligt at kontrollere, at beregningerne er korrekte: summen af værdierne for hver variabel skal nu være lig med 1:

> Anvend(sck, 2, sum)

derefter er de opnåede proportioner simpelthen repræsenteret i et enkelt diagram, i betragtning af værdierne for hver variabel variablen CK som abscissae:

> hyplot(pres.prop + abs.prop ck, data = sck, type = c (“b”,”g”),

auto.key = TRUE, ylab = “frekvens”)

instruktionstypen = c (“b”,” g”) betyder, at vi ønsker, at punkterne skal vises forbundet med linjer (“b” = ” o ” + “l’) samt et gitter (“g”).

regressionsmodellen for grupperede data er skrevet som:

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

> resume (glm.res)

forudsigelserne, udtrykt i form af sandsynligheder og ikke på log-odds skalaen, opnås ved hjælp af kommandoen predict ved at angive option type = “response” som følger:

> glm.pred < – forudsige (glm.res, type = “svar”)

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

Ved at overveje, at sandsynlighederne for 0,5 udpeger personer, der er syge, opnås følgende fordeling således:

> glm.pred

det konkluderes, at folk vil blive betragtet som syge, ifølge denne model, for værdier ck på 80 eller mere. Dette verificeres let med et diagram, hvor sandsynlighederne forudsagt baseret på værdierne for variablen ck er repræsenteret:

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

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

ylab = “Sandsynlighed”, slab = “ck”,

panel = funktion(…) {

panel.( … )

panel.19,col = “grå”)

})

i et første trin er det nødvendigt at” dekomprimere ” de grupperede data og oprette en tabel med to kolonner: den første repræsenterer variablen ck og den anden repræsenterer tilstedeværelsen eller fraværet af sygdom. Vi bruger tællingerne efter undergruppe tidligere beregnet og tilgængelig i variablen ni:

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

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

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

> tabel (sck.Udvid$sick)

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

de sidste to kommandoer er beregnet til at sikre, at de samme tællinger opstår, og at fordelingen af syge mennesker i henhold til værdien af ck er korrekt. Det kan også verificeres, at de samme resultater opnås med hensyn til den logistiske regression, og på samme tid kan de værdier, der er forudsagt i den foregående datatabel, tilføjes (den samme procedure som tidligere kunne følges, og forudsigelserne kunne replikeres, men dette er enklere på denne måde):

> glm.res2 < – glm (sick ck, data = sck.Udvid, familie = binomial)

> sck.Udvid $ forudsigelse < – ifelse(forudsige(glm.res2,

type = “svar”) > = 0.5,

1, 0)

> med(sck.Udvid, tabel (syg, forudsigelse))

den sidste kommando gør det muligt at vise en forvirringsmatrice, hvor den virkelige diagnostik og dem, der forudsiges af regressionsmodellen, krydses. De korrekte klassificeringshastigheder kan sammenlignes, når vi varierer referencetærsklen:

> classif.tab < – med (sck.Udvid, tabel(syg, forudsigelse))

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

for at vise ROC-kurven bruger vi pakken ROCR:

> library(ROCR)

> pred< – prediction(predict(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”