Zum Datensatz

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.

Aufgabe

Importiere den Datensatz und verschaffe dir einen Überblick. Hinterlge die Variablen low, smoke und socclass als Faktoren.

Lösungen

Mit str() können wir die Struktur vom Datensatz anschauen:

library(rio)
bwt <- import("../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" ...

Wir hinterlegen die Variablen low, smoke und socclass als Faktoren.

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

Bei socclass würde bwt$socclass <- factor(bwt$socclass) genügen, da die Labels bereits vorliegen und die Ordnung hier richtig wäre. Mit str() kontrollieren ist immer eine gute Idee:

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

Analyse von 2x2 Tabellen

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.

Aufgaben

  1. Erstelle eine 2x2 Tabelle, mit der Variable smoke in den Zeilen und der Variable low in den Spalten.
  2. Welche Möglichkeiten kennst du bereits aus dem Unterricht, um zu untersuchen, ob zwischen den beiden Variablen ein Zusammenhang besteht?
  3. Berechne das Odds Ratio für den Event “untergewichtiges Kind” für die Raucherinnen versus die Nicht-Raucherinnen. Nutze die oddsratio() Funktion aus dem Package epitools. Wie ist das Resultat zu interpretieren?

Lösungen

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  30

Alleine 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 logistischer 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. Wir setzen hier das Argument method="wald", damit wir dieses Resultat vergleichen können mit einem glm()-Fit weiter unten.

library(epitools)
oddsratio(smoke_2by2,method="wald")$measure
##      odds ratio with 95% C.I.
## smoke estimate   lower    upper
##   no  1.000000      NA       NA
##   yes 2.021944 1.08066 3.783112

Wie ist das Resultat zu interpretieren?

Das OR beträgt 2.022. Die Odds (“Chance”) für ein untergewichtiges Kind bei Geburt ist bei Raucherinnen demnach um den Faktor 2.022 höher (verglichen mit Nicht-Raucherinnen). Mit einem 95% CI von 1.081 - 3.783 ist dieses OR statistisch signifikant unterschiedlich von 1. Es besteht also Evidenz dafür, dass die Kinder von Raucherinnen bei der Geburt eher untergewichtig sind als von Nichtraucherinnen.

Einfache logistische Regression

Aufgaben

Es werden nun die Ergebnisse von oben via logistische Regression berechnet.

  1. Berechne den Zusammenhang zwischen low (abhängige Variable) und smoke (unabhängige Variable) mit der glm() Funktion. Speichere das Ergebnis als glm_smoke ab. Lass dir am Schluss eine Zusammenfassung der gerechneten logistischen Regression anzeigen.
  2. Welches sind die Regressionskoeffizienten?
  3. Was bedeuten die Regressionskoeffizienten?
  4. Lass dir die beiden Regressionskoeffizienten in exponierter Form anzeigen.
  5. Lass dir die 95% Konfidenzintervalle der beiden Regressionskoeffizienten in exponierter Form anzeigen.
  6. Berechne die Odds für tiefes Geburtsgewicht bei Müttern, welche rauchen.
  7. Berechne die Wahrscheinlichkeit für ein tiefes Geburtsgewicht für die Raucherinnen und die Nichtraucherinnen.
  8. Überprüfe anhand der 2x2 Tabelle die in 7. berechneten Wahrscheinlichkeiten.

Lösungen

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)
## 
## 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: 4

Die 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.7040592

Die Regressionskoeffizienten finden sich unter Estimate. Das geschätzte Intercept \(\hat{\beta}_0\) beträgt -1.0871 und der geschätzte Koeffizient für smoke == yes (\(\hat{\beta}_1\)) beträgt 0.7041.

\(\hat{\beta}_0\) (Intercept) steht für die logits für ein untergewichtiges Kind der nicht exponierten Gruppe (in diesem Fall der Nichtraucherinnen). \(\hat{\beta}_1\) entspricht dem log(Odds Ratio) für das Ereignis “untergewichtiges Kind” der exponierten versus der nicht exponierten Gruppe.

Mit exp() lassen sich die Koeffizienten exponieren:

exp(glm_smoke$coefficients)
## (Intercept)    smokeyes 
##   0.3372093   2.0219436

(Profile)-Likelihood-Konfidenzintervalle für das wahre \(OR\) bekommen wir mit confint(). Das ist die Default-Methode für glm-Objekte. Will man explizit Wald-Intervalle (approximative Normalverteilung), kann man confint.default() brauchen.

