Confounding bei binären Daten

Zum Datensatz

Der Datensatz (“hospAB.csv”) für diese Übung enthält die folgenden Variablen:

death: die abhängige, binäre Variable, ob der Patient/die Patientin gestorben ist (1) oder nicht (0)
risk: eine binäre Variable, welche einschätzt, wie gefährdet der Patient/die Patientin ist (1 = low risk, 2 = high risk)
age: eine binäre Variable für das Alter (0 = <60J., 1 = >=60J.)
sex: eine binäre Variable für das Geschlecht (1 = female, 2 = male)
hospital: eine binäre Variable, um welches Spital es sich handelt (0 = A, 1 = B)

Ziel ist es herauszufinden, ob sich die Odds um zu Sterben zwischen den beiden Spitälern unterscheiden.

Übung 1

Aufgaben

Importiere den Datensatz.
Schau dir die Struktur des Datensatzes an. Erstelle wo nötig Faktoren mit entsprechenden Labels.
Erstelle geeignete Grafiken, um einen Eindruck zur Verteilung der einzelnen Variablen zu erhalten.

Lösungen

Importiere den Datensatz.

library(rio)
data <- import("hospAB.csv")

Schau dir die Struktur des Datensatzes an. Erstelle wo nötig Faktoren mit entsprechenden Labels.

str(data)

## 'data.frame':    1600 obs. of  5 variables:
##  $ risk    : int  1 2 1 2 1 2 1 2 1 2 ...
##  $ sex     : int  2 2 2 2 1 1 1 1 2 2 ...
##  $ age     : int  0 0 1 1 0 0 1 1 0 0 ...
##  $ death   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ hospital: int  0 0 0 0 0 0 0 0 1 1 ...

data$risk <- factor(data$risk, levels = c(1,2), labels = c("low risk", "high risk"))
data$sex <- factor(data$sex, levels = c(1,2), labels = c("female", "male"))
data$age <- factor(data$age, levels = c(0,1), labels = c("<60", "60+"))
data$death <- factor(data$death, levels = c(0,1), labels = c("survived", "died"))
data$hospital <- factor(data$hospital, levels = c(0,1), labels = c("A", "B"))
str(data)

## 'data.frame':    1600 obs. of  5 variables:
##  $ risk    : Factor w/ 2 levels "low risk","high risk": 1 2 1 2 1 2 1 2 1 2 ...
##  $ sex     : Factor w/ 2 levels "female","male": 2 2 2 2 1 1 1 1 2 2 ...
##  $ age     : Factor w/ 2 levels "<60","60+": 1 1 2 2 1 1 2 2 1 1 ...
##  $ death   : Factor w/ 2 levels "survived","died": 2 2 2 2 2 2 2 2 2 2 ...
##  $ hospital: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 2 2 ...

Erstelle geeignete Grafiken, um einen Eindruck zur Verteilung der einzelnen Variablen zu erhalten.

barplot(table(data$risk, data$hospital), legend = TRUE, main = "Risk", xlab = "Hospital")

barplot(table(data$sex, data$hospital), legend = TRUE, main = "Sex", xlab = "Hospital")

barplot(table(data$age, data$hospital), legend = TRUE, main = "Age", xlab = "Hospital")

barplot(table(data$death, data$hospital), legend = TRUE, main = "Death", xlab = "Hospital")

Für die Motivierten hier noch eine etwas elegantere Lösung…

# Code ist nicht prüfungsrelevant
library(tidyverse)
datalong <- data %>% 
  pivot_longer(risk:death, names_to = "profile", values_to = "values")
ggplot(datalong, aes(x = values, fill = hospital)) + geom_bar()

Übung 2

Aufgaben

Erstelle eine 2x2 Tabelle mit dem Spital in den Zeilen und den Todesfällen in den Spalten. Was stellst du fest?
Berechne die Differenz der zwei Proportionen der gestorbenen Personen. Berechne ebenfalls das dazugehörige 95% CI und den p-Wert. Was stellst du fest?
Berechne das Odds Ratio, dessen 95% CI und den p-Wert. Was stellst du fest?

Lösungen

Erstelle eine 2x2 Tabelle mit dem Spital in den Zeilen und den Todesfällen in den Spalten. Was stellst du fest?

tab <- table(data$hospital, data$death) 
tab

##    
##     survived died
##   A      760   40
##   B      719   81

In Spital B gibt es mehr Todesfälle.

Berechne die Differenz der zwei Proportionen der gestorbenen Personen. Berechne ebenfalls das dazugehörige 95% CI und den p-Wert. Was stellst du fest?

