6.6 toepassingen

het aanmaken van de gegevenstabel kan worden gedaan zoals hierna aangegeven, uit een tabel gemaakt met de opdracht matrix() :

> tab < – as.tabel (matrix(C (815,115,208,327), nrow = 2, byrow = TRUE,

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

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

> tab

opgemerkt moet worden dat de tabel werd gereorganiseerd om overeen te komen met de gebruikelijke notaties, met “positieve” gebeurtenissen op de eerste rij (typisch, de blootstelling) en de eerste kolom (typisch, de ziekte) van de tabel.

Het zou dan mogelijk zijn om alle vereiste hoeveelheden te berekenen met behulp van zeer weinig rekenkundige bewerkingen met R. Het epir-pakket bevat echter alle instrumenten die nodig zijn om te reageren op epidemiologische vragen in de 2 × 2-tabellen. Dus, het commando epi.tests() geeft prevalentie, gevoeligheid/specificiteit waarden, evenals positieve en negatieve voorspellende waarden:

> library(epiR)

> epi.tests(tab)

Voor het gebied onder de curve kan ook een extern Commando worden gebruikt, roc.van.table (), in het pakket epicalc:

> library (epicalc)

> roc.van.tabel (tab, graph = FALSE)

importeren van de data hdis.dat wordt bereikt als volgt:

> bp < – read.tabel (“hdis.dat”, header = TRUE)

> str(bp)

aangezien het kan worden geverifieerd, is er geen label geassocieerd met de niveaus van de variabelen van belang (bpress voor bloeddruk, chol voor het cholesterolpercentage). Om labels te genereren en te koppelen, kunnen de volgende commando ‘ s worden ingevoegd:

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

Het Laatste Commando zet de oorspronkelijke variabelen om in kwalitatieve variabelen en associeert met hun modaliteiten de labels gedefinieerd door blab en clab. Om te controleren of de database nu in de gewenste vorm is, zijn de volgende commando ’s nuttig:

> str (bp)

> summary (bp)

de tabel van de relatieve frequenties die in het statement worden blootgesteld, kan nu worden weergegeven:

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

aangezien we ons alleen gaan richten op de relatie tussen bloeddruk en myocardinfarct, is het noodzakelijk om gegevens over cholesterolspiegels samen te voegen. Met andere woorden, het is noodzakelijk om de waarden voor elk bloeddrukniveau te accumuleren, ongeacht de cholesterolniveaus. De mogelijkheid moet ook worden aangegrepen om de niveaus van de variabele bpress te hernoemen met behulp van de centra van de klassenintervallen:

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

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

> dfrm < – aggregate(bp,

lijst(bpress = bp), som)

Het is het laatste opdracht, de aggregate(), die het mogelijk maakt om de aggregatie van de gegevens: de bedragen van de waarden (som) van de variabele chol die niet zijn opgenomen in de lijst van variabelen worden weerhouden voor analyse. Tegelijkertijd worden de resultaten opgeslagen in een nieuwe database met de naam dfrm. Merk op dat we in dit geval geen gebruik maken van de formulenotatie maar van de alternatieve syntaxis voor aggregate() (responsvariabele gevolgd door een lijst die bestaat uit een of meer verklarende variabelen). Een overzicht van de geaggregeerde gegevens kan worden verkregen door head() als input te geven:

> head(dfrm, 5)

wanneer een verhouding p beschikbaar is, wordt de waarde op een schaal waarvan de eenheden logits zijn, gegeven door het relatielogboek(p/(1 − p)). Daaruit kunnen de volgende commando ‘ s worden afgeleid om de verhoudingen om te zetten, berekend als hdis/totaal in logit eenheden:

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

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

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

Het moet worden opgemerkt dat een kleine functie is gedefinieerd voor het omzetten van waarden van x, die hier wordt uitgegaan van een verhouding, in het equivalent log(x/(1 − x)). Ook zou men kunnen schrijven:

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

de berekende waarden zijn direct toegevoegd aan het gegevensframe onder de naam prop en logit.

het logistische regressiemodel is als volgt geschreven:

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

de gebruikte formulering, cbind(hdis, total-hdis) bpress, houdt rekening met het feit dat we toegang hebben tot gegroepeerde gegevens en niet tot individuele reacties. Het commando glm () met de optie family = binomial komt overeen met een logistische regressie, die, zonder veel in detail te gaan, standaard de logit functie gebruikt als een canonieke link.

de regressiecoëfficiënten kunnen worden weergegeven door gebruik te maken van summary ():

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

family = binomial))

