Exemples pour la causalité : appariement et score de propension (chapitre 2)

library(tidyverse)
library(tableone) #create 'Table 1' to describe baseline characteristics
library(Matching) #multivariate and propensity score matching with balance optimization

Dataset : right heart catheterization

Une description complète est disponible à cette adresse http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.html. Il s’agit de données sur des patients (2184 traités et 3551 contrôles) admis aux urgences dans 5 hôpitaux, la variable de traitement est swang1 (right heat catheterization vs. non) et l’outcome est death (yes or no). Nous allons considérer les variables de confusion suivantes : - cat1 : Primary disease category - age - sex`` -meanbp1`: Mean blood pressure

Matching (greedy)

Questions :

• Charger les données http://biostat.mc.vanderbilt.edu/wiki/Main/DataSets,

rhc = read_csv("./Data/rhc_data.csv")
rhc_small = rhc %>% mutate(treatment = swang1 ) %>% 
                dplyr::select(treatment,death, cat1, age,sex,meanbp1,-swang1)
                
glimpse(rhc_small)
unique(rhc_small$cat1)

• Créer la “Table 1”. Attention le package tableone ne supporte pas les variables discrètes à plusieurs facteurs, il faut donc d’abord construire la table disjonctive.

Correction

rhc_small_disjunctive = as_tibble(model.matrix(~ . , data=rhc_small)[,-1])
glimpse(rhc_small_disjunctive)
names(rhc_small_disjunctive) = c("treatment","Death","CHF","Cirr",
                                 "colcan","Coma","COPD",
                                 "lungcan","MOSF","sepsis",
                                 "age","male","meanbp1")

Correction

xvars = c("CHF","Cirr","colcan","Coma","COPD",
         "lungcan","MOSF","sepsis",
         "age","male","meanbp1") # variables de confusion
table1 = CreateTableOne(vars=xvars,strata="treatment", data=rhc_small_disjunctive, test=FALSE)
print(table1,smd=TRUE)

• Faire un matching “greedy” basé sur la distance de Mahalanobis.

Correction

greedymatch = Match(Tr=rhc_small_disjunctive$treatment,
                    M=1,
                    X=rhc_small_disjunctive[xvars],replace=FALSE)

rhc_small_disjunctive_matched = rhc_small_disjunctive %>%
                                  slice(c(greedymatch$index.control,
                                          greedymatch$index.treated)) %>%
                                  mutate( match_id = rep(c(1:2184),2))
                                

matchedtab1 = CreateTableOne(vars=xvars,strata="treatment", data=rhc_small_disjunctive_matched, test=FALSE)
print(matchedtab1, smd = TRUE)

• Sur les données “matchées” étudiez l’outcome

Correction

y_trt = rhc_small_disjunctive_matched %>% filter(treatment==1) %>% 
                                      dplyr::select(Death) %>% pull()
y_con = rhc_small_disjunctive_matched %>% filter(treatment==0) %>%
                                      dplyr::select(Death)%>% pull()
table(y_trt,y_con)
mcnemar.test(matrix(table(y_trt,y_con),2,2))

Propensity score matching

• Estimer le score propension grâce à une régression logistique et calculer les scores de propension.

Correction

psmodel = glm(treatment ~ . - Death ,
              family=binomial(),data=rhc_small_disjunctive)

summary(psmodel)
#create propensity scores
pscore = predict(psmodel,type="response")

• Représenter les distributions des scores de propension chez les contrôles et les traités.

Correction

ggplot(rhc_small_disjunctive %>% mutate(propensity = pscore) ,
       aes(as.factor(treatment),pscore,fill = treatment) )+ geom_violin()

• Calculer le logit des scores de propension, puis leur écart-type et le “caliper”

Correction

logit = function(p) {log(p)-log(1-p)}
logit_pscore = logit(pscore) # or predict(psmodel)
caliper = 0.2 *sd(logit_pscore)

• Appliquer le greedy matching sur le logitdu score de propension avec le “capiler” calculé à la question précédente

Correction

psmatch = Match(Tr=rhc_small_disjunctive$treatment,M=1,X=logit_pscore,
             replace=FALSE,caliper=caliper)
n_matched = length(psmatch$index.control)

rhc_small_disjunctive_matched_propensity =rhc_small_disjunctive %>%
                                  slice(c(psmatch$index.control,
                                          psmatch$index.treated)) %>%
                                  mutate( match_id = rep(c(1:n_matched),2))

• Calculer la “Table 1”

Correction

matchedtab1<-CreateTableOne(vars=xvars, strata ="treatment", 
                            data=rhc_small_disjunctive_matched_propensity, test = FALSE)
print(matchedtab1, smd = TRUE)

• Sur les données “matchées” étudiez l’outcome

Correction

y_trt = rhc_small_disjunctive_matched_propensity %>% filter(treatment==1) %>% 
                                      dplyr::select(Death) %>% pull()
y_con = rhc_small_disjunctive_matched_propensity %>% filter(treatment==0) %>%
                                      dplyr::select(Death)%>% pull()
table(y_trt,y_con)
mcnemar.test(matrix(table(y_trt,y_con),2,2))

• Faire varier la valeur du “caliper” et observer les changements dans la “Table 1” et le test sur l’outcome

Correction

for (cal_const in c(0.0001,0.01,0.1,0.5,1,5)){
  
  caliper = cal_const *sd(logit_pscore)
  print(paste(c("Le caliper vaut ",caliper),collapse=""))
  psmatch = Match(Tr=rhc_small_disjunctive$treatment,M=1,X=logit_pscore,
               replace=FALSE,caliper=caliper)
  n_matched = length(psmatch$index.control)
  
  rhc_small_disjunctive_matched_propensity =rhc_small_disjunctive %>%
                                    slice(c(psmatch$index.control,
                                            psmatch$index.treated)) %>%
                                    mutate( match_id = rep(c(1:n_matched),2))
  print(paste(c("Il y a ",length(psmatch$index.treated)," données matchées"),collapse=""))
  matchedtab1<-CreateTableOne(vars=xvars, strata ="treatment", 
                              data=rhc_small_disjunctive_matched_propensity, test = FALSE)
  print(matchedtab1, smd = TRUE)
  y_trt = rhc_small_disjunctive_matched_propensity %>% filter(treatment==1) %>% 
                                        dplyr::select(Death) %>% pull()
  y_con = rhc_small_disjunctive_matched_propensity %>% filter(treatment==0) %>%
                                        dplyr::select(Death)%>% pull()
  print(table(y_trt,y_con))
  print(mcnemar.test(matrix(table(y_trt,y_con),2,2)))
  }