prop.test(c(40,81), c(800, 800)) # Der erste Vektor steht für die "Events", der zweite für das n pro Gruppe

## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  c(40, 81) out of c(800, 800)
## X-squared = 14.305, df = 1, p-value = 0.0001555
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.07828847 -0.02421153
## sample estimates:
##  prop 1  prop 2 
## 0.05000 0.10125

Die Proportionsdifferenz beträgt 0.10125- 0.05000 = 0.05125 = 5.125%. Das 95% CI schneidet die H0 nicht. Der p-Wert ist mit 0.016% klein. Es scheint also Evidenz dafür zu geben, dass in Spital B wirklich mehr Menschen sterben…

Alternativ könnte auch ein Chi-Quadrat Test durchgeführt werden:

chisq.test(tab)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  tab
## X-squared = 14.305, df = 1, p-value = 0.0001555

Der Chi-squared Test kommt genau zum gleichen p-Wert, liefert jedoch keine Effektgrösse.

Berechne das Odds Ratio, dessen 95% CI und den p-Wert. Was stellst du fest?

library(epitools)
oddsratio(tab, correction = TRUE)

## $data
##        
##         survived died Total
##   A          760   40   800
##   B          719   81   800
##   Total     1479  121  1600
## 
## $measure
##    odds ratio with 95% C.I.
##     estimate    lower    upper
##   A 1.000000       NA       NA
##   B 2.135312 1.450752 3.191248
## 
## $p.value
##    two-sided
##       midp.exact fisher.exact   chi.square
##   A           NA           NA           NA
##   B 9.931952e-05 0.0001372591 0.0001554553
## 
## $correction
## [1] TRUE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

Das OR beträgt 2.14. Das heisst die Chance, um zu sterben, ist in Spital B mehr als doppelt so hoch wie in Spital A. Das 95% CI schneidet H0 deutlich nicht, der p-Wert ist wiederum gleich wie bei den anderen Analysen.

Fazit: Spital B meiden?!

Übung 3

Aufgabe

Nun geht es darum, mögliche Confounding-Variablen zu berücksichtigen. Berechne das Odds Ratio erneut, allerdings stratifiziert nach Risikogruppe. Was stellst du fest?

Tipp: Du brauchst die subset() Funktion

Lösung

Zuerst müssen Teildatensätze und je eine neue 2x2 Tabelle erstellt werden:

lowrisk <- subset(data, risk == "low risk")
highrisk <- subset(data, risk == "high risk")
tab_lowrisk <- table(lowrisk$hospital, lowrisk$death)
tab_highrisk <- table(highrisk$hospital, highrisk$death)
tab_lowrisk

##    
##     survived died
##   A      585   15
##   B      386   14

tab_highrisk

##    
##     survived died
##   A      175   25
##   B      333   67

Nun kann für jedes Stratum das Odds Ratio berechnet werden:

oddsratio(tab_highrisk)

## $data
##        
##         survived died Total
##   A          175   25   200
##   B          333   67   400
##   Total      508   92   600
## 
## $measure
##    odds ratio with 95% C.I.
##     estimate     lower    upper
##   A 1.000000        NA       NA
##   B 1.402734 0.8644944 2.340432
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   A         NA           NA         NA
##   B   0.173498    0.1876631  0.1731913
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

oddsratio(tab_lowrisk)

## $data
##        
##         survived died Total
##   A          585   15   600
##   B          386   14   400
##   Total      971   29  1000
## 
## $measure
##    odds ratio with 95% C.I.
##     estimate     lower    upper
##   A 1.000000        NA       NA
##   B 1.415058 0.6637963 2.998005
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   A         NA           NA         NA
##   B  0.3636991    0.4420656  0.3559016
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

In beiden Strata sinkt das Odds Ratio beträchtlich. Auffällig ist auch, dass beide 95% CI’s nun H0 schneiden. Die Variable “risk” scheint also ein wichtiger Confounder zu sein.

Übung 4

Aufgabe

Versuche nun, alle potentiellen Confounder zu berücksichtigen. Erstelle dazu einen Teildatensatz, in welchem nur Männer über 60 Jahre und einem hohen Risiko vertreten sind. Was findest du heraus, wenn du die Odds Ratio erneut berechnest?

Tipp: innerhalb der subset() Funktion kann man mehrere Bedingungen mit einem & aneinander knüpfen.

Lösung

Zuerst wird wieder der Teildatensatz und die 2x2 Tabelle erstellt

high_risk_profile <- subset(data, sex == "male" & risk == "high risk" & age == "60+")
tab_high_risk_profile <- table(high_risk_profile$hospital, high_risk_profile$death)
tab_high_risk_profile

