6.6 Anwendungen

Die Erstellung der Datentabelle kann wie nachfolgend angegeben aus einer mit dem Befehl matrix() erstellten Tabelle erfolgen:

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

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

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

> tab

Es ist zu beachten, dass die Tabelle so umstrukturiert wurde, dass sie den üblichen Notationen entspricht, mit „positiven“ Ereignissen in der ersten Zeile (typischerweise die Exposition) und der ersten Spalte (typischerweise die Krankheit) der Tabelle.

Es wäre dann möglich, alle benötigten Größen mit sehr wenigen Rechenoperationen mit R zu berechnen. Das Paket epiR enthält jedoch alle Instrumente, die zur Beantwortung epidemiologischer Fragen in den 2 × 2-Tabellen erforderlich sind. Somit ist der Befehl epi.tests() liefert Prävalenz-, Sensitivitäts- / Spezifitätswerte sowie positive und negative Vorhersagewerte:

> Bibliothek(epiR)

> epi.tests(tab)

In Bezug auf die Fläche unter der Kurve kann auch ein externer Befehl verwendet werden, roc.von.tabelle(), im Paket epicalc:

> Bibliothek(epicalc)

> roc.von.tabelle(tab, graph = FALSE)

Importieren der Daten hdis.dat wird wie folgt erreicht:

> bp < – read .tabelle(„hdis.dat“, header = TRUE)

> str(bp)

Wie überprüft werden kann, ist den Pegeln der interessierenden Variablen (bpress für Blutdruck, chol für die Cholesterinrate) keine Bezeichnung zugeordnet. Um Beschriftungen zu generieren und zuzuordnen, können die folgenden Befehle eingefügt werden:

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

Der letzte Befehl konvertiert die ursprünglichen Variablen in qualitative Variablen und ordnet ihren Modalitäten die durch blab und clab definierten Labels zu. Um zu überprüfen, ob die Datenbank nun in der gewünschten Form vorliegt, sind die folgenden Befehle hilfreich:

> str(bp)

> summary(bp)

Die Tabelle der relativen Häufigkeiten, die in der Anweisung angezeigt werden, kann nun:

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

Da wir uns nur auf die Beziehung zwischen Blutdruck und Myokardinfarkt konzentrieren werden, müssen Daten zum Cholesterinspiegel aggregiert werden. Mit anderen Worten, es ist notwendig, die Werte für jedes Blutdruckniveau unabhängig vom Cholesterinspiegel zu akkumulieren. Es sollte auch die Möglichkeit genutzt werden, die Ebenen der Variablen bpress anhand der Zentren der Klassenintervalle umzubenennen:

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

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

> dfrm < – Aggregat(bp,

Liste (bpress = bp), sum)

Es ist der letzte Befehl aggregate(), der es ermöglicht, die Daten zu aggregieren: Er summiert alle Werte (sum) der Variablen chol, die nicht in der Liste der Variablen enthalten sind, die für die Analyse beibehalten werden sollen. Gleichzeitig werden die Ergebnisse in einer neuen Datenbank namens dfrm gespeichert. Beachten Sie, dass wir in diesem Fall nicht die Formelnotation verwenden, sondern die alternative Syntax für aggregate() (Antwortvariable gefolgt von einer Liste, die aus einer oder mehreren erklärenden Variablen besteht). Eine Übersicht über die aggregierten Daten erhalten Sie, indem Sie head() als Eingabe angeben:

> head(dfrm, 5)

Wenn ein Anteil p verfügbar ist, wird sein Wert auf einer Skala, deren Einheiten logits sind, durch die Beziehung log(p/(1 − p)) . Daraus können die folgenden Befehle abgeleitet werden, um die Proportionen, berechnet wie hdis / total, in Logit-Einheiten umzurechnen:

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

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

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

Es ist zu beachten, dass eine kleine Funktion definiert wurde, um Werte x, die hier als Wert angenommen werden, in ihr äquivalentes log(x/(1 − x)) umzuwandeln. Ebenso könnte man schreiben:

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

Die berechneten Werte wurden unter dem Namen prop und logit direkt zum Datenrahmen hinzugefügt.

Das logistische Regressionsmodell ist wie folgt geschrieben:

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

Die verwendete Formulierung, cbind(hdis, total-hdis) bpress, berücksichtigt die Tatsache, dass wir auf gruppierte Daten und nicht auf einzelne Antworten zugreifen. Der Befehl glm() mit der Option family = binomial entspricht einer logistischen Regression, die, ohne viel ins Detail zu gehen, standardmäßig die Funktion logit als kanonische Verknüpfung verwendet.

Die Regressionskoeffizienten können mit summary() angezeigt werden:

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

family = binomial))

Das vorherige Ergebnis enthält die wesentlichen Informationen zur Beantwortung der Frage nach der statistischen Signifikanz des Zusammenhangs zwischen Blutdruck und Herzinfarktwahrscheinlichkeit: der Regressionskoeffizient (auf der Log-Odds-Skala) ist gleich 0,024 und ist bei der üblichen Schwelle von 5% signifikant (siehe Spalte Pr(>|z|)).

