6.6 applikationer

skapandet av datatabellen kan göras som anges nedan, från en tabell som skapats med kommandomatrisen ():

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

dimnes = lista (EST = c(”+”,”−”),

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

> tab

det bör noteras att tabellen omorganiserades för att motsvara de vanliga notationerna, med” positiva ” händelser på första raden (vanligtvis exponeringen) och den första kolumnen (vanligtvis sjukdomen) i tabellen.

det skulle då vara möjligt att beräkna alla erforderliga kvantiteter med mycket få aritmetiska operationer med R. Paketet epiR innehåller dock alla verktyg som behövs för att svara på epidemiologiska frågor i tabellerna 2 oz 2. Således kommandot epi.tester () ger prevalens, känslighets – /specificitetsvärden samt positiva och negativa prediktiva värden:

> bibliotek(epiR)

> epi.tester (tab)

När det gäller området under kurvan kan ett externt kommando också användas, roc.från.tabell (), i paketet epicalc:

> bibliotek(epicalc)

> roc.från.tabell (tab, graph = FALSE)

importera data hdis.dat uppnås enligt följande:

> bp < – läs.tabell (”hdis.dat”, header = TRUE)

> str(bp)

som det kan verifieras finns det ingen etikett associerad med nivåerna av variablerna av intresse (bpress för blodtryck, chol för kolesterolhastigheten). För att generera och associera etiketter kan följande kommandon infogas:

> 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, etiketter = clab)

det sista kommandot konverterar de ursprungliga variablerna till kvalitativa variabler och associerar till deras modaliteter etiketterna definierade av blab och clab. För att verifiera att databasen nu är i önskad form är följande kommandon användbara:

> str(bp)

> sammanfattning (bp)

tabellen över de relativa frekvenserna som exponeras i uttalandet kan nu reproduceras:

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

eftersom vi bara ska fokusera på förhållandet mellan blodtryck och hjärtinfarkt, är det nödvändigt att aggregera data om kolesterolnivåer. Med andra ord är det nödvändigt att ackumulera värdena för varje blodtrycksnivå, oavsett kolesterolnivåerna. Möjligheten bör också tas för att byta namn på nivåerna för variabeln bpress med hjälp av centren för klassintervallen:

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

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

> dfrm < – aggregat(bp,

lista(bpress = BP), summa)

det är det sista kommandot, aggregat(), vilket gör det möjligt att aggregera data: det summerar alla värden (summa) för variabeln chol som inte ingår i listan över variabler som ska behållas för analys. Samtidigt lagras resultaten i en ny databas med namnet dfrm. Observera att vi i detta fall inte använder formelnotationen utan den alternativa syntaxen för aggregat() (svarvariabel följt av en lista som består av en eller flera förklarande variabler). En översikt över de aggregerade data kan erhållas genom att ge head () som ingång:

> head(dfrm, 5)

När en andel p är tillgänglig, dess värde på en skala vars enheter är logi ges av relationsloggen(p/(1 − p)). Därifrån kan följande kommandon härledas för att konvertera proportionerna, beräknade som hdis / total i logit-enheter:

> logit < – funktion(X) log(x/(1-x))

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

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

det bör observeras att en liten funktion har definierats för att konvertera värden x, som här antas vara en proportioner, i dess ekvivalenta logg(dfrm $ hdis/dfrm $ totalt)

det bör observeras att en liten funktion har definierats för att konvertera värden x, som här antas vara en proportioner, i dess ekvivalenta logg(x / (1 − x)). På samma sätt kan man skriva:

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

de beräknade värdena har lagts direkt till dataramen under namnet prop och logit.

den logistiska regressionsmodellen är skriven på följande sätt:

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

formuleringen som används, cbind(hdis, total-hdis) bpress, tar hänsyn till det faktum att vi har tillgång till grupperade data och inte individuella svar. Kommandot glm () med alternativet family = binomial motsvarar en logistisk regression, som utan att gå mycket i detalj använder som standard logit-funktionen som en kanonisk länk.

regressionskoefficienterna kan visas genom att använda sammanfattning ():

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

familj = binomial))