##    
##     survived died
##   A       30   10
##   B      120   40

Nun das adjustierte Odds Ratio:

oddsratio(tab_high_risk_profile)

## $data
##        
##         survived died Total
##   A           30   10    40
##   B          120   40   160
##   Total      150   50   200
## 
## $measure
##    odds ratio with 95% C.I.
##      estimate     lower  upper
##   A 1.0000000        NA     NA
##   B 0.9919928 0.4547018 2.3173
## 
## $p.value
##    two-sided
##     midp.exact fisher.exact chi.square
##   A         NA           NA         NA
##   B   0.984425            1          1
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

Das Odds Ratio beträgt 1. Das heisst, wenn man für alle Confounder kontrolliert, gibt es keinen Unterschied zwischen Spital A und B in Bezug auf die Chance, zu sterben.

Übung 5

Aufgabe

Abschliessend soll die Mantel-Haenszel Methode verwendet werden. Die Funktion mhor() befindet sich im Package epiDisplay. Die Funktion braucht drei Argumente:

Outcome: die abhängige Variable
Gruppe: eine Variable welche Exponierte und nicht-Exponierte identifiziert
Confounder: eine potentielle Confounding Variable

Was ist das Resultat dieser Analyse?

Tipp: Wir müssen alle Confounding-Variablen (sex, age, risk) in eine Variable bringen (am Schluss sollten \(2*2*2 = 8\) Ausprägungen existieren). Um dies zu erreichen ist z.B. die paste() Funktion hilfreich. Hier ein Beispiel:

a <- "A"
b <- "B"
c <- "C"
abc <- paste(a,b,c, sep = " ")
abc

## [1] "A B C"

Das sep = " " definiert, wie die einzelnen Elemente getrennt werden (in diesem Fall ein Leerschlag)

Lösung

Zuerst werden alle Confounding-Variablen (sex, age, risk) in eine Variable gebracht:

data$profile <- factor(paste(data$risk,data$sex, data$age))
tab <- table(data$profile, data$death)
tab

##                       
##                        survived died
##   high risk female <60       95    5
##   high risk female 60+      119   21
##   high risk male <60        144   16
##   high risk male 60+        150   50
##   low risk female <60       297    3
##   low risk female 60+       144    6
##   low risk male <60         245    5
##   low risk male 60+         285   15

Dann können wir die mhor() Funktion anwenden:

library(epiDisplay)
mhor(data$death, data$hospital, data$profile)

## 
## Stratified analysis by  Var3 
##                            OR lower lim. upper lim. P value
## Var3 high risk female <60   1     0.0800       9.16       1
## Var3 high risk female 60+   1     0.3313       3.42       1
## Var3 high risk male <60     1     0.3086       3.55       1
## Var3 high risk male 60+     1     0.4280       2.50       1
## Var3 low risk female <60    1     0.0168      19.43       1
## Var3 low risk female 60+    1     0.0875       7.26       1
## Var3 low risk male <60      1     0.0199      10.40       1
## Var3 low risk male 60+      1     0.3011       3.84       1
## M-H combined                1     0.6455       1.55       1
## 
## M-H Chi2(1) = 0 , P value = 1 
## Homogeneity test, chi-squared 7 d.f. = 0 , P value = 1

Wir sehen, dass für jede separate Strata das Odds Ratio 1 beträgt. Wie oben kommen wir zum Schluss, dass es keinen Unterschied zwischen Spital A und B in Bezug auf die Chance zu sterben gibt, wenn man alle Confounding-Variablen berücksichtigt.

“Confounding” in einer linearen Regression

Zu den Daten

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).
smoke: smoking status during pregnancy.
bwt: Birthweight in grams.
socclass: socioeconomic status (I = highest, III = lowest).

In dieser Übung soll der kausale Effekte des Rauchens auf das Geburtsgewicht des Kindes analysiert werden. Der Datensatz ist als “bwt.csv” auf Moodle abgelegt.

Übung 1

Aufgabe

Importiere den Datensatz “bwt.csv”.
Definiere Faktoren, falls indiziert.
Erstelle Boxplots, um einen möglichen Zusammenhang zwischen smoke, socclass und der abhängigen Variable bw zu untersuchen.

Lösung

Importiere den Datensatz “bwt.csv”.
Definiere Faktoren, falls indiziert.

library(rio)
bwt <- import("bwt.csv")
bwt$low <- factor(bwt$low, levels = c(0,1), labels = c("no", "yes")) # Diese Variable wird für diese Übung nicht gebraucht.
bwt$smoke <- factor(as.numeric(bwt$smoke), levels = c(0,1), labels = c("no", "yes")) # Eine Umwandlung in einen Faktor ist hier nicht zwingend nötig.
bwt$socclass <- factor(bwt$socclass) # Die Levels müssen nicht explizit definiert werden, da nominale Variable. Die existieren Labels werden übernommen.
str(bwt)