het vorige resultaat bevat de essentiële informatie om de vraag te beantwoorden van de statistische significantie van het verband tussen bloeddruk en kans op een hartaanval.: de regressiecoëfficiënt (op de log-odds schaal) is gelijk aan 0,024 en is significant bij de gebruikelijke drempelwaarde van 5% (Zie kolom Pr (>|z|)).

de kans op een hartaanval volgens de verschillende niveaus van de bloeddruk in kwestie wordt als volgt verkregen:

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

family = binomial)

> predict (glm.res)

opgemerkt moet worden dat de tussentijdse resultaten die door R zijn gegenereerd, zijn opgeslagen met de naam glm.res voor het gebruik van het commando predict (). De door R gegenereerde voorspellingen worden uitgedrukt in logit-vorm en het is mogelijk om de waargenomen en voorspelde logits te vergelijken:

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

vrijwel alle elementen zijn beschikbaar om de waargenomen en voorspelde proporties van myocardinfarct grafisch weer te geven aan de hand van het niveau van de bloeddruk. De voorspellingen van het model uitgedrukt in proporties en niet op de logodds schaal ontbreken. Het kan dus wenselijk zijn om de voorspellingscurve te tekenen, dat wil zeggen de kans op een hartaanval volgens de bloeddruk, zonder deze te beperken tot de acht waargenomen waarden voor de variabele bpress. Hier volgt een mogelijke oplossing:

> 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”,…)

})

aangezien de (individuele) ruwe gegevens niet beschikbaar zijn, is het noodzakelijk om direct te werken met de frequentietabel in de probleemstelling:

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

de optie log = FALSE zorgt ervoor dat het resultaat correct overeenkomt met een odds ratio en niet met de log van de odds ratio. Om een asymptotisch betrouwbaarheidsinterval te verkrijgen wordt het volgende gebruikt:

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

evenzo kan summary(oddsratio (alcohol)) worden gebruikt om een hypothese test uit te voeren op de log odds ratio(H0 : log (θ) = 0).

het commando prop.test() kan worden gebruikt om de hypothese te testen dat het percentage personen met een dagelijkse consumptie ≥ 80 g bedraagt en zowel in de gevallen als in de controles identiek is, met vermelding van de waarden die zijn waargenomen op basis van de kruistabel in de probleemstelling:

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

Deze test is exact gelijk aan een Z-test om het verschil te testen tussen twee verhoudingen die op basis van de gegevens worden geschat (als er geen continuïteitscorrectie wordt toegepast).

wat het logistische regressiemodel betreft, moet de contingententabel worden omgezet in een gegevenstabel waarin de twee kwalitatieve variabelen (ziekte en blootstelling) goed worden weergegeven, dat wil zeggen met een gegevenskader. Dit kan worden gedaan met behulp van het commando melt() van het pakket reshape2 die het mogelijk maakt om de gegevens in tabelvorm om te zetten in “lang” formaat. In dit stadium moeten de niveaus van de variabele ziekte opnieuw worden gecodeerd met 0 (controle) en 1 (gevallen), hoewel dit niet echt nodig is en het niveau 0 moet worden beschouwd als de referentiecategorie (die de interpretatie van de resultaten vergemakkelijkt):

> library(reshape2)

> alcohol.df < – smelt (alcohol)

> namen(alcohol.df) < – c(“disease”, “exposure”, “n”)

> spiegels(alcohol.DF$disease) < – C (1,0)

> alcohol.DF$ziekte < – niveau (alcohol.DF$disease, “0”)

het logistische regressiemodel wordt dan als volgt geschreven:

> glm.res < – glm (blootstelling aan ziekten, gegevens = alcohol.df,

familie = binomiaal, gewichten = n)

> samenvatting(glm.res)

het betrokken resultaat is de rij die geassocieerd is met expositie > = 80, aangezien deze informatie verschaft over de waarde van de met de blootstellingsvariabele verbonden regressiecoëfficiënt, zoals geschat door R, met de standaardfout, en de waarde van de teststatistiek. Hier wordt de regressiecoëfficiënt geïnterpreteerd als de log van de odds ratio. Merk op dat men precies dezelfde resultaten zou verkrijgen door het wisselen van de rol van de variabelen in de vorige formulering, blootstelling ziekte.

de hierboven berekende odds ratio kan op de volgende wijze worden gevonden aan de hand van de regressiecoëfficiënt geassocieerd met de factor of interest (blootstelling) en het 95% betrouwbaarheidsinterval ervan:

> exp(coef(glm.res))

> exp (confint (glm.res))