föregående resultat innehåller den väsentliga informationen för att svara på frågan om den statistiska betydelsen av sambandet mellan blodtryck och hjärtinfarkt Sannolikhet: regressionskoefficienten (på log-oddsskalan) är lika med 0,024 och är signifikant vid den vanliga tröskeln på 5% (se kolumn Pr(>|z|)).

sannolikheten för hjärtinfarkt enligt de olika nivåerna av blodtryck som behandlas erhålls enligt följande:

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

familj = binomial)

> förutsäga(glm.res)

det bör noteras att de mellanliggande resultaten som genereras av R har lagrats med namnet glm.res innan du använder kommandot förutsäga (). Förutsägelserna som genereras av R uttrycks i logit-form och det är möjligt att jämföra de observerade och förutsagda logiterna:

> cbind(dfrm, logit.predit = förutsäga (glm.res))

praktiskt taget alla element är tillgängliga för att representera proportionerna av hjärtinfarkt observerat och förutsagt grafiskt enligt blodtrycksnivån. Modellens förutsägelser uttryckta i form av proportioner och inte på logodds-skalan saknas. Det kan således vara önskvärt att rita prediktionskurvan, det vill säga sannolikheten för att få en hjärtattack enligt blodtrycket, utan att begränsa den senare till de åtta värden som observerats för variabeln bpress. Här följer en möjlig lösning:

> dfrm$prop.predit < – förutsäga(glm.res, Typ = ”svar”)

>f < – funktion(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”,…)

})

eftersom (enskilda) rådata inte är tillgängliga är det nödvändigt att arbeta direkt med frekvenstabellen i problemutlåtandet:

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

alternativet log = FALSE säkerställer att resultatet som returneras motsvarar korrekt ett oddsförhållande och inte loggen för oddsförhållandet. För att få ett asymptotiskt konfidensintervall kommer följande att användas:

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

på samma sätt kan sammanfattning(oddsratio(alkohol)) användas för att utföra ett hypotesprov på logoddsförhållandet (H0 : logg(Kubi) = 0).

kommandot prop.test () kan användas för att testa hypotesen att andelen personer med en daglig konsumtion är 80 g 80 g och är identisk i Fallen såväl som i kontrollerna, vilket indikerar de värden som observerats från korstabuleringen som ges i problemutlåtandet:

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

detta test motsvarar exakt ett Z-test för att testa skillnaden mellan två proportioner uppskattade från data (om ingen kontinuitetskorrigering används).

När det gäller den logistiska regressionsmodellen måste beredskapstabellen omvandlas till en datatabell där de två kvalitativa variablerna (sjukdom och exponering) representeras ordentligt, det vill säga med en dataram. Detta kan göras genom att använda kommandot smälta () från paketet reshape2 som gör det möjligt att omvandla data i tabellform till ”långt” format. I detta skede bör nivåerna av den variabla sjukdomen kodas om med 0 (kontroll) och 1 (fall), även om detta inte är nödvändigt och nivån 0 bör betraktas som referenskategori(vilket underlättar tolkningen av resultaten):

> bibliotek (reshape2)

> alkohol.df < – smälta(alkohol)

> namn(alkohol.df) < – c (”sjukdom”,” exponering”,”n”)

> nivåer (alkohol.DF $ sjukdom) < – c(1,0)

> alkohol.DF$sjukdom < – relevans(alkohol.DF$disease, ”0”)

den logistiska regressionsmodellen skrivs sedan på följande sätt:

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

familj = binomial, vikter = n)

> sammanfattning(glm.Res)

resultatet av intresse är raden associerad med exposition > = 80 eftersom det ger information om värdet av regressionskoefficienten associerad med exponeringsvariabeln som uppskattas av R, med dess standardfel, liksom värdet av teststatistiken. Här tolkas regressionskoefficienten som loggen för oddsförhållandet. Observera att man skulle få exakt samma resultat genom att byta variablernas roll i den tidigare formuleringen, exponeringssjukdom.

oddsförhållandet beräknat ovan kan hittas från regressionskoefficienten associerad med räntefaktorn (exponering), liksom dess 95% konfidensintervall, på följande sätt:

> exp(coef(glm.res))

> exp(confint(glm.res))

den andra lösningen består av att överväga antalet fall och det totala antalet individer:

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

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

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

