library(tidyverse)
library(tableone) #create 'Table 1' to describe baseline characteristics
library(Matching) #multivariate and propensity score matching with balance optimization
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
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)
tableone
ne supporte pas les variables discrètes à plusieurs facteurs, il faut donc d’abord construire la table disjonctive.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")
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)
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)
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))
psmodel = glm(treatment ~ . - Death ,
family=binomial(),data=rhc_small_disjunctive)
summary(psmodel)
#create propensity scores
pscore = predict(psmodel,type="response")
ggplot(rhc_small_disjunctive %>% mutate(propensity = pscore) ,
aes(as.factor(treatment),pscore,fill = treatment) )+ geom_violin()
logit = function(p) {log(p)-log(1-p)}
logit_pscore = logit(pscore) # or predict(psmodel)
caliper = 0.2 *sd(logit_pscore)
logit
du score de propension avec le “capiler” calculé à la question précédentepsmatch = 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))
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()
table(y_trt,y_con)
mcnemar.test(matrix(table(y_trt,y_con),2,2))
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)))
}