exp(confint(glm_smoke)) ##Likelhood-Konfidenzintervalle (default-Einstellung für glm())
##                 2.5 %    97.5 %
## (Intercept) 0.2177709 0.5070199
## smokeyes    1.0818724 3.8005817
exp(confint.default(glm_smoke)) ## Wald-Konfidenzintervalle 
##                 2.5 %    97.5 %
## (Intercept) 0.2213695 0.5136665
## smokeyes    1.0806599 3.7831106

Die Wald-Konfidenzintervalle sind äquivalent zum Resultat mit oddsratio() weiter oben.

Es ist immer gut, mit verschiedenen Methoden/Funktionen das Ergebnis zu reproduzieren. Wir bekommen diese Resultate auch mit parameters() aus dem gleichnamigen Packet. Mit exponentiate=TRUE haben wir Schätzwerte und 95% Konfidenzintervalle für das wahre \(OR\), hier wieder mit entweder Likelihood (method="profile") oder Wald-Konfidenzintervallen (method="wald").

parameters::parameters(glm_smoke,ci_method="profile",exponentiate=TRUE,digits=3)
## Parameter   | Odds Ratio |    SE |         95% CI |      z |      p
## -------------------------------------------------------------------
## (Intercept) |      0.337 | 0.072 | [0.218, 0.507] | -5.062 | < .001
## smoke [yes] |      2.022 | 0.646 | [1.082, 3.801] |  2.203 | 0.028
parameters::parameters(glm_smoke,ci_method="wald",exponentiate=TRUE,digits=3)
## Parameter   | Odds Ratio |    SE |         95% CI |      z |      p
## -------------------------------------------------------------------
## (Intercept) |      0.337 | 0.072 | [0.221, 0.514] | -5.062 | < .001
## smoke [yes] |      2.022 | 0.646 | [1.081, 3.783] |  2.203 | 0.028

Analog zur linearen Regression kann man bei der logistischen Regression den geschätzten linearen Prädiktor verwenden, um Odds (bzw. logits) zu schätzen. Hier berechnen wir die Odds für tiefes Geburtsgewicht bei Müttern, welche rauchen. Der lineare Prädiktor ist:

logit_smokers <-  glm_smoke$coefficients[1] + glm_smoke$coefficients[2] # analog lineares Model

Der lineare Prädiktor ist auf \(logit\)-Skala. Mit Exponieren erhalten wir die Odds:

unname(exp(logit_smokers))
## [1] 0.6818182

Die Wahrscheinlichkeit bekommen wir mit der logistischen Funktion:

logit_smokers <- glm_smoke$coefficients[1] + glm_smoke$coefficients[2]
p_smokers <- exp(logit_smokers)/(1+exp(logit_smokers)) # p = Odds/(1+Odds)
unname(p_smokers)
## [1] 0.4054054
logit_non_smokers <- glm_smoke$coefficients[1]
p_nonsmokers <- exp(logit_non_smokers)/(1+exp(logit_non_smokers)) # p = Odds/(1+Odds)
unname(p_nonsmokers)
## [1] 0.2521739

In dieser (einfachen) Situation mit einem zweiwertigen Prädiktor kann man diese Wahrscheinlichkeiten gut anhand der 2x2 Tabelle nachrechnen und nachvollziehen:

smoke_2by2
##      low
## smoke no yes
##   no  86  29
##   yes 44  30
smoke_2by2[,2]/(smoke_2by2[,1] + smoke_2by2[,2]) # stimmt mit der log. Reg. überein
##        no       yes 
## 0.2521739 0.4054054

Diese bedingten Wahrscheinlichkeiten (gegeben Raucher oder Nicht-Raucher) bekommen wir auch mit proportions(), wenn man auf Zeilen bedingt.

proportions(smoke_2by2,margin=1)##gegeben Zeilen
##      low
## smoke        no       yes
##   no  0.7478261 0.2521739
##   yes 0.5945946 0.4054054


Multiple logistische Regression

Aufgaben

  1. Füge dem Model glm_smoke die Variable socclass hinzu. Speichere das neue Modell im Objekt glm_multi_1. Lass dir die Zusammenfassung des Modells anzeigen.
  2. Vergleiche das erste und das zweite Modell mit einem Likelihood Ratio Test. Verbessert die Hinzunahme der Variable socclass den Fit des Modells?
  3. Hat sich der Effekt des Rauchens im multiplen Modell im Vergleich zum ersten Model (glm_smoke) verändert? Warum ist das so?
  4. Was sind die Odds und die Wahrscheinlichkeit vom Ereignis für ein untergewichtiges Kind einer Mutter mit sozioökonomischem Status = III, welche nicht raucht?
  5. Was sind die entsprechenden Odds und Wahrscheinlichkeit, wenn eine Mutter mit sozioökonomischen Status = III raucht?

