library(tidyverse)
library(tableone) #create 'Table 1' to describe baseline characteristics
library(ipw) #inverse probability weighting
library(sandwich) #for robust variance estimation
library(survey)

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

Questions :

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.
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)
  • Calculer les scores de propension, puis les poids IPTW

Correction

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)))
  • Réprésenter la distribution des poids

Correction

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()
  • Calculer la moyenne pondérée pour l’âge dans le groupe traité, puis dans le groupe non-traité

Correction

#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)))
  • Calculer la Table 1 avec la pondération IPTW

Correction

#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)
  • Estimer le log odds ratio causal via une régression quasi-binomiale pondérée par les poids

Correction

#get causal risk difference
glm.obj<-glm(Death~treatment,weights=weight,family=quasibinomial,data = rhc_small_disjunctive)
summary(glm.obj)
betaiptw<-coef(glm.obj)
  • Calculer un estimateur sandwich et un intervalle de confiance pour le log odds ratio causal

Correction

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)
  • Faire de même pour l’odds ratio causal

Correction

#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)
  • Nous avions des poids très élévés, appliquer la méthode de troncation (avec un maximum fixé à 10) pour calculer le log-odds ratio causal

Correction

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)
  • Refaire la même étude avec le package ipw

Correction

#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)
LS0tCnRpdGxlOiAiRXhlbXBsZXMgcG91ciBsYSBjYXVzYWxpdMOpIDogSVBUVyBldCBNU00iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KYGBge3J9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KHRhYmxlb25lKSAjY3JlYXRlICdUYWJsZSAxJyB0byBkZXNjcmliZSBiYXNlbGluZSBjaGFyYWN0ZXJpc3RpY3MKbGlicmFyeShpcHcpICNpbnZlcnNlIHByb2JhYmlsaXR5IHdlaWdodGluZwpsaWJyYXJ5KHNhbmR3aWNoKSAjZm9yIHJvYnVzdCB2YXJpYW5jZSBlc3RpbWF0aW9uCmxpYnJhcnkoc3VydmV5KQpgYGAKCiMgRGF0YXNldCA6IHJpZ2h0IGhlYXJ0IGNhdGhldGVyaXphdGlvbgpVbmUgZGVzY3JpcHRpb24gY29tcGzDqHRlIGVzdCBkaXNwb25pYmxlIMOgIGNldHRlIGFkcmVzc2UgaHR0cDovL2Jpb3N0YXQubWMudmFuZGVyYmlsdC5lZHUvd2lraS9wdWIvTWFpbi9EYXRhU2V0cy9yaGMuaHRtbC4gSWwgcydhZ2l0IGRlIGRvbm7DqWVzIHN1ciBkZXMgcGF0aWVudHMgKDIxODQgdHJhaXTDqXMgZXQgMzU1MSBjb250csO0bGVzKSBhZG1pcyBhdXggdXJnZW5jZXMgZGFucyA1IGjDtHBpdGF1eCwgbGEgdmFyaWFibGUgZGUgdHJhaXRlbWVudCBlc3QgYHN3YW5nMWAgKHJpZ2h0IGhlYXQgY2F0aGV0ZXJpemF0aW9uIHZzLiBub24pIGV0IGwnb3V0Y29tZSBlc3QgYGRlYXRoYCAoeWVzIG9yIG5vKS4gTm91cyBhbGxvbnMgY29uc2lkw6lyZXIgbGVzIHZhcmlhYmxlcyBkZSBjb25mdXNpb24gc3VpdmFudGVzIDoKLSBgY2F0MWAgOiBQcmltYXJ5IGRpc2Vhc2UgY2F0ZWdvcnkKLSBgYWdlYAotIGBzZXhgYAotIGBtZWFuYnAxYDogTWVhbiBibG9vZCBwcmVzc3VyZQoKIyMgUXVlc3Rpb25zIDogCi0gQ2hhcmdlciBsZXMgZG9ubsOpZXMgaHR0cDovL2Jpb3N0YXQubWMudmFuZGVyYmlsdC5lZHUvd2lraS9NYWluL0RhdGFTZXRzLCAKYGBge3J9CnJoYyA9IHJlYWRfY3N2KCIuL0RhdGEvcmhjX2RhdGEuY3N2IikKcmhjX3NtYWxsID0gcmhjICU+JSBtdXRhdGUodHJlYXRtZW50ID0gc3dhbmcxICkgJT4lIAogICAgICAgICAgICAgICAgZHBseXI6OnNlbGVjdCh0cmVhdG1lbnQsZGVhdGgsIGNhdDEsIGFnZSxzZXgsbWVhbmJwMSwtc3dhbmcxKQogICAgICAgICAgICAgICAgCmdsaW1wc2UocmhjX3NtYWxsKQp1bmlxdWUocmhjX3NtYWxsJGNhdDEpCmBgYAotIENyw6llciBsYSAiVGFibGUgMSIuIEF0dGVudGlvbiBsZSBwYWNrYWdlIGB0YWJsZW9uZWAgbmUgc3VwcG9ydGUgcGFzIGxlcyB2YXJpYWJsZXMgZGlzY3LDqHRlcyDDoCBwbHVzaWV1cnMgZmFjdGV1cnMsIGlsIGZhdXQgZG9uYyBkJ2Fib3JkIGNvbnN0cnVpcmUgbGEgdGFibGUgZGlzam9uY3RpdmUuCmBgYHtyfQpyaGNfc21hbGxfZGlzanVuY3RpdmUgPSBhc190aWJibGUobW9kZWwubWF0cml4KH4gLiAsIGRhdGE9cmhjX3NtYWxsKVssLTFdKQpnbGltcHNlKHJoY19zbWFsbF9kaXNqdW5jdGl2ZSkKbmFtZXMocmhjX3NtYWxsX2Rpc2p1bmN0aXZlKSA9IGMoInRyZWF0bWVudCIsIkRlYXRoIiwiQ0hGIiwiQ2lyciIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJjb2xjYW4iLCJDb21hIiwiQ09QRCIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJsdW5nY2FuIiwiTU9TRiIsInNlcHNpcyIsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICJhZ2UiLCJtYWxlIiwibWVhbmJwMSIpCmBgYAoKYGBge3J9Cnh2YXJzID0gYygiQ0hGIiwiQ2lyciIsImNvbGNhbiIsIkNvbWEiLCJDT1BEIiwKICAgICAgICAgImx1bmdjYW4iLCJNT1NGIiwic2Vwc2lzIiwKICAgICAgICAgImFnZSIsIm1hbGUiLCJtZWFuYnAxIikgIyB2YXJpYWJsZXMgZGUgY29uZnVzaW9uCnRhYmxlMSA9IENyZWF0ZVRhYmxlT25lKHZhcnM9eHZhcnMsc3RyYXRhPSJ0cmVhdG1lbnQiLCBkYXRhPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSwgdGVzdD1GQUxTRSkKcHJpbnQodGFibGUxLHNtZD1UUlVFKQpgYGAKCi0gQ2FsY3VsZXIgbGVzIHNjb3JlcyBkZSBwcm9wZW5zaW9uLCBwdWlzIGxlcyBwb2lkcyBJUFRXCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KcHNtb2RlbCA9IGdsbSh0cmVhdG1lbnQgfiAuIC0gRGVhdGggLAogICAgICAgICAgICAgIGZhbWlseT1iaW5vbWlhbCgpLGRhdGE9cmhjX3NtYWxsX2Rpc2p1bmN0aXZlKQoKc3VtbWFyeShwc21vZGVsKQojY3JlYXRlIHByb3BlbnNpdHkgc2NvcmVzCnBzY29yZSA9IHByZWRpY3QocHNtb2RlbCx0eXBlPSJyZXNwb25zZSIpCgpyaGNfc21hbGxfZGlzanVuY3RpdmUgPSByaGNfc21hbGxfZGlzanVuY3RpdmUgJT4lIG11dGF0ZShwc2NvcmUgPSBwc2NvcmUpCgojY3JlYXRlIHdlaWdodHMKcmhjX3NtYWxsX2Rpc2p1bmN0aXZlID0gcmhjX3NtYWxsX2Rpc2p1bmN0aXZlICU+JSBtdXRhdGUod2VpZ2h0ID0gaWZlbHNlKHRyZWF0bWVudD09MSwxLyhwc2NvcmUpLDEvKDEtcHNjb3JlKSkpCmBgYAotIFLDqXByw6lzZW50ZXIgbGEgZGlzdHJpYnV0aW9uIGRlcyBwb2lkcwpgYGB7cn0KCmBgYAojIyMjIENvcnJlY3Rpb24KYGBge3J9CmdncGxvdChyaGNfc21hbGxfZGlzanVuY3RpdmUsIGFlcyh3ZWlnaHQpKSArIGdlb21fZGVuc2l0eSgpCmdncGxvdChyaGNfc21hbGxfZGlzanVuY3RpdmUsIGFlcyh4PWMoMTpucm93KHJoY19zbWFsbF9kaXNqdW5jdGl2ZSkpLHk9c29ydCh3ZWlnaHQpKSkgKyBnZW9tX3BvaW50KCkKYGBgCi0gQ2FsY3VsZXIgbGEgbW95ZW5uZSBwb25kw6lyw6llIHBvdXIgbCfDomdlIGRhbnMgbGUgZ3JvdXBlIHRyYWl0w6ksIHB1aXMgZGFucyBsZSBncm91cGUgbm9uLXRyYWl0w6kKYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQojdG8gZ2V0IGEgd2VpZ2h0ZWQgbWVhbiBmb3IgYSBzaW5nbGUgY292YXJpYXRlIGRpcmVjdGx5Ogp0cmVhdGVkID1yaGNfc21hbGxfZGlzanVuY3RpdmUgJT4lIGZpbHRlcih0cmVhdG1lbnQ9PTEpCnByaW50KG1lYW4odHJlYXRlZCR3ZWlnaHQqdHJlYXRlZCRhZ2UpLyhtZWFuKHRyZWF0ZWQkd2VpZ2h0KSkpCmNvbnRyb2wgPXJoY19zbWFsbF9kaXNqdW5jdGl2ZSAlPiUgZmlsdGVyKHRyZWF0bWVudD09MCkKcHJpbnQobWVhbihjb250cm9sJHdlaWdodCpjb250cm9sJGFnZSkvKG1lYW4oY29udHJvbCR3ZWlnaHQpKSkKYGBgCgotIENhbGN1bGVyIGxhIFRhYmxlIDEgYXZlYyBsYSBwb25kw6lyYXRpb24gSVBUVwpgYGB7cn0KCmBgYAojIyMjIENvcnJlY3Rpb24KYGBge3J9CiNhcHBseSB3ZWlnaHRzIHRvIGRhdGEKd2VpZ2h0ZWRkYXRhPC1zdnlkZXNpZ24oaWRzID0gfiAxLCBkYXRhID1yaGNfc21hbGxfZGlzanVuY3RpdmUsIHdlaWdodHMgPSB+IHdlaWdodCkKCiN3ZWlnaHRlZCB0YWJsZSAxCndlaWdodGVkdGFibGUgPC1zdnlDcmVhdGVUYWJsZU9uZSh2YXJzID0geHZhcnMsIHN0cmF0YSA9ICJ0cmVhdG1lbnQiLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgZGF0YSA9IHdlaWdodGVkZGF0YSwgdGVzdCA9IEZBTFNFKQojIyBTaG93IHRhYmxlIHdpdGggU01ECnByaW50KHdlaWdodGVkdGFibGUsIHNtZCA9IFRSVUUpCgpgYGAKLSBFc3RpbWVyIGxlIGxvZyBvZGRzIHJhdGlvIGNhdXNhbCB2aWEgdW5lIHLDqWdyZXNzaW9uIHF1YXNpLWJpbm9taWFsZSBwb25kw6lyw6llIHBhciBsZXMgcG9pZHMKYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCmBgYHtyfQojZ2V0IGNhdXNhbCByaXNrIGRpZmZlcmVuY2UKZ2xtLm9iajwtZ2xtKERlYXRofnRyZWF0bWVudCx3ZWlnaHRzPXdlaWdodCxmYW1pbHk9cXVhc2liaW5vbWlhbCxkYXRhID0gcmhjX3NtYWxsX2Rpc2p1bmN0aXZlKQpzdW1tYXJ5KGdsbS5vYmopCmJldGFpcHR3PC1jb2VmKGdsbS5vYmopCmBgYAotIENhbGN1bGVyIHVuIGVzdGltYXRldXIgc2FuZHdpY2ggZXQgdW4gaW50ZXJ2YWxsZSBkZSBjb25maWFuY2UgcG91ciBsZSBsb2cgb2RkcyByYXRpbyBjYXVzYWwKYGBge3J9CgpgYGAKIyMjIyBDb3JyZWN0aW9uCgpgYGB7cn0KU0U8LXNxcnQoZGlhZyh2Y292SEMoZ2xtLm9iaiwgdHlwZT0iSEMwIikpKQoKY2F1c2FscmQ8LShiZXRhaXB0d1syXSkKbGNsPC0oYmV0YWlwdHdbMl0tMS45NipTRVsyXSkKdWNsPC0oYmV0YWlwdHdbMl0rMS45NipTRVsyXSkKYyhsY2wsY2F1c2FscmQsdWNsKQpgYGAKCi0gRmFpcmUgZGUgbcOqbWUgcG91ciBsJ29kZHMgcmF0aW8gY2F1c2FsCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KI2dldCBjYXVzYWwgcmVsYXRpdmUgcmlzay4gV2VpZ2h0ZWQgR0xNCmdsbS5vYmo8LWdsbShEZWF0aH50cmVhdG1lbnQsd2VpZ2h0cz13ZWlnaHQsZmFtaWx5PXF1YXNpYmlub21pYWwsZGF0YT1yaGNfc21hbGxfZGlzanVuY3RpdmUpCiNzdW1tYXJ5KGdsbS5vYmopCmJldGFpcHR3PC1jb2VmKGdsbS5vYmopCiN0byBwcm9wZXJseSBhY2NvdW50IGZvciB3ZWlnaHRpbmcsIHVzZSBhc3ltcHRvdGljIChzYW5kd2ljaCkgdmFyaWFuY2UKU0U8LXNxcnQoZGlhZyh2Y292SEMoZ2xtLm9iaiwgdHlwZT0iSEMwIikpKQoKI2dldCBwb2ludCBlc3RpbWF0ZSBhbmQgQ0kgZm9yIHJlbGF0aXZlIHJpc2sgKG5lZWQgdG8gZXhwb25lbnRpYXRlKQpjYXVzYWxycjwtZXhwKGJldGFpcHR3WzJdKQpsY2w8LWV4cChiZXRhaXB0d1syXS0xLjk2KlNFWzJdKQp1Y2w8LWV4cChiZXRhaXB0d1syXSsxLjk2KlNFWzJdKQpjKGxjbCxjYXVzYWxycix1Y2wpCmBgYAoKLSBOb3VzIGF2aW9ucyBkZXMgcG9pZHMgdHLDqHMgw6lsw6l2w6lzLCBhcHBsaXF1ZXIgbGEgbcOpdGhvZGUgZGUgdHJvbmNhdGlvbiAoYXZlYyB1biBtYXhpbXVtIGZpeMOpIMOgIDEwKSBwb3VyIGNhbGN1bGVyIGxlIGxvZy1vZGRzIHJhdGlvIGNhdXNhbApgYGB7cn0KCmBgYAojIyMjIENvcnJlY3Rpb24KYGBge3J9CnJoY19zbWFsbF9kaXNqdW5jdGl2ZSA9IHJoY19zbWFsbF9kaXNqdW5jdGl2ZSAlPiUgbXV0YXRlKHRydW5jYXRlZF93ZWlnaHQgPSBpZmVsc2Uod2VpZ2h0PjEwLDEwLHdlaWdodCkpCgojZ2V0IGNhdXNhbCByaXNrIGRpZmZlcmVuY2UKZ2xtLm9iajwtZ2xtKERlYXRofnRyZWF0bWVudCx3ZWlnaHRzPXRydW5jYXRlZF93ZWlnaHQsZmFtaWx5PXF1YXNpYmlub21pYWwsZGF0YSA9IHJoY19zbWFsbF9kaXNqdW5jdGl2ZSkKI3N1bW1hcnkoZ2xtLm9iaikKYmV0YWlwdHc8LWNvZWYoZ2xtLm9iaikKU0U8LXNxcnQoZGlhZyh2Y292SEMoZ2xtLm9iaiwgdHlwZT0iSEMwIikpKQoKY2F1c2FscmQ8LShiZXRhaXB0d1syXSkKbGNsPC0oYmV0YWlwdHdbMl0tMS45NipTRVsyXSkKdWNsPC0oYmV0YWlwdHdbMl0rMS45NipTRVsyXSkKYyhsY2wsY2F1c2FscmQsdWNsKQpgYGAKCi0gUmVmYWlyZSBsYSBtw6ptZSDDqXR1ZGUgYXZlYyBsZSBwYWNrYWdlIGBpcHdgCmBgYHtyfQoKYGBgCiMjIyMgQ29ycmVjdGlvbgpgYGB7cn0KI3Bsb3Qgb2Ygd2VpZ2h0cwppcHdwbG90KHdlaWdodHMgPSByaGNfc21hbGxfZGlzanVuY3RpdmUkd2VpZ2h0LCBsb2dzY2FsZSA9IEZBTFNFLAogICAgICAgICBtYWluID0gIndlaWdodHMiLCB4bGltID0gYygwLCAyMikpCmlwd3Bsb3Qod2VpZ2h0cyA9IHJoY19zbWFsbF9kaXNqdW5jdGl2ZSR0cnVuY2F0ZWRfd2VpZ2h0LCBsb2dzY2FsZSA9IEZBTFNFLAogICAgICAgICBtYWluID0gIndlaWdodHMiLCB4bGltID0gYygwLCAyMikpCiNmaXQgYSBtYXJnaW5hbCBzdHJ1Y3R1cmFsIG1vZGVsIChyaXNrIGRpZmZlcmVuY2UpCm1zbSA8LSBzdnlnbG0oRGVhdGggfiB0cmVhdG1lbnQsIAogICAgICAgICAgICAgIGRlc2lnbiA9IHN2eWRlc2lnbih+IDEsIHdlaWdodHMgPSB+IHJoY19zbWFsbF9kaXNqdW5jdGl2ZSR0cnVuY2F0ZWRfd2VpZ2h0LCBkYXRhID1yaGNfc21hbGxfZGlzanVuY3RpdmUpLAogICAgICAgICAgICAgIGZhbWlseT1zdGF0czo6cXVhc2liaW5vbWlhbCgpKQpjb2VmKG1zbSkKY29uZmludChtc20pCmBgYAoK