Die Wahrscheinlichkeit eines Herzinfarkts entsprechend den verschiedenen betrachteten Blutdruckwerten ergibt sich wie folgt:

> glm.res < – glm(cbind(hdis, Gesamt-hdis) bpress, Daten = dfrm,

Familie = binomial)

> vorhersagen(glm.res)

Es ist zu beachten, dass die von R generierten Zwischenergebnisse unter dem Namen glm gespeichert wurden.res vor der Verwendung des Befehls predict(). Die von R erzeugten Vorhersagen werden in Logit-Form ausgedrückt und es ist möglich, die beobachteten und vorhergesagten Logits zu vergleichen:

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

Praktisch alle Elemente stehen zur Verfügung, um die beobachteten und vorhergesagten Anteile des Myokardinfarkts grafisch entsprechend dem Blutdruckniveau darzustellen. Die Vorhersagen des Modells, ausgedrückt in Form von Proportionen und nicht auf der Logodds-Skala, fehlen. Es kann daher wünschenswert sein, die Vorhersagekurve, d.h. die Wahrscheinlichkeit eines Herzinfarkts gemäß dem Blutdruck, zu zeichnen, ohne diesen auf die acht für die Variable bpress beobachteten Werte zu beschränken. Hier folgt eine mögliche Lösung:

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

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

})

Da die (einzelnen) Rohdaten nicht verfügbar sind, ist es notwendig, direkt mit der in der Problemstellung bereitgestellten Häufigkeitstabelle zu arbeiten:

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

> Bibliothek(vcd)

> oddsratio(alcohol, log = FALSE)

Die Option log = FALSE stellt sicher, dass das zurückgegebene Ergebnis korrekt einer Odds Ratio und nicht dem Log der Odds ratio entspricht. Um ein asymptotisches Konfidenzintervall zu erhalten, wird Folgendes verwendet:

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

In ähnlicher Weise könnte summary(oddsratio(alcohol)) verwendet werden, um einen Hypothesentest für das log odds ratio (H0 : log(θ) = 0) .

Der Befehl prop.test() kann verwendet werden, um die Hypothese zu testen, dass der Anteil der Personen mit einem täglichen Verbrauch ≥ 80 g beträgt und sowohl in den Fällen als auch in den Kontrollen identisch ist, wobei die Werte angegeben werden, die aus der in der Problemstellung angegebenen Kreuztabellierung beobachtet wurden:

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

Dieser Test entspricht genau einem Z-Test, um die Differenz zwischen zwei aus den Daten geschätzten Anteilen zu testen (wenn keine Kontinuitätskorrektur verwendet wird).

Im Hinblick auf das logistische Regressionsmodell muss die Kontingenztabelle in eine Datentabelle umgewandelt werden, in der die beiden qualitativen Variablen (Krankheit und Exposition) richtig, d.h. mit einem Datenrahmen, dargestellt werden. Dies kann mit dem Befehl melt() aus dem Paket reshape2 erfolgen, der es ermöglicht, die Daten in tabellarischer Form in ein „langes“ Format umzuwandeln. In diesem Stadium sollten die Ebenen der Variablen Krankheit mit 0 (Kontrolle) und 1 (Fälle) umcodiert werden, obwohl dies nicht wirklich notwendig ist und die Ebene 0 als Referenzkategorie betrachtet werden sollte (was die Interpretation der Ergebnisse erleichtert):

> library(reshape2)

> Alkohol.df < – Schmelze(Alkohol)

> Namen(Alkohol.df) < – c(„Krankheit“, „Exposition“, „n“)

> Ebenen(Alkohol.df$Krankheit) < – c(1,0)

> Alkohol.df$disease < – relevel(Alkohol.df$0, „0“)

Das logistische Regressionsmodell wird dann folgendermaßen geschrieben:

> glm.res < – glm(Krankheit Exposition, Daten = Alkohol.df,

Familie = Binomial, Gewichte = n)

> Zusammenfassung(glm.res)

Das interessierende Ergebnis ist die Zeile, die der Exposition > = 80 zugeordnet ist, da sie Informationen über den Wert des Regressionskoeffizienten liefert, der der Expositionsvariablen zugeordnet ist, wie von R geschätzt, mit seinem Standardfehler sowie dem Wert der Teststatistik. Hier wird der Regressionskoeffizient als Log der Odds Ratio interpretiert. Beachten Sie, dass man genau die gleichen Ergebnisse erzielen würde, wenn man die Rolle der Variablen in der vorherigen Formulierung Expositionskrankheit vertauscht.

Die oben berechnete Odds Ratio kann aus dem Regressionskoeffizienten, der dem Faktor von Interesse (Exposition) zugeordnet ist, sowie seinem 95% -Konfidenzintervall auf folgende Weise ermittelt werden:

> exp(coef(glm.res))

> exp(confint(glm.res))

Die zweite Lösung besteht darin, die Anzahl der Fälle und die Gesamtzahl der Personen zu berücksichtigen:

> alkohol2 < – Daten.rahmen (expos = c(„<> = 80“), Fall = c (104,96),

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

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

family = binomial))

