I datasættet sugar2.robj ses de samme data som fra tidligere, men nu tilføjet en kolonne som angiver om personen ryger eller ej. Indlæs først data:

#indlæsning af udvidelsespakker
library(jtools)
library(ggfortify)
library(ggplot2)
#indlæs data
df_sugar <- dget("https://statepi.statnoter.dk/data/sugar2.robj")
df_sugar[sample(nrow(df_sugar), 10), ]#vis 10 tilfældige rækker fra datasættet
     BMI sugar smoke
323 25.7     0     0
197 24.8     0     0
205 29.3     0     0
478 29.6     1     0
318 28.5     0     0
367 25.9     0     0
42  26.5     0     0
328 28.0     0     0
177 27.1     0     0
597 27.1     1     1

Modellen

Først skal modellen defineres. Denne gang har vi to prediktorer, da vi ønsker at undersøge om den tidligere observerede sammenhæng mellem sukker og BMI eventuelt er confoundet af rygning. Modellen defineres næsten som tidligere:

#bemærk at der nu er angivet to prediktorer med et "+" imellem.
model_sugar <- lm(BMI ~ sugar + smoke, data = df_sugar)

Forudsætninger

Plottene til at undersøge de relevante forudsætninger laves som altid, og afslører ikke nogle kritiske afvigelser:

autoplot(object = model_sugar, which = 1:3)

Resultater

Tidligere så vi at dem der bruger sukker i kaffen i gennemsnit har et BMI der er ca. 1,2 kg/m2 lavere end dem der ikke bruger sukker i kaffen (95%-Konfidensinterval [-1,8 ; -0,6]. Spørgsmålet er altså om denne forskel faktisk skyldes brugen af sukker, eller om den i stedet skyldes en confounder - f.eks. rygning som kunne være ulige fordelt i de to grupper, og som i sig selv kan have en effekt på BMI.

Dette kan undersøges ved at ændre modellen, så rygning tilføjes som en prediktor, hvilket vi har gjort ovenfor. Resultat af beregningerne kan ses her, idet vi igen beregner et konfidensinterval for middelforskellen mellem sukker- og ikke-sukker-grupperne:

summary(model_sugar)

Call:
lm(formula = BMI ~ sugar + smoke, data = df_sugar)

Residuals:
    Min      1Q  Median      3Q     Max 
-10.187  -2.114   0.013   2.193   9.586 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  27.2870     0.1695 160.988  < 2e-16 ***
sugar1       -0.3158     0.3040  -1.039    0.299    
smoke1       -2.1573     0.3040  -7.097 3.64e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.182 on 597 degrees of freedom
Multiple R-squared:  0.1052,    Adjusted R-squared:  0.1022 
F-statistic:  35.1 on 2 and 597 DF,  p-value: 3.866e-15
confint(object = model_sugar, parm = c("sugar1", "smoke1"))
           2.5 %     97.5 %
sugar1 -0.912798  0.2812458
smoke1 -2.754356 -1.5603127

Effektplot

Man kan naturligvis også lave effektplot, idet disse laves separat for hver prediktor:

effect_plot(model = model_sugar, pred = "sugar", cat.interval.geom = "line")

effect_plot(model = model_sugar, pred = "smoke", cat.interval.geom = "line")

Fortolkning

Tidligere, hvor der alene blev kigget på effekten af “sugar”, blev det beregnet at sukker-brugere i gennemsnit har et BMI som er 1,2 kg/m2 lavere end dem som ikke bruger sukker. I dette eksempel er “smoker” tilføjet som mulig forklaring på forskellen i BMI, og resultatet viser to interessante oplysninger:

“smoke” har i sig selv en signifikant sammenhæng med BMI, idet rygere i gennemsnit har et BMI der er ca. 2,2 kg/m2 lavere end ikke-rygere (95%-konfidensinterval = [-2,8 ; -1,6].

Effekten af “sugar” er blevet mindre, eller helt forsvundet, idet den nu kun er -0,3 kg/m2 (95%-konfidensinterval = [-0,9 ; 0,3]. Dette indikerer at rygning er en confounder. Hvis rygning faktisk er en confounder, betyder det at en del af den effekt der blev konstateret tidligere faktisk skyldtes at rygning påvirker BMI, og at rygning samtidig er korreleret til brugen af sukker i kaffen, dvs. at der f.eks. er flere rygere blandt dem der benytter sukker, end blandt dem der ikke benytter sukker.

Forklaring

En yderligere forståelse for årsagen til denne confounding kan fås ved at illustrere data med boxplots. Betragt først et boxplot med BMI vist for dem der benytter sukker og dem der ikke benytter sukker:

ggplot(data = df_sugar) +
    geom_boxplot(aes(y = BMI, x = sugar))

Dette plot viser tilsyneladende at sukker-brugerne har et lavere BMI.

Men lad os nu prøve at splitte dette plot op, så man også kan se rygerne. Det er heldigvis nemt i R:

ggplot(data = df_sugar) +
    geom_boxplot(aes(y = BMI, x = sugar, color = smoke))

Kigger man nu på de 4 boxplot, ser man at BMI tydeligt afhænger af om man ryger eller ej. Kigger man kun på dem der ryger, har sukker derimod ikke længere nogen betydning, og det samme gælder for dem der ikke ryger. Confoundingen opstår fordi der er langt flere rygere blandt dem der bruger sukker - og hvis vi ignorerer rygning, som i det forrige plot, så vil disse rygere trække middel-BMI for sukkerbrugerne ned, og dermed kommer det til at se ud som om sukkerbrugerne har et lavere BMI. Dvs. confoundingen er opstået fordi rygning både er korreleret til outcome (BMI) og til prediktoren (sukker).

Hvis man gerne vil se antal rygere og sukkerbrugere så kan man nemt lave en såkaldt antalstabel som vist på side 1.5, f.eks. således

xtabs(formula =  ~ sugar + smoke, data = df_sugar)
     smoke
sugar   0   1
    0 323  77
    1  77 123

Alternativt kan man eventuelt lave et søjlediagram, som her:

ggplot(data = df_sugar) +
    geom_bar(mapping = aes(x = sugar, fill = smoke))

Eller

ggplot(data = df_sugar) +
    geom_bar(mapping = aes(x = sugar, fill = smoke), position = position_dodge())