familj = binomial))

4) dessa data kommer från en studie som syftar till att fastställa den prognostiska giltigheten av kreatinkinas koncentration i kroppen på förebyggande av förekomsten av hjärtinfarkt .

eftersom namnet på variablerna inte finns i datafilen är det nödvändigt att tilldela dem omedelbart efter att data importerats:

> sck < – läs.tabell (”sck.dat”, header = FALSE)

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

> sammanfattning(SCK)

total antal ämnen motsvarar summan av frekvenserna för de två variablerna pres och abs, det vill säga:

> summa(SCK)

eller likvärdigt: summa (sck $ pres) + summa(sck$abs) (men inte sck$pres + sck$abs).

för att beräkna de relativa frekvenserna för dessa två variabler är det nödvändigt att känna till de totala kvantiteterna med variabel. Dessa kan erhållas genom att använda kommandot tillämpa och arbeta med kolumner:

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

baserat på detta räcker det att dela varje värde av variablerna pres och abs med de värden som beräknats tidigare. De erhållna värdena lagras i två nya variabler i samma datatabell:

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

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

det är möjligt att verifiera att beräkningarna är korrekta: summan av värdena för varje variabel måste nu vara lika med 1:

> apply(sck, 2, sum)

då är de erhållna proportionerna helt enkelt representerade i ett enda diagram, med tanke på värdena av variabeln ck som abscissae:

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

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

instruktionstypen = c (”b”,” g”) betyder att vi vill att punkterna ska visas anslutna med linjer (”b” = ” o ” + ”l’) samt ett rutnät (”g”).

regressionsmodellen för grupperade data skrivs som:

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

> sammanfattning(glm.Res)

förutsägelserna, uttryckta i form av sannolikheter och inte på logoddsskalan, erhålls med kommandot predict genom att ange alternativet type = ”response”enligt följande:

> glm.pred < – förutsäga(glm.res, Typ = ”svar”)

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

genom att överväga att sannolikheten för att 0,5 xnumx xnumx utse individer som är sjuka erhålls således följande fördelning:

> glm.pred

det dras slutsatsen att människor kommer att betraktas som sjuka, enligt denna modell, för värden ck på 80 eller mer. Detta verifieras enkelt med ett diagram där sannolikheterna som förutses baserat på värdena för variabeln ck representeras:

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

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

ylab = ”Sannolikhet”, xlab = ”ck”,

panel = funktion(…) {

panel.xyplot ( … )

panel.xyplot (sck$ck, sck$sick, pch = 19, col = ”grå”)

})

i ett första steg är det nödvändigt att ”dekomprimera” de grupperade data och skapa en tabell med två kolumner: den första som representerar variabeln ck och den andra som representerar närvaron eller frånvaron av sjukdom. Vi kommer att använda räkningarna per undergrupp som tidigare beräknats och finns i variabeln ni:

> sck.expandera < – data.ram(ck = C (rep(sck $ ck,sck $ pres),

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

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

> tabell(sck.expandera $ sick)

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

de två sista kommandona är avsedda att säkerställa att samma räkningar påträffas och att fördelningen av sjuka personer enligt värdet på ck är korrekt. Det kan också verifieras att samma resultat erhålls angående den logistiska regressionen och samtidigt kan de värden som förutses i föregående datatabell läggas till (samma procedur som tidigare kunde följas och förutsägelserna kunde replikeras, men detta är enklare på detta sätt):

> glm.res2 < – glm(sjuk ck, data = sck.expandera, familj = binomial)

> sck.expandera$prediction < – ifelse(förutsäga(glm.res2,

typ = ”svar”) > = 0.5,

1, 0)

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

det sista kommandot gör det möjligt att visa en förvirringsmatris där den verkliga diagnostiken och de som förutses av regressionsmodellen korsas. De korrekta klassificeringshastigheterna kan jämföras när vi varierar referenströskeln:

> classif.fliken < – med (sck.expandera, tabell (sjuk, förutsägelse))

> summa(diag(classif.tab))/summa(klassificera.för att visa Roc-kurvan använder vi paketet ROCR:

> bibliotek(ROCR)

> pred < – prediction(predict(glm.res2, Typ = ”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”