4) Diese Daten stammen aus einer Studie, die die prognostische Validität der Kreatinkinase-Konzentration im Körper auf die Prävention des Auftretens von Myokardinfarkt untersucht .

Da der Name der Variablen nicht in der Datendatei enthalten ist, müssen sie unmittelbar nach dem Importieren der Daten zugewiesen werden:

> sck < – lesen.tabelle(„sck.dat“, header = FALSE)

> Namen(sck) < – c(„ck“, „pres“, „abs“)

> Zusammenfassung(sck)

Die Gesamtzahl der Probanden entspricht der summe der Häufigkeiten für die beiden Variablen pres und abs, also:

> sum(sck)

oder äquivalent: summe(sck $pres) + Summe(sck $abs) (aber nicht sck $pres + sck $abs).

Um die relativen Häufigkeiten dieser beiden Variablen zu berechnen, ist es notwendig, die Gesamtmengen nach Variablen zu kennen. Diese können mit dem Befehl apply und spaltenweise erhalten werden:

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

Auf dieser Grundlage genügt es, jeden Wert der Variablen pres und abs durch die zuvor berechneten Werte zu dividieren. Die erhaltenen Werte werden in zwei neuen Variablen in derselben Datentabelle gespeichert:

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

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

Es ist möglich zu überprüfen, ob die Berechnungen korrekt sind: Die Summe der Werte für jede Variable muss jetzt gleich 1 sein:

> Anwenden(sck, 2, Summe)

Dann werden die erhaltenen Anteile einfach in einem einzigen Diagramm dargestellt, wobei die Werte der Variablen ck als die Abszissen:

> xyplot(pres.stütze + abs.prop ck, Daten = sck, Typ = c(„b“, „g“),

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

Der Befehlstyp = c(„b“, „g“) bedeutet, dass die Punkte durch Linien („b“ = „o“ + „l“) sowie ein Raster („g“) verbunden angezeigt werden sollen.

Das Regressionsmodell für gruppierte Daten wird wie folgt geschrieben:

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

> Zusammenfassung(glm.res)

Die Vorhersagen, ausgedrückt in Form von Wahrscheinlichkeiten und nicht auf der Log-Odds-Skala, werden mittels des Befehls predict erhalten, indem die Option type = „response“ wie folgt angegeben wird:

> glm.pred < – vorhersagen(glm.res, type = „response“)

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

Unter Berücksichtigung, dass die Wahrscheinlichkeiten ≥ 0,5 Personen bezeichnen, die krank sind, ergibt sich somit die folgende Verteilung:

> glm.pred

Es wird der Schluss gezogen, dass Menschen nach diesem Modell für Werte ck von 80 oder mehr als krank angesehen werden. Dies lässt sich leicht mit einem Diagramm überprüfen, in dem die auf der Grundlage der Werte der Variablen ck vorhergesagten Wahrscheinlichkeiten dargestellt werden:

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

> xyplot(glm.pred ~ sck$ck, type = „l“,

ylab = „Wahrscheinlichkeit“, xlab = „ck“,

panel = Funktion(…) {

panel.xyplot(…)

Bedienfeld.xyplot(sck$ck, sck$sick, pch = 19, col = „grey“)

})

In einem ersten Schritt ist es notwendig, die gruppierten Daten zu „dekomprimieren“ und eine Tabelle mit zwei Spalten zu erstellen: die erste repräsentiert die Variable ck und die zweite repräsentiert das Vorhandensein oder Fehlen einer Krankheit. Wir verwenden die zuvor berechneten und in der Variablen ni verfügbaren Zählungen nach Untergruppen:

> sck.erweitern Sie < – Daten.rahmen(ck = c(rep(sck$ck, sck$pres),

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

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

> Tabelle(sck.expand$sick)

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

Die letzten beiden Befehle sollen sicherstellen, dass die gleichen Zählungen angetroffen werden und dass die Verteilung der Kranken nach dem Wert von ck korrekt ist. Es kann auch überprüft werden, dass die gleichen Ergebnisse in Bezug auf die logistische Regression erhalten werden und gleichzeitig die in der vorherigen Datentabelle vorhergesagten Werte hinzugefügt werden können (das gleiche Verfahren wie zuvor könnte befolgt und die Vorhersagen könnten repliziert werden, dies ist jedoch auf diese Weise einfacher):

> glm.res2 < – glm(krank ck, Daten = sck.erweitern, Familie = binomial)

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

type = „response“) > = 0.5,

1, 0)

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

Der letzte Befehl ermöglicht die Anzeige einer Verwirrungsmatrix, in der die realen Diagnosen und die vom Regressionsmodell vorhergesagten gekreuzt werden. Die korrekten Klassifizierungsraten können verglichen werden, wenn wir den Referenzschwellenwert variieren:

> classif.tab < – mit(sck.erweitern, Tabelle(krank, Vorhersage))

> Summe(diag(classif.tab))/Summe(classif.um die ROC-Kurve anzuzeigen, verwenden wir das Paket ROCR:

> library(ROCR)

> pred < – prediction(predict(glm.res2, type = „Antwort“),

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”