de tweede oplossing bestaat uit het in aanmerking nemen van het aantal gevallen en het totale aantal personen:

> alcohol2 < – gegevens.frame (expos = c (“<> = 80”), case = c(104,96),

totaal = C(104 + 666, 96 + 109))

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

familie = binomiaal))

4) deze gegevens zijn afkomstig van een studie naar de prognostische validiteit van creatine kinase concentratie in het lichaam ter preventie van het optreden van myocardinfarct .

aangezien de naam van de variabelen niet in het gegevensbestand staat, is het noodzakelijk om ze onmiddellijk na het importeren toe te wijzen:

> sck < – read.tabel (“sck.dat”, header = FALSE)

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

> samenvatting(sck)

Het totale aantal onderwerpen komt overeen met de som van de frequenties voor de twee variabelen pres en abs, dat is:

> sum(sck)

of, equivalent: Som(SCK $ pres) + Som(SCK$abs) (maar niet sck$pres + SCK$abs).

om de relatieve frequenties van deze twee variabelen te berekenen, is het noodzakelijk de totale hoeveelheden per variabele te kennen. Deze kunnen worden verkregen met het commando toepassen en werken op kolommen:

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

Op basis hiervan volstaat het om elke waarde van de variabelen pres en abs te delen door de eerder berekende waarden. De verkregen waarden worden opgeslagen in twee nieuwe variabelen in dezelfde gegevenstabel:

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

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

Het is mogelijk om te controleren of de berekeningen juist zijn: de som van de waarden voor elke variabele moet nu gelijk aan 1:

> toepassen(sck, 2, som)

Dan de verhoudingen verkregen zijn gewoon vertegenwoordigd in een enkele grafiek gezien de waarden van de variabele ck als abscis:

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

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

het instructietype = c(“b”, “g”) betekent dat we willen dat de punten verbonden worden weergegeven door lijnen (“b” = “o” + “l’) en een raster (“g”).

het regressiemodel voor gegroepeerde gegevens wordt geschreven als:

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

> samenvatting (glm.res)

de voorspellingen, uitgedrukt in de vorm van waarschijnlijkheden en niet op de log-odds schaal, worden verkregen door middel van het commando predict door option type = “response” als volgt te specificeren:

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

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

Door ervan uit te gaan dat de waarschijnlijkheden ≥ 0,5 personen aanwijzen die ziek zijn, wordt de volgende verdeling verkregen:

> glm.pred

Er wordt geconcludeerd dat mensen volgens dit model als ziek worden beschouwd voor waarden ck van 80 of meer. Dit kan gemakkelijk worden geverifieerd met een grafiek waarin de voorspelde waarschijnlijkheden op basis van de waarden van de variabele ck worden weergegeven:

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

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

ylab = “Probability”, xlab = “ck”,

paneel = Functie(…) {

Paneel.xyplot ( … )

Paneel.xyplot (sck$ck, sck$sick, pch = 19, col = “grey”)

})

In een eerste stap is het noodzakelijk om de gegroepeerde gegevens te “decomprimeren” en een tabel te maken met twee kolommen: de eerste die de variabele ck voorstelt en de tweede die de aanwezigheid of afwezigheid van ziekte voorstelt. We zullen de tellingen per subgroep gebruiken die eerder zijn berekend en beschikbaar zijn in de variabele ni:

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

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

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

> table(SCK.expand$sick)

> met(sck.expand, tapply (sick, ck, sum))

de laatste twee commando ‘ s zijn bedoeld om ervoor te zorgen dat dezelfde tellingen worden aangetroffen en dat de verdeling van zieke mensen volgens de waarde van ck correct is. Er kan ook worden gecontroleerd of dezelfde resultaten worden verkregen met betrekking tot de logistische regressie en tegelijkertijd kunnen de voorspelde waarden in de vorige gegevenstabel worden toegevoegd (dezelfde procedure als voorheen zou kunnen worden gevolgd en de voorspellingen zouden kunnen worden gerepliceerd, maar dit is eenvoudiger op deze manier):

> glm.res2 < – glm (sick ck, data = sck.expand, family = binomial)

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

type = “respons”) > = 0.5,

1, 0)

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

Het Laatste commando maakt het mogelijk om een verwarmingsmatrix weer te geven waarin de werkelijke diagnostiek en die voorspeld door het regressiemodel worden gekruist. De juiste classificatiepercentages kunnen worden vergeleken wanneer we de referentiedrempel variëren:

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

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

om de ROC-curve weer te geven, zullen we gebruik maken van het pakket 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”