## 'data.frame':    189 obs. of  4 variables:
##  $ low     : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ smoke   : Factor w/ 2 levels "no","yes": 1 1 2 2 2 1 1 1 2 2 ...
##  $ 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 ...

Erstelle Boxplots, um einen möglichen Zusammenhang zwischen smoke, socclass und der abhängigen Variable bw zu untersuchen.

boxplot(bwt ~ socclass, data = bwt)

Es sist eine Tendenz erkennbar, dass das Geburtsgewicht bei socclass == "I" am höchsten ist. Eine Interpretation ist aber schwierig, weil das \(n\) pro Gruppe der Grafik nicht entnommen werden kann. Ggf. ist das Geburtsgewicht in der Gruppe socclass == "II" am tiefsten.

boxplot(bwt ~ smoke, data = bwt)

Ggf. ist das Geburtsgewicht in der Gruppe smoke == "yes" etwas tiefer. Eine Interpretation ist aber auch hier schwierig, weil das \(n\) pro Gruppe der Grafik nicht entnommen werden kann.

# Für die Motivierten, Code ist nicht prüfungsrelevant ;-)
library(tidyverse)
ggplot(bwt, aes(y = bwt, x = socclass, fill = smoke)) + geom_boxplot()

Übung 2

Aufgabe

Erstelle die folgenden linearen Regressionsmodelle:
1. bw als abhängige und smoke als unabhängige Variable.
2. bw als abhängige und `socclass`` als unabhängige Variable.
3. bw als abhängige und smoke und socclass als unabhängige Variablen (= multiple lineare Regression).
Interpretiere die drei Outputs mit einem speziellen Augenmerk auf das Thema Confounding:
1. Wie unterscheiden sich die Steigungsparameter in den verschiedenen Modellen?
2. Wie beeinflussen sich die Variablen gegenseitig?
3. Welches ist die Confounding-Variable?

Lösung

Erstelle die folgenden linearen Regressionsmodelle:

lm1 <- lm(bwt ~ smoke, data = bwt)
summary(lm1)

## 
## Call:
## lm(formula = bwt ~ smoke, data = bwt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2062.9  -475.9    34.3   545.1  1934.3 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3055.70      66.93  45.653  < 2e-16 ***
## smokeyes     -283.78     106.97  -2.653  0.00867 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 717.8 on 187 degrees of freedom
## Multiple R-squared:  0.03627,    Adjusted R-squared:  0.03112 
## F-statistic: 7.038 on 1 and 187 DF,  p-value: 0.008667

Kinder von Raucherinnnen haben durchschnittlich ein 283.78 Gramm tieferes Geburtsgewicht als Frauen, welche während der Schwangerschaft nicht geraucht haben.

confint(lm1)

##                 2.5 %     97.5 %
## (Intercept) 2923.6543 3187.73697
## smokeyes    -494.7973  -72.75612

Das 95% CI beinhaltet den 0-Wert deutlich nicht, weshalb der p-Wert deutlich kleiner als 5% ist.

lm2 <- lm(bwt ~ socclass, data = bwt)
summary(lm2)

## 
## Call:
## lm(formula = bwt ~ socclass, data = bwt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2096.28  -502.72   -12.72   526.28  1887.28 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3102.72      72.92  42.548  < 2e-16 ***
## socclassII   -383.03     157.96  -2.425  0.01627 *  
## socclassIII  -297.44     113.74  -2.615  0.00965 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 714.5 on 186 degrees of freedom
## Multiple R-squared:  0.05017,    Adjusted R-squared:  0.03996 
## F-statistic: 4.913 on 2 and 186 DF,  p-value: 0.008336

Kinder von Müttern mit dem höchsten sozio-ökonoimschen Status haben durchschnittlich das höchste Geburtsgewicht. Das durchschnittliche Gewicht ist bei den anderen Gruppen 383.03, bzw. 297.44 Gramm tiefer.

confint(lm2)

##                 2.5 %     97.5 %
## (Intercept) 2958.8563 3246.58121
## socclassII  -694.6575  -71.39540
## socclassIII -521.8254  -73.04498

Auch hier schneiden die 95% CIs die 0 nicht, weshalb die p-Werte klein sind.

lm3 <- lm(bwt ~ smoke + socclass, data = bwt)
summary(lm3)

## 
## Call:
## lm(formula = bwt ~ smoke + socclass, data = bwt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2313.95  -440.22    15.78   492.14  1655.05 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3334.95      91.78  36.338  < 2e-16 ***
## smokeyes     -428.73     109.04  -3.932 0.000119 ***
## socclassII   -450.36     153.12  -2.941 0.003687 ** 
## socclassIII  -452.88     116.48  -3.888 0.000141 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 688.2 on 185 degrees of freedom
## Multiple R-squared:  0.1234, Adjusted R-squared:  0.1092 
## F-statistic: 8.683 on 3 and 185 DF,  p-value: 2.027e-05

Bevor das Modell interpretiert wird, noch etwas Modelldiagnostik (wurde oben übersprungen, sollte aber immer gemacht werden, wenn man ein Modell interpretiert):

plot(lm3$fitted.values, lm3$residuals)
abline(h = 0, col = "blue")

hist(lm3$residuals)

Die Punkte im Residuenplot sind in 6 Kolonnen aufgereiht (man muss genau hinschauen, um das zu erkennen), weil es für die Variablen smoke und socclass \(2 * 3 = 6\) mögliche Kombinationen gibt. Die Residuen sind homogen um 0 und ungefähr normal verteilt.

summary(lm3)

## 
## Call:
## lm(formula = bwt ~ smoke + socclass, data = bwt)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2313.95  -440.22    15.78   492.14  1655.05 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3334.95      91.78  36.338  < 2e-16 ***
## smokeyes     -428.73     109.04  -3.932 0.000119 ***
## socclassII   -450.36     153.12  -2.941 0.003687 ** 
## socclassIII  -452.88     116.48  -3.888 0.000141 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 688.2 on 185 degrees of freedom
## Multiple R-squared:  0.1234, Adjusted R-squared:  0.1092 
## F-statistic: 8.683 on 3 and 185 DF,  p-value: 2.027e-05

Interpretation:

Wenn smoke und socclass in das gleiche Modell genommen werden, ändern sich die Steigungsparameter im Vergelich zu den einfachen Regressionsmodellen. Bei gleichem sozio-ökonomischem Status ist das Geburtsgewicht von Kindern durschnittlich 428.73 Gramm tiefer, wenn die Mutter während der Schwangerschaft raucht. Diese Differenz ist grösser und der p-Wert kleiner, wenn man für socclass korrigiert (verlgleiche mit dem einfachen Modell oben). Die Differenz hat sich durch das Berücksichtigen des sozio-ökonomischen Status erheblich (genauer um (428.73-283.78)/283.78 = 0.510783 = 51.0783%) verändert.
Auch die Steigungsparameter bzgl. socclass sind im multiplen Modell, also unter Berücksichtigung des Rauchens, deutlich grösser und die korrespondierenden p-Werte kleiner.

Mögliche Erklärung:

Es ist denkbar, dass der sozio-ökonomische Status sowohl das Rauchverhalten sowie die restlichen Umstände während der Schwangerschaft (z.B. Ernährung, medizinische Versorgung) und somit das Geburtsgewicht des Kindes beeinflusst. Ein kausaler Zusammenhang zwischen Rauchstatus während der Schwangerschaft und dem Geburtsgewicht ist ebenfalls plausibel. Man könnte also den folgenden DAG zeichnen:

library(dagitty)
dag <- dagitty("dag{
  socclass -> birthweight
  socclass -> smoke
  smoke -> birthweight
}")

plot(graphLayout(dag))

Folglich ist der sozio-ökonomische Status ein Confounder. Mütter mit dem höchsten sozio-ökonomischen Status rauchen deutlich mehr als die anderen (die Studie stammt von 1986…). Das lässt sich einfach mit einem Chi-squared Test und einem Blick auf die standardisierten Residuen zeigen:

chisq.test(table(bwt$smoke, bwt$socclass))

## 
##  Pearson's Chi-squared test
## 
## data:  table(bwt$smoke, bwt$socclass)
## X-squared = 21.779, df = 2, p-value = 1.865e-05

chisq.test(table(bwt$smoke, bwt$socclass))$residuals

##      
##                 I          II         III
##   no  -1.88578277  0.04522852  2.22912825
##   yes  2.35084895 -0.05638265 -2.77886927

Mütter mit dem höchsten sozio-ökonomischen Status rauchen also mehr, haben gleichzeitig aber bessere Bedingungen während der Schwangerschaft. Dies fürht dazu, dass ohne Berücksichtigung der Variable socclass der Effekt des Rauchens unterschätzt wird. Gleichzeitig sind die Unterschiede in den drei sozialen Schichten deutlicher, wenn man das Rauchen berücksichtigt.