Lösungen

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)
## 
## 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: 4

Mit dem Likelihood Ratio Test werden die beiden Modelle verglichen. Zur Erinnerung: zwei lineare Modelle vergleichen wir mit dem Befehl anova(). (Es gibt wie bei lm-Objekten auch eine anova() Methode für glm-Objekte).

anova(glm_smoke, glm_multi_1)
## Analysis of Deviance Table
## 
## Model 1: low ~ smoke
## Model 2: low ~ smoke + socclass
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)   
## 1       187     229.81                        
## 2       185     219.97  2   9.8299 0.007336 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In diesem Test ist die Nullhypothese, dass das reduzierte Modell stimmt, dass also die Parameter von socclass Null sind. Diese Hypothese kann verworfen werden. Der Fit verbessert sich statistisch signifikant (p = 0.007). In generalisierten linearen Modellen wie der logistischen Regression werden Residual-Quadratsummen ersetzt durch die Residuen-Devianz. Der Unterschied in Residuen-Devianz zwischen dem reduzierten (\(H_0\)-Modell, “small”) und dem geschätzten vollen Modell (“large”) ist

\[ \begin{align} D_{small}-D_{large}&=-2(l_{small}-l_{saturated})-(-2(l_{large}-l_{saturated}))\\ &=-2(l_{small}-l_{large}) \end{align} \]

Diese Grösse ist \(\chi^2\)-verteilt ist mit \(p_{large}-p_{small}=2\) Freiheitsgraden. Wie reproduzieren das Ergebnis von oben mit

## LRT statistik: Unterschied der Devianz
LR<-deviance(glm_smoke)-deviance(glm_multi_1)
LR
## [1] 9.829889
## Verteilung der Statistik
curve(dchisq(x,df=2),from=0,to=15)
rug(LR,add=TRUE)

## p-Wert
1-pchisq(LR,2)
## [1] 0.007336125

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  12

Der Effekt des Rauchens verstärkt sich, wenn die Variable socclass berücksichtig wird.

Was sind die Odds und die Wahrscheinlichkeit 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 den linearen Prädiktor ein. Wir erhalten so die logits. Um die Odds zu erhalten, werde diese noch exponiert:

logit_III_NS <- glm_multi_1$coefficients[1] + glm_multi_1$coefficients[4]
unname(logit_III_NS) ## logit
## [1] -0.731976
unname(exp(logit_III_NS))## odds 
## [1] 0.4809577
unname(exp(logit_III_NS)/(1+exp(logit_III_NS))) ## probability
## [1] 0.3247613

Mit predict() können wir Vorhersagen machen auf logit-Skala und auf Response-Skala (Wahrscheinlichkeiten):

predict(glm_multi_1,newdata=data.frame(smoke="no",socclass="III"))#logit
##         1 
## -0.731976
predict(glm_multi_1,newdata=data.frame(smoke="no",socclass="III"),type="response")#prob
##         1 
## 0.3247613

Was sind die entsprechende Odds und Wahrscheinlichkeit, wenn eine Mutter mit sozioökonomischen Status = III raucht? Das Vorgehen erfolgt analog wie oben:

logit_III_S <- glm_multi_1$coefficients[1] + glm_multi_1$coefficients[2] + glm_multi_1$coefficients[4]
unname(logit_III_S) ## logit
## [1] 0.3840279
unname(exp(logit_III_S)) ## odds
## [1] 1.468186
unname(exp(logit_III_S)/(1+exp(logit_III_S))) ## probability
## [1] 0.5948442
predict(glm_multi_1,newdata=data.frame(smoke="yes",socclass="III"))#logit
##         1 
## 0.3840279
predict(glm_multi_1,newdata=data.frame(smoke="yes",socclass="III"),type="response")#prob
##         1 
## 0.5948442

Logistische Regression mit quantitativen Prädiktorvariablen

