library(tidyverse)
library(tableone) #create 'Table 1' to describe baseline characteristics
library(ipw) #inverse probability weighting
library(sandwich) #for robust variance estimation
library(survey)
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)
psmodel = glm(treatment ~ . - Death ,
family=binomial(),data=rhc_small_disjunctive)
summary(psmodel)
#create propensity scores
pscore = predict(psmodel,type="response")
rhc_small_disjunctive = rhc_small_disjunctive %>% mutate(pscore = pscore)
#create weights
rhc_small_disjunctive = rhc_small_disjunctive %>% mutate(weight = ifelse(treatment==1,1/(pscore),1/(1-pscore)))
ggplot(rhc_small_disjunctive, aes(weight)) + geom_density()
ggplot(rhc_small_disjunctive, aes(x=c(1:nrow(rhc_small_disjunctive)),y=sort(weight))) + geom_point()
#to get a weighted mean for a single covariate directly:
treated =rhc_small_disjunctive %>% filter(treatment==1)
print(mean(treated$weight*treated$age)/(mean(treated$weight)))
control =rhc_small_disjunctive %>% filter(treatment==0)
print(mean(control$weight*control$age)/(mean(control$weight)))
#apply weights to data
weighteddata<-svydesign(ids = ~ 1, data =rhc_small_disjunctive, weights = ~ weight)
#weighted table 1
weightedtable <-svyCreateTableOne(vars = xvars, strata = "treatment",
data = weighteddata, test = FALSE)
## Show table with SMD
print(weightedtable, smd = TRUE)
#get causal risk difference
glm.obj<-glm(Death~treatment,weights=weight,family=quasibinomial,data = rhc_small_disjunctive)
summary(glm.obj)
betaiptw<-coef(glm.obj)
SE<-sqrt(diag(vcovHC(glm.obj, type="HC0")))
causalrd<-(betaiptw[2])
lcl<-(betaiptw[2]-1.96*SE[2])
ucl<-(betaiptw[2]+1.96*SE[2])
c(lcl,causalrd,ucl)
#get causal relative risk. Weighted GLM
glm.obj<-glm(Death~treatment,weights=weight,family=quasibinomial,data=rhc_small_disjunctive)
#summary(glm.obj)
betaiptw<-coef(glm.obj)
#to properly account for weighting, use asymptotic (sandwich) variance
SE<-sqrt(diag(vcovHC(glm.obj, type="HC0")))
#get point estimate and CI for relative risk (need to exponentiate)
causalrr<-exp(betaiptw[2])
lcl<-exp(betaiptw[2]-1.96*SE[2])
ucl<-exp(betaiptw[2]+1.96*SE[2])
c(lcl,causalrr,ucl)
rhc_small_disjunctive = rhc_small_disjunctive %>% mutate(truncated_weight = ifelse(weight>10,10,weight))
#get causal risk difference
glm.obj<-glm(Death~treatment,weights=truncated_weight,family=quasibinomial,data = rhc_small_disjunctive)
#summary(glm.obj)
betaiptw<-coef(glm.obj)
SE<-sqrt(diag(vcovHC(glm.obj, type="HC0")))
causalrd<-(betaiptw[2])
lcl<-(betaiptw[2]-1.96*SE[2])
ucl<-(betaiptw[2]+1.96*SE[2])
c(lcl,causalrd,ucl)
ipw
#plot of weights
ipwplot(weights = rhc_small_disjunctive$weight, logscale = FALSE,
main = "weights", xlim = c(0, 22))
ipwplot(weights = rhc_small_disjunctive$truncated_weight, logscale = FALSE,
main = "weights", xlim = c(0, 22))
#fit a marginal structural model (risk difference)
msm <- svyglm(Death ~ treatment,
design = svydesign(~ 1, weights = ~ rhc_small_disjunctive$truncated_weight, data =rhc_small_disjunctive),
family=stats::quasibinomial())
coef(msm)
confint(msm)