Die Daten zu “birthweight” wurden am Baystate Medical Center, Springfield, 1986 gesammelt (Ref).
Folgende Variablen sind enthalten:
low: indicator of birth weight less than 2.5 kg (0 =
no, 1 = yes).age: mother’s age in years.lwt: mother’s weight in pounds at last menstrual
period.smoke: smoking status during pregnancy (0 = no, 1 =
yes).ht: history of hypertension (0 = no, 1 = yes).ui: presence of uterine irritability (0 = no, 1 =
yes).ftv: number of physician visits during the first
trimester.bwt: Birthweight in grams.socclass: socioeconomic status (I = highest, III =
lowest).Der Datensatz ist als “bwt_full.csv” auf Moodle abgelegt.
Importiere den Datensatz und verschaffe dir einen Überblick.
Hinterlge die Variablen low, smoke und
socclass als Faktoren.
Hier wird mit str() kontrolliert, ob der Datensatz
korrekt eingelesen wurde:
library(rio)
bwt <- import("~/Documents/Git/stats-msc-bfh/Data/bwt_full.csv")
str(bwt)## 'data.frame': 189 obs. of 9 variables:
## $ low : int 0 0 0 0 0 0 0 0 0 0 ...
## $ age : int 19 33 20 21 18 21 22 17 29 26 ...
## $ lwt : int 182 155 105 108 107 124 118 103 123 113 ...
## $ smoke : int 0 0 1 1 1 0 0 0 1 1 ...
## $ ht : int 0 0 0 0 0 0 0 0 0 0 ...
## $ ui : int 1 0 0 1 1 0 0 0 0 0 ...
## $ ftv : int 0 3 1 2 0 0 1 1 1 0 ...
## $ bwt : int 2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
## $ socclass: chr "II" "III" "I" "I" ...Hinterlegen der Variablen low, smoke und
socclass als Faktoren. Bei socclass wird
bewusst der Zusatz levelsI() gebraucht, das diese Variable
ordinal skaliert ist. R würde aber auf Grund der
alphabetischen Reihenfolge die gleichen Levels bestimmen.
bwt$low <- factor(bwt$low, levels = c(0,1), labels = c("no", "yes"))
bwt$smoke <- factor(bwt$smoke, levels = c(0,1), labels = c("no", "yes"))
bwt$socclass <- factor(bwt$socclass, levels = c("I", "II", "III"), labels = c("I", "II", "III"))
str(bwt)## 'data.frame': 189 obs. of 9 variables:
## $ low : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ age : int 19 33 20 21 18 21 22 17 29 26 ...
## $ lwt : int 182 155 105 108 107 124 118 103 123 113 ...
## $ smoke : Factor w/ 2 levels "no","yes": 1 1 2 2 2 1 1 1 2 2 ...
## $ ht : int 0 0 0 0 0 0 0 0 0 0 ...
## $ ui : int 1 0 0 1 1 0 0 0 0 0 ...
## $ ftv : int 0 3 1 2 0 0 1 1 1 0 ...
## $ bwt : int 2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...
## $ socclass: Factor w/ 3 levels "I","II","III": 2 3 1 1 1 3 1 3 1 1 ...In dieser Übung ist die abhängige Variable low
(Untergewicht des Kindes bei Geburt, das heisst < 2.5 kg) und die
unabhängige Variable smoke. Zu Beginn sollst du eine
Analyse durchführen, um einen möglichen Zusammenhang zwischen diesen
beiden Variablen zu prüfen.
smoke in
den Zeilen und der Variable low in den Spalten.oddsratio() Funktion aus dem Package epitools.
Wie ist das Resultat zu interpretieren?Erstellen der 2x2 Tabelle, mit der Variable low in den
Spalten und der Variable smoke in den Zeilen:
smoke_2by2 <- table(bwt$smoke, bwt$low, dnn = c("smoke", "low"))
smoke_2by2## low
## smoke no yes
## no 86 29
## yes 44 30Alleine die Betrachtung der Häufigkeiten lässt erahnen, dass es einen
Zusammenhang zwischen smoke und low geben
könnte.
Man könnte einen \(\chi^2\) Test
machen (chisq.test()), die Proportion der untergewichtigen
Kinder der beiden Gruppen vergleichen (prop.test()), das
Odds Ratio (oddsratio) oder das Risk Ratio
(riskratio) (bei longitudinalen Daten) analysieren. Heute
hast du gelernt, dass die Daten auch via logistische Regression
analysiert werden können.
Berechnung des OR mit oddsratio():
Zuerst sollte mit ?oddsratio überprüft werden, ob die
Tabelle richtig ausgerichtet ist. Dies ist der Fall.
library(epitools)
oddsratio(smoke_2by2)## $data
## low
## smoke no yes Total
## no 86 29 115
## yes 44 30 74
## Total 130 59 189
##
## $measure
## odds ratio with 95% C.I.
## smoke estimate lower upper
## no 1.00000 NA NA
## yes 2.01268 1.073703 3.794579
##
## $p.value
## two-sided
## smoke midp.exact fisher.exact chi.square
## no NA NA NA
## yes 0.02914865 0.0361765 0.02649064
##
## $correction
## [1] FALSE
##
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"Wie ist das Resultat zu interpretieren?
Das OR beträgt 2.01. Die Odds für ein untergewichtiges Kind bei Geburt ist bei Raucherinnen demnach um den Faktor 2.01 höher (verglichen mit nicht-Raucherinnen). Mit einem 95% CI von 1.07 - 3.79 ist dieses OR statistisch signifikant unterschiedlichen von 1 (p = 0.026). Es besteht also Evidenz dafür, dass die Kinder von Raucherinnen bei der Geburt eher Untergewichtig sind als von Nichtraucherinnen.
Es werden nun die Ergebnisse von oben via logistische Regression berechnet.
low und
smoke mit der glm() Funktion. Speichere das
Ergebnis als glm_smoke ab. Lass dir am Schluss eine
Zusammenfassung der gerechneten logistischen Regression anzeigen.Es wird das folgende Modell angepasst:
glm_smoke <- glm(low ~ smoke, data = bwt, family = binomial)
summary(glm_smoke)##
## Call:
## glm(formula = low ~ smoke, family = binomial, data = bwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.0197 -0.7623 -0.7623 1.3438 1.6599
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
## smokeyes 0.7041 0.3196 2.203 0.0276 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 229.80 on 187 degrees of freedom
## AIC: 233.8
##
## Number of Fisher Scoring iterations: 4Die Regressionskoeffizienten sind nun im Objekt
glm_smoke gespeichert. Wir können sie mit dem Zusatz
$coefficients abrufen:
glm_smoke$coefficients## (Intercept) smokeyes
## -1.0870515 0.7040592Die Regressionskoeffizienten finden sich unter Estimate.
Das Intercept \(b_0\) beträgt -1.0871
und der Steigungskoeffizient für smoke == yes (\(b_1\)) beträgt 0.7041.
\(b_0\) (Intercept) steht für die log(Odds) für ein untergewichtiges Kind der nicht exponierten Gruppe (in diesem Fall der Nichtraucherinnen). \(b_1\) entspricht dem log(Odds Ratio), für den Event “untergewichtiges Kind” der exponierten versus der nicht exponierten Gruppe. Das log(Odds Ratio) im logistischen Regressionsmodell ist also die Steigung.
Mit exp() lassen sich die Koeffizienten, bzw. deren 95%
Konfidenzintervalle exponieren:
exp(glm_smoke$coefficients)## (Intercept) smokeyes
## 0.3372093 2.0219436exp(confint(glm_smoke))## 2.5 % 97.5 %
## (Intercept) 0.2177709 0.5070199
## smokeyes 1.0818724 3.8005817Analog zur linearen Regression kann man bei der logistischen Regression die Regressionsgleichenung verwenden, um Odds (bzw. log(Odds)) zu schätzen. Hier berechnen wir die Odds für tiefes Geburtsgewicht bei Müttern, welche Rauchen.
Es gibt zwei Möglichkeiten:
log_odds_smokers <- glm_smoke$coefficients[1] + glm_smoke$coefficients[2] # analog lineares Model
exp(log_odds_smokers) # exponieren nicht vergessen!## (Intercept)
## 0.6818182oder 2) Man exponiert die Koeffizienten bereits zu Beginn. Aus \(log(Intercept) + log(Odds Ratio)\) wird \(Intercept * Odds Ratio\)
odds_smokers <- exp(glm_smoke$coefficients[1]) * exp(glm_smoke$coefficients[2])
odds_smokers## (Intercept)
## 0.6818182Ich empfehle den ersten Ansatz, weil dieser grundsätzlich einfacher ist (insbesondere bei multiplen Modellen).
Mit \(P = \frac{Odds}{1+Odds}\) können Odds in Wahrscheinlichkeiten umgerechnet werden:
log_odds_smokers <- glm_smoke$coefficients[1] + glm_smoke$coefficients[2]
p_smokers <- exp(log_odds_smokers)/(1+exp(log_odds_smokers)) # p = Odds/(1+Odds)
cat("p smokers =", p_smokers)## p smokers = 0.4054054log_odd_non_smokers <- glm_smoke$coefficients[1]
p_nonsmokers <- exp(log_odd_non_smokers)/(1+exp(log_odd_non_smokers)) # p = Odds/(1+Odds)
cat("p non-smokers =", p_nonsmokers)## p non-smokers = 0.2521739In dieser (einfachen) Situation, kann man diese Wahrscheinlichkeiten gut anhand der 2x2 Tabelle nachrechnen und nachvollziehen:
smoke_2by2## low
## smoke no yes
## no 86 29
## yes 44 30smoke_2by2[,2]/(smoke_2by2[,1] + smoke_2by2[,2]) # stimmt mit der log. Reg. überein## no yes
## 0.2521739 0.4054054Wir kommen auf die genau gleichen Werte.
glm_smoke die Variable
socclass hinzu. Speichere das neue Modell im Objekt
glm_multi_1. Lass dir die Zusammenfassung des Modells
anzeigen.socclass
den Fit des Modells? Tipp: lrtest() aus dem Paket
lmtestglm_smoke) verändert? Warum ist das
so?Das neue Modell wird wie folgt angepasst:
glm_multi_1 <- glm(low ~ smoke + socclass, data = bwt, family = binomial)
summary(glm_multi_1)##
## Call:
## glm(formula = low ~ smoke + socclass, family = binomial, data = bwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.3442 -0.8862 -0.5428 1.4964 1.9939
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8405 0.3529 -5.216 1.83e-07 ***
## smokeyes 1.1160 0.3692 3.023 0.00251 **
## socclassII 1.0841 0.4900 2.212 0.02693 *
## socclassIII 1.1086 0.4003 2.769 0.00562 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 219.97 on 185 degrees of freedom
## AIC: 227.97
##
## Number of Fisher Scoring iterations: 4Mit dem Likelihood Ratio Test werden die beiden Modelle verglichen.
Zur Erinnerung: zwei lineare Modelle vergleichen wir mit dem Befehl
anova(). Das würde hier auch gehen, man müsste aber den
P-Wert noch selber ausrechnen.
library("lmtest")
lrtest(glm_smoke, glm_multi_1)## Likelihood ratio test
##
## Model 1: low ~ smoke
## Model 2: low ~ smoke + socclass
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 2 -114.90
## 2 4 -109.99 2 9.8299 0.007336 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(glm_smoke, glm_multi_1)## Analysis of Deviance Table
##
## Model 1: low ~ smoke
## Model 2: low ~ smoke + socclass
## Resid. Df Resid. Dev Df Deviance
## 1 187 229.81
## 2 185 219.97 2 9.8299Der Fit verbessert sich statistisch signifikant (p = 0.007). Leider steht im Output der logistischen Regression kein intuitives Mass zur Quantifizierung der Güte zur Verfügung (wie beispielsweise das \(R^2\) bei linearen Modellen). Wir müssen uns deshalb (nebst inhaltlichen Überlegungen) auf den P-Wert fokussieren.
Im multiplen Modell wird der Effekt des Rauchens unter
Berücksichtigung der anderen Variablen berechnet (adjustment). In diesem
Fall korreliert smoke mit socclass: Mütter in
socclass == I rauchen öfters als die anderen). Das zeigt
auch die folgende Tabelle:
table(bwt$socclass, bwt$smoke)##
## no yes
## I 44 52
## II 16 10
## III 55 12Der Effekt des Rauchens verstärkt sich, wenn die
Variable socclass berücksichtig wird.
Was sind die Odds für ein untergewichtiges Kind einer Mutter mit
sozioökonomischen Status = III, welche nicht raucht? Um diese Frage zu
beantworten, setzten wir die entsprechenden Regressionskoeffizienten in
die Regressionsgleichung ein. Wir erhalten so die log(Odds). Um die Odds
zu erhalten, werde diese noch exponiert:
log_odds_III <- glm_multi_1$coefficients[1] + glm_multi_1$coefficients[4]
exp(log_odds_III)## (Intercept)
## 0.4809577Was sind die entsprechenden Odds, wenn eine Mutter mit sozioökonomischen Status = III raucht? Das Vorgehen erfolgt analog wie oben:
log_odds_III_smoke <- glm_multi_1$coefficients[1] + glm_multi_1$coefficients[2] + glm_multi_1$coefficients[4]
exp(log_odds_III_smoke)## (Intercept)
## 1.468186Merke: durch die log-Transformation können die Regressionskoeffizienten analog zur linearen Regression in die Regressionsgleichung eingesetzt werden, weil auf der Log-Skala die Effekte additiv sind:
\[ log(a*b) = log(a) + log(b) \]
lwt (mother’s weight in
pounds at last menstrual period) ins Modell. Speichere das Modell im
Objekt glm_multi_2 ab. Vergleiche anschliessend das neue
Modell mit glm_mulit_1 mittels Likelihood Ratio Test. Lohnt
es sich aus statistischer Sicht, die Variable lwt in das
Modell zu nehmen?glm_multi_2 wird wie folgt angepasst:
glm_multi_2 <- glm(low ~ smoke + socclass + lwt, data = bwt, family = binomial)Alternativ kann man die update() Funktion brauchen, welche etwas Schreibarbeit erspart.
glm_multi_2 <- update(glm_multi_1, . ~ . + lwt, data = bwt, family = binomial)Nun werden die Modelle mit dem Likelihood Ratio Test verglichen:
lrtest(glm_multi_1, glm_multi_2)## Likelihood ratio test
##
## Model 1: low ~ smoke + socclass
## Model 2: low ~ smoke + socclass + lwt
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 -109.99
## 2 5 -107.51 1 4.9601 0.02594 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Auch hier zeigt sich eine signifikante Verbesserung des Fits. Aus
rein statistischen Überlegungen lohnt es sich, die Varialbe
lwt zu berücksichtigen.
Wir schauen uns nun die Zusammenfassung von glm_multi_2
an:
summary(glm_multi_2)##
## Call:
## glm(formula = low ~ smoke + socclass + lwt, family = binomial,
## data = bwt)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5278 -0.9053 -0.5863 1.2878 2.0364
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.10922 0.88211 -0.124 0.90146
## smokeyes 1.06001 0.37832 2.802 0.00508 **
## socclassII 1.29009 0.51087 2.525 0.01156 *
## socclassIII 0.97052 0.41224 2.354 0.01856 *
## lwt -0.01326 0.00631 -2.101 0.03562 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 234.67 on 188 degrees of freedom
## Residual deviance: 215.01 on 184 degrees of freedom
## AIC: 225.01
##
## Number of Fisher Scoring iterations: 4Weil in der Zusammenfassung alle Koeffizienten auf der Log-Skala angezeigt werden, exponieren wir diese (inkl. 95% CIs):
cbind(OR = exp(glm_multi_2$coefficients), exp(confint(glm_multi_2))) # Das cbind() wird gebraucht, um mehrere Spalten zu kombinieren. Somit sieht der Output etwas schöner aus. Mann kann die Koeffizienten und deren 95% CIs natürlich auch in zwei Schritten ausgeben lassen.## OR 2.5 % 97.5 %
## (Intercept) 0.8965324 0.1658131 5.3953378
## smokeyes 2.8863880 1.3945161 6.1980517
## socclassII 3.6331297 1.3380529 10.0823263
## socclassIII 2.6393030 1.1927957 6.0578281
## lwt 0.9868281 0.9738926 0.9984753Interpretation:
socclass == I, haben Mütter mit socclass == II
oder socclass == III 3.6 resp. 2.6 mal höhere Odds für ein
untergewichtiges Kind. Dies, bei gleichem Rauchstatus und gleichem
Gewicht der Mutter.lwt ist unter 1. Das heisst, die
Odds für ein untergewichtiges Kind nehmen mit zunehmendem Gewicht der
Mutter ab (bei gleichem Rauchstatus und gleicher sozialen Klasse).Anmerkung: Der “kleine” Effekt von lwt ist täuschend.
Beachte, dass lwt eine kontinuierliche Variable ist und
deshalb dieser Effekt pro Pfund gilt. Ist eine Mutter 20 Pfund schwerer,
dann ist der Effekt, bzw. das OR bereits \(exp(-0.01326 \times 20) = 0.77\)!
Was sind die Odds für ein untergewichtiges Kind bei einer Mutter mit sozioökonomischem Status = II, welche raucht und 110 Pfund wiegt?
Zuerst nochmals die logarithmierten Koeffizienten:
glm_multi_2$coefficients## (Intercept) smokeyes socclassII socclassIII lwt
## -0.10922083 1.06000589 1.29009446 0.97051488 -0.01325945Nun werden die entsprechenden Werte in die Regressionsgleichung eingesetzt.
\[log(Odds) = b_{0} + b_{1}*smoke + b_{2}
* socclassII + b_{3} * socclassIII + b_{4} * lwt\]
\[=\]
\[log(Odds) = b_{0} + b_{1}*1 + b_{2} * 1
+ b_{3} * 0 + b_{4} * 110\]
Am Schluss werden die log(Odds) exponiert:
log_odds <- glm_multi_2$coefficients[1] + glm_multi_2$coefficients[2] + glm_multi_2$coefficients[3] + glm_multi_2$coefficients[5] * 110
odds <- exp(log_odds)
odds## (Intercept)
## 2.186582Diese Odds könnte man auch noch in eine Wahrscheinlichkeit umrechnen:
p <- odds/(odds+1)
p## (Intercept)
## 0.6861842Die in der Aufgabe 3 genannten Werte treffen genau auf die Frau mit der 165. Zeile zu:
bwt[165,]## low age lwt smoke ht ui ftv bwt socclass
## 165 yes 18 110 yes 0 0 0 2296 IIWir erhalten die oben berechnete Wahrscheinlichkeit auch über die fitted values:
glm_multi_2$fitted.values[165]## 165
## 0.6861842