Aufgaben

  1. Nimm zusätzlich die Variable 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?
  2. Schaue dir den Output als Zusammenfassung an. Lass dir alle Regressionskoeffizienten sowie deren 95% CI’s in exponierter Form anzeigen, damit du diese einfacher interpretieren kannst. Was sagen dir die vier Koeffizienten?
  3. Was sind die Odds und die Wahrscheinlichkeit für ein untergewichtiges Kind bei einer Mutter mit sozioökonomischem Status = II, welche raucht und 110 Pfund wiegt?

Lösungen

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:

anova(glm_multi_1, glm_multi_2)
## Analysis of Deviance Table
## 
## Model 1: low ~ smoke + socclass
## Model 2: low ~ smoke + socclass + lwt
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)  
## 1       185     219.97                       
## 2       184     215.01  1   4.9601  0.02594 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Auch 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)
## 
## 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: 4

Weil in der Zusammenfassung alle Koeffizienten auf Logit-Skala angezeigt werden, exponieren wir diese (inkl. 95% CIs). Wir können beides machen mit

cbind(OR = exp(glm_multi_2$coefficients), exp(confint(glm_multi_2)))
##                    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.9984753
## 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.

oder dann wieder einfacher mit

parameters::parameters(glm_multi_2,exponentiate=TRUE,digits=4)
## Parameter      | Odds Ratio |     SE |            95% CI |       z |     p
## --------------------------------------------------------------------------
## (Intercept)    |     0.8965 | 0.7908 | [0.1658,  5.3953] | -0.1238 | 0.901
## smoke [yes]    |     2.8864 | 1.0920 | [1.3945,  6.1981] |  2.8019 | 0.005
## socclass [II]  |     3.6331 | 1.8561 | [1.3381, 10.0823] |  2.5253 | 0.012
## socclass [III] |     2.6393 | 1.0880 | [1.1928,  6.0578] |  2.3543 | 0.019
## lwt            |     0.9868 | 0.0062 | [0.9739,  0.9985] | -2.1013 | 0.036

Interpretation:

  • Zeile 1: Achtung: das Intercept ist kein OR sondern die Odds für ein untergewichtiges Kind bei Geburt, wenn die Mutter ein Gewicht von 0 Pfund hat, den sozioökonomische Stauts = I hat und nicht raucht.
  • Zeile 2: Das Rauchen erhöht die Odds für ein untergewichtiges Kind um den Faktor 2.9, verglichen mit nicht rauchen (bei gleichem sozioökonomischem Status und gleichem Gewicht der Mutter).
  • Zeilen 3 & 4: Im Vergleich zu Müttern mit 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.
  • Zeile 5: Das OR für 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).
  • Präzision der Schätzungen: Keines der 95% CI’s beinhaltet \(H_0\) (OR = 1). Somit unterscheiden sich alle Odds Ratios signifikant von 1. Über die klinische Relevanz sagen die p-Werte nichts aus.

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.01325945

Nun werden die entsprechenden Werte in den linearen Prädiktor eingesetzt (\(I_k\) ist eine Indikatorvariable, die 1 ist, wenn die \(i\) in der entsprechenden Kategorie \(k\) ist).

\[ logit = \hat{\beta}_{0} + \hat{\beta}_{1} I_{smoke} + \hat{\beta}_{2} I_{socclassII} + \hat{\beta}_{3} I_{socclassIII} + \hat{\beta}_{4} lwt \]

\[ logit = \hat{\beta}_{0}+ \hat{\beta}_{1} + \hat{\beta}_{2}+\hat{\beta}_4 \times 110 \]

Am Schluss werden die log(Odds) exponiert:

logit <- glm_multi_2$coefficients[1] + glm_multi_2$coefficients[2] + glm_multi_2$coefficients[3] + glm_multi_2$coefficients[5] * 110
unname(logit)
## [1] 0.7823397
odds <- exp(logit)
unname(odds)
## [1] 2.186582


Diese Odds könnte man auch noch in eine Wahrscheinlichkeit umrechnen:

p <- odds/(odds+1)
unname(p)
## [1] 0.6861842


Die 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       II


Wir erhalten die oben berechnete Wahrscheinlichkeit auch über die fitted values für diese Person

glm_multi_2$fitted.values[165]
##       165 
## 0.6861842

oder allgemein mit predict()

predict(glm_multi_2,newdata=data.frame(smoke="yes",socclass="II",lwt=110))#logit
##         1 
## 0.7823397
predict(glm_multi_2,newdata=data.frame(smoke="yes",socclass="II",lwt=110),type="response")#prob
##         1 
## 0.6861842