Kapitel 29 Collider Bias

In letzter Zeit gab es viele Meldungen aus der Wissenschaft, die uns aufhorchen liessen. So gab es verschiedene Studien, die angeblich zeigten, dass Raucher ein geringeres Risiko haben, sich mit SARS-CoV-2 zu infizieren als Nichtraucher. Sie finden hier einen Kommentar zu diesem Thema Are smokers protected against SARS-CoV-2 infection (COVID-19)? The origins of the myth.

Die Interpration solcher Studien ist nicht einfach und es können sich verschiedene Bias-Formen einschleichen. Siehe zum Beispiel hier Griffith et al. Collider bias undermines our understanding of COVID-19 disease risk and severity.

29.1 Pakete

library(tidyverse)
library(bindata)
library(psych)
library(Statamarkdown) # Würde nur benötigt, wenn wir die Analysen auch in Stata durchführen möchten
library(dagitty)
library(lavaan)
library(finalfit) # benötigen wir nicht in diesem Kapitel, wäre jedoch eine alternative Möglichkeit, um Modelle darzustellen
library(jtools)
library(sjPlot)

29.2 Daten simulieren

Wir simulieren Daten in denen es keinen Zusammenhang zwischen Rauchen und Covid-Infektion gibt, aber jeweils einen positiven Zusammenhang zwischen Rauchen und Hospitalisierung, sowie zwischen Covid-Infektion und Hospitalisierung. Dank an https://stackoverflow.com/questions/16089178/how-to-simulate-correlated-binary-data-with-r.

set.seed(666)
a<-sample(0:1, rep=TRUE, size=100)
b<-sample(0:1, rep=TRUE, size=100)
c<-sample(0:1, rep=TRUE, size=100)

data<-data.frame(a,b,c)
## Construct a binary correlation matrix
rho <- 0.3
rho2<-0.2
rho3<-0.0
m <- matrix(c(1,rho,rho2,rho,1,rho3, rho2, rho3, 1), ncol=3)   
m

     [,1] [,2] [,3]
[1,]  1.0  0.3  0.2
[2,]  0.3  1.0  0.0
[3,]  0.2  0.0  1.0

## Simulate 10000 x-y pairs, and check that they have the specified
## correlation structure
x <- rmvbin(1e5, margprob = c(0.07, 0.1, 0.3), bincorr = m) 
data<-data.frame(x)
#           [,1]      [,2]
# [1,] 1.0000000 0.7889613
# [2,] 0.7889613 1.0000000

colnames(data)<-c("Hospitalisation", "Smoker", "Covid")

rio::export(data, "smoker_covid_hospitalisation.dta")

Wir können nun in R das relative Risiko berechnen, als Raucher eine Covid-Infektion zu erleiden, im Vergleich zu Nichtraucher:innen. Wir sehen, dass es kein erhöhtes Risiko für Raucher:innen gibt, sich mit Covid zu infizieren (Achtung: Das sind fiktive Daten, ich habe keine Ahnung, wie es in der Realität ist).

glm.fitSmokerCovid<-glm(Covid~Smoker, data=data,  family = quasipoisson)
tab_model(glm.fitSmokerCovid, string.est="Risk Ratio")

Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

	Covid
Predictors	Risk Ratio	CI	p
(Intercept)	0.30	0.30 – 0.30	<0.001
Smoker	1.00	0.97 – 1.03	0.884
Observations	100000
R2 Nagelkerke	0.000

Wenn wir nun jedoch nur die Patient:innen einschliessen, die hospitalisiert wurden, führen wir einen Selektions-Bias ein, da wir auf einen Colider konditionieren. Rauchen erhöhrt die Wahrscheinlichkeit einer Hospitalisation und Covid erhöht die Wahrscheinlichkeit einer Hospitalisation. Deswegen ist Hospitalisation ein Collider.

testImplications <- function( covariance.matrix, sample.size ){
    library(ggm)
    tst <- function(i){ pcor.test( pcor(i,covariance.matrix), length(i)-2, sample.size )$pvalue }
tos <- function(i){ paste(i,collapse=" ") }
implications <- list(c("Covid","Smoker"))
    data.frame( implication=unlist(lapply(implications,tos)),
        pvalue=unlist( lapply( implications, tst ) ) )

}


dag<-dagitty('dag {
bb="0,0,1,1"
Covid [outcome,pos="0.248,0.497"]
Hospitalisation [adjusted,pos="0.462,0.408"]
Smoker [exposure,pos="0.265,0.156"]
Covid -> Hospitalisation
Smoker -> Hospitalisation
}')
plot(dag)

Abbildung 5.5: Directed Acyclic Graph. Der Zusammenhang zwischen Smoker und Covid ist verfälscht, falls man für die Hospitalisierung in einem Modell kontrolliert (adjusted), oder wenn man nur Patient:innen einschliesst, die hospitalisiert sind.

Wir sehen jetzt, dass die Raucher:innen eine absolute Risikoreduktion von 18% und eine relative Risikoreduktion von 25%. Wir sehen dies an der Risk Ratio, die 0.75 ist, d.h. das Risiko an Covid zu erkranken ist für die Raucher 0.75 Mal kleiner als für die Nichtraucher:innen.

use smoker_covid_hospitalisation.dta
keep if Hospitalisation==1 

cs  Covid Smoker

(93,057 observations deleted)

                 |         Smoker         |
                 |   Exposed   Unexposed  |      Total
-----------------+------------------------+-----------
           Cases |      1606        2863  |       4469
        Noncases |      1376        1098  |       2474
-----------------+------------------------+-----------
           Total |      2982        3961  |       6943
                 |                        |
            Risk |  .5385647    .7227973  |   .6436699
                 |                        |
                 |      Point estimate    |    [95% conf. interval]
                 |------------------------+------------------------
 Risk difference |        -.1842326       |   -.2069141    -.161551 
      Risk ratio |         .7451117       |    .7170315    .7742916 
 Prev. frac. ex. |         .2548883       |    .2257084    .2829685 
 Prev. frac. pop |         .1094738       |
                 +-------------------------------------------------
                               chi2(1) =   251.76  Pr>chi2 = 0.0000

data_Hospitalisiert<-data %>% 
  filter(Hospitalisation==1)

glm.fitSmokerCovidHosp<-glm(Covid~Smoker, data=data_Hospitalisiert, family=quasipoisson)
tab_model(glm.fitSmokerCovidHosp, string.est="Risk Ratio")

	Covid
Predictors	Risk Ratio	CI	p
(Intercept)	0.72	0.71 – 0.74	<0.001
Smoker	0.75	0.72 – 0.77	<0.001
Observations	6943
R2 Nagelkerke	0.030

Wir sehen, dass wir durch das Konditionieren auf die Hospitalisierung, d.h. wenn die Hospitalisierung ein Einschlusskriterium für die Studie wird, ein Bias eingeführt wird.

Das gleiche passiert, wenn wir die Hospitalisierung als eine Variable in das Regressionsmodell einfügen

Im Modell mit allen Teilnehmenden, dass heisst inklusive nicht-hospitalisierten Personen, sehen wir keinen Bias.

glm.fitSmokerCovid<-glm(Covid~Smoker, data=data, family=quasipoisson)
tab_model(glm.fitSmokerCovid, string.est="Risk Ratio")

Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

	Covid
Predictors	Risk Ratio	CI	p
(Intercept)	0.30	0.30 – 0.30	<0.001
Smoker	1.00	0.97 – 1.03	0.884
Observations	100000
R2 Nagelkerke	0.000

Der Bias erscheint, wenn wir nur die Hospitalisierten Personen in die Studie einschliessen:

glm.fitSmokerCovidHosp<-glm(Covid~Smoker, data=data_Hospitalisiert, family=quasipoisson)
tab_model(glm.fitSmokerCovidHosp, string.est="Risk Ratio")

	Covid
Predictors	Risk Ratio	CI	p
(Intercept)	0.72	0.71 – 0.74	<0.001
Smoker	0.75	0.72 – 0.77	<0.001
Observations	6943
R2 Nagelkerke	0.030

Ein Bias entsteht auch, wenn wir die Variable Hospitalisation in das Regressionsmodell einschliessen:

glm.fitSmokerCovidAdjustedHosp<-glm(Covid~Smoker+Hospitalisation, data=data, family=quasipoisson)
tab_model(glm.fitSmokerCovidAdjustedHosp, string.est="Risk Ratio")

Profiled confidence intervals may take longer time to compute.
  Use `ci_method="wald"` for faster computation of CIs.

	Covid
Predictors	Risk Ratio	CI	p
(Intercept)	0.28	0.28 – 0.28	<0.001
Smoker	0.72	0.69 – 0.74	<0.001
Hospitalisation	2.62	2.55 – 2.70	<0.001
Observations	100000
R2 Nagelkerke	0.049

Wir können dies auch noch Graphisch anschauen:

plot_summs(glm.fitSmokerCovid, glm.fitSmokerCovidHosp,glm.fitSmokerCovidAdjustedHosp, exp=TRUE,scale=TRUE, point.size=4, colors=c("#D55E00", "#56B4E9", "#009E73"))

Abbildung 22.3: Model 1: Alle Teilnehmenden, Model 2: Nur Hospitalisierte Patienten, Model 3: Hospitalisierung als Kovariable im Regressionsmodell.

Wir sehen, dass der Bias sowohl in der Studie der hospitalisierten Patient:innen besteht, als auch in der Studie aller Teilnehmenden, wenn für die Variable Hospitalisierung kontrolliert wird.

Mehr zu diesem Thema finden Sie auch hier https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7665028/pdf/41467_2020_Article_19478.pdf.