Causal Models : Propensity Score Methods
Estimating causal effects from observational studies using the propensity score approach
Fist of all, we will upload the required libraries:
library(twang)
library(magrittr)
library(tidyverse)
library(gtsummary)
library(stddiff)
library(ggplot2)
library(data.table)
library(boot)
library(splines)
library(PSAgraphics)
library(Matching)
library(sandwich)
library(survey)
library(rms)We will use a dataset coming from an observational study of 996 patients receiving an initial Percutaneous Coronary Intervention (PCI) at Ohio Heart Health, Christ Hospital, Cincinnati in 1997 and followed for at least 6 months by the staff of the Lindner Center.
The patients thought to be more severely diseased were assigned to treatment with abciximab (an expensive, high-molecular-weight IIb/IIIa cascade blocker); in fact, only 298 (29.9 percent) of patients received usual-care-alone with their initial PCI.
Our research question aims at estimating the treatment effect of abciximab+PCI (abcix) vs the standard care on the probability of of being deceased at 6 months.
In this practical, we will apply the methods based on the propensity score.
Measured pre-treatment characteristics that could confound the treatment-outcome relationship are:
- acutemi: Recent acute myocardial infarction ( 0 No, 1 Yes)
- ejecfrac: Left Ejection Fraction (%)
- ves1proc: number of vessels involved ( from 1 to 5)
- stent: 1 indicates coronary stent inserted
- diabetic: 1 indicates the subject is diabetic
- height: height of the subject in cm
- female: 1 indicates female subjects
Exploring the data
data("lindner",package="twang")
set.seed(123)summary(lindner)## lifepres cardbill abcix stent
## Min. : 0.0 Min. : 2216 Min. :0.0000 Min. :0.0000
## 1st Qu.:11.6 1st Qu.: 10219 1st Qu.:0.0000 1st Qu.:0.0000
## Median :11.6 Median : 12458 Median :1.0000 Median :1.0000
## Mean :11.3 Mean : 15674 Mean :0.7008 Mean :0.6687
## 3rd Qu.:11.6 3rd Qu.: 16660 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :11.6 Max. :178534 Max. :1.0000 Max. :1.0000
## height female diabetic acutemi
## Min. :108.0 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:165.0 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000
## Median :173.0 Median :0.0000 Median :0.0000 Median :0.0000
## Mean :171.4 Mean :0.3474 Mean :0.2239 Mean :0.1436
## 3rd Qu.:178.0 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.0000
## Max. :196.0 Max. :1.0000 Max. :1.0000 Max. :1.0000
## ejecfrac ves1proc sixMonthSurvive
## Min. : 0.00 Min. :0.000 Mode :logical
## 1st Qu.:45.00 1st Qu.:1.000 FALSE:26
## Median :55.00 Median :1.000 TRUE :970
## Mean :50.97 Mean :1.386
## 3rd Qu.:56.00 3rd Qu.:2.000
## Max. :90.00 Max. :5.000How is the treatment variable distributed in the population?
# Exposure
lindner %>%
dplyr::select(abcix) %>%
tbl_summary()| Characteristic | N = 9961 |
|---|---|
| abcix | 698 (70%) |
| 1 n (%) | |
ggplot(lindner)+
geom_bar(aes(x=abcix,fill=as.factor(abcix)),stat="count")+
scale_fill_discrete("Treatment",labels = c("PCI", "PCI+abciximab")) +
theme_classic()
Is the outcome rare?
# Outcome
lindner$sixMonthDeath <- 1-lindner$sixMonthSurvive
lindner %>%
dplyr::select(sixMonthDeath) %>%
tbl_summary()| Characteristic | N = 9961 |
|---|---|
| sixMonthDeath | 26 (2.6%) |
| 1 n (%) | |
What is the crude odds ratio for mortality?
# Crude OR
fit.crude <- glm(sixMonthDeath~abcix,family = binomial,data=lindner)
tbl_regression(fit.crude,exponentiate=T)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| abcix | 0.30 | 0.13, 0.66 | 0.003 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
Let’s now summarize the confounding variables by treatment group and by outcome, to have a general idea about the observed associations:
# Possible confounders
lindner %>%
dplyr::select(acutemi,
ejecfrac,
ves1proc,
stent,
diabetic,
female,
height) %>%
tbl_summary()| Characteristic | N = 9961 |
|---|---|
| acutemi | 143 (14%) |
| ejecfrac | 55 (45, 56) |
| ves1proc | |
| 0 | 4 (0.4%) |
| 1 | 680 (68%) |
| 2 | 252 (25%) |
| 3 | 45 (4.5%) |
| 4 | 14 (1.4%) |
| 5 | 1 (0.1%) |
| stent | 666 (67%) |
| diabetic | 223 (22%) |
| female | 346 (35%) |
| height | 173 (165, 178) |
| 1 n (%); Median (IQR) | |
# Descriptive statistics of patients'characteristics by treatment group
lindner %>%
dplyr::select(acutemi,
ejecfrac,
ves1proc,
stent,
diabetic,
female,
height,
abcix) %>%
tbl_summary(by=abcix) %>%
add_overall() %>%
add_p()| Characteristic | Overall, N = 9961 | 0, N = 2981 | 1, N = 6981 | p-value2 |
|---|---|---|---|---|
| acutemi | 143 (14%) | 18 (6.0%) | 125 (18%) | <0.001 |
| ejecfrac | 55 (45, 56) | 55 (50, 60) | 55 (45, 56) | 0.004 |
| ves1proc | <0.001 | |||
| 0 | 4 (0.4%) | 1 (0.3%) | 3 (0.4%) | |
| 1 | 680 (68%) | 243 (82%) | 437 (63%) | |
| 2 | 252 (25%) | 47 (16%) | 205 (29%) | |
| 3 | 45 (4.5%) | 6 (2.0%) | 39 (5.6%) | |
| 4 | 14 (1.4%) | 1 (0.3%) | 13 (1.9%) | |
| 5 | 1 (0.1%) | 0 (0%) | 1 (0.1%) | |
| stent | 666 (67%) | 174 (58%) | 492 (70%) | <0.001 |
| diabetic | 223 (22%) | 80 (27%) | 143 (20%) | 0.027 |
| female | 346 (35%) | 115 (39%) | 231 (33%) | 0.10 |
| height | 173 (165, 178) | 173 (165, 178) | 173 (165, 180) | 0.8 |
| 1 n (%); Median (IQR) | ||||
| 2 Pearson’s Chi-squared test; Wilcoxon rank sum test; Fisher’s exact test | ||||
# Descriptive statistics of patients'characteristics by outcome
lindner %>%
dplyr::select(acutemi,
ejecfrac,
ves1proc,
stent,
diabetic,
female,
height,
sixMonthDeath) %>%
tbl_summary(by=sixMonthDeath) %>%
add_overall() %>%
add_p()| Characteristic | Overall, N = 9961 | 0, N = 9701 | 1, N = 261 | p-value2 |
|---|---|---|---|---|
| acutemi | 143 (14%) | 137 (14%) | 6 (23%) | 0.2 |
| ejecfrac | 55 (45, 56) | 55 (45, 56) | 50 (33, 54) | 0.002 |
| ves1proc | 0.10 | |||
| 0 | 4 (0.4%) | 4 (0.4%) | 0 (0%) | |
| 1 | 680 (68%) | 661 (68%) | 19 (73%) | |
| 2 | 252 (25%) | 249 (26%) | 3 (12%) | |
| 3 | 45 (4.5%) | 41 (4.2%) | 4 (15%) | |
| 4 | 14 (1.4%) | 14 (1.4%) | 0 (0%) | |
| 5 | 1 (0.1%) | 1 (0.1%) | 0 (0%) | |
| stent | 666 (67%) | 649 (67%) | 17 (65%) | 0.9 |
| diabetic | 223 (22%) | 212 (22%) | 11 (42%) | 0.014 |
| female | 346 (35%) | 335 (35%) | 11 (42%) | 0.4 |
| height | 173 (165, 178) | 173 (165, 178) | 168 (163, 177) | 0.10 |
| 1 n (%); Median (IQR) | ||||
| 2 Fisher’s exact test; Wilcoxon rank sum test; Pearson’s Chi-squared test | ||||
We can now calculate the Standardized Difference,which can be use as a measure of balance in the treatment groups. It is a measure of difference between groups that is independent from statistical testing (remember that p values always depend on sample size !!). It is very similar to the definition of effect size that we discussed in Block 2.
It can be defined for a continuous covariate as:
\(SD_{c}=\frac{\overline{x_{1}}-\overline{x_{0}}} {\sqrt{\frac{s_1^2+s_0^2}{2}}}\)
and for a dichotomous covariate as:
\(SD_{d}=\frac{\overline{p_{1}}-\overline{p_{0}}} {\sqrt{\frac{\overline{p_{1}}(1-\overline{p_{1})}+\overline{p_{0}}(1-\overline{p_{0}})}{2}}}\)
The rough interpretation is that imbalance is present if the standardized difference is greater than 0.1 or 0.2.
s1 <- stddiff.numeric(vcol="height",gcol="abcix",data=lindner)
s2 <- stddiff.numeric(vcol="ejecfrac",gcol="abcix",data=lindner)
s3 <- stddiff.numeric(vcol="ves1proc",gcol="abcix",data=lindner)
s4 <- stddiff.binary(vcol="stent",gcol="abcix",data=lindner)
s5 <- stddiff.binary(vcol="female",gcol="abcix",data=lindner)
s6 <- stddiff.binary(vcol="diabetic",gcol="abcix",data=lindner)
s7 <- stddiff.binary(vcol="acutemi",gcol="abcix",data=lindner)
cont.var <- as.data.frame(rbind(s1,s2,s3))
rownames(cont.var) <- c("height", "ejecfrac", "ves1proc")
cont.var## mean.c sd.c mean.t sd.t missing.c missing.t stddiff stddiff.l
## height 171.446 10.589 171.443 10.695 0 0 0.000 -0.135
## ejecfrac 52.289 10.297 50.403 10.419 0 0 0.182 0.046
## ves1proc 1.205 0.480 1.463 0.706 0 0 0.427 0.290
## stddiff.u
## height 0.136
## ejecfrac 0.318
## ves1proc 0.564bin.var <- as.data.frame(rbind(s4,s5,s6, s7))
rownames(bin.var) <- c("stent", "female", "diabetic", "acutemi")
bin.var## p.c p.t missing.c missing.t stddiff stddiff.l stddiff.u
## stent 0.584 0.705 0 0 0.255 0.119 0.391
## female 0.386 0.331 0 0 0.115 -0.021 0.251
## diabetic 0.268 0.205 0 0 0.150 0.014 0.286
## acutemi 0.060 0.179 0 0 0.372 0.235 0.508Estimating the propensity score
Fit now a propensity score model — a logistic regression model with abciximab+PCI (vs. PCI) as the outcome, and the confounders listed in the table above included as covariates. We exclude from the list the variable height, since there was not a relevant difference between the groups.
# Fit a propensity score model
fit.ps<- glm(abcix~
acutemi+
ejecfrac+
ves1proc+
stent+
diabetic+
female,
data=lindner,family = binomial)
summary(fit.ps)##
## Call:
## glm(formula = abcix ~ acutemi + ejecfrac + ves1proc + stent +
## diabetic + female, family = binomial, data = lindner)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.5479 -1.2311 0.6309 0.8799 1.6053
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.272481 0.440711 0.618 0.536394
## acutemi 1.178158 0.269479 4.372 1.23e-05 ***
## ejecfrac -0.015076 0.007391 -2.040 0.041374 *
## ves1proc 0.757365 0.138279 5.477 4.32e-08 ***
## stent 0.571165 0.150208 3.803 0.000143 ***
## diabetic -0.409191 0.170371 -2.402 0.016316 *
## female -0.133346 0.151377 -0.881 0.378379
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1215.5 on 995 degrees of freedom
## Residual deviance: 1127.0 on 989 degrees of freedom
## AIC: 1141
##
## Number of Fisher Scoring iterations: 4# Save the estimated propensity score
lindner$ps <- fitted(fit.ps)# Plot estimated ps
ggplot(lindner) +
geom_boxplot(aes(y = ps,group = as.factor(abcix),col = as.factor(abcix))) +
scale_y_continuous("Estimated PS") +
scale_color_discrete("Treatment",labels = c("PCI", "PCI+abciximab")) +
theme_classic()
ggplot(lindner) +
geom_histogram(aes(x = ps,group = as.factor(abcix),fill = as.factor(abcix))) +
scale_y_continuous("Estimated PS") +
facet_grid(cols=vars(abcix))+
scale_fill_discrete("Treatment",labels = c("PCI", "PCI+abciximab")) +
theme_classic()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Assess whether there are non-overlapping scores (positivity violation) in the two exposure groups:
lindner %>%
dplyr::select(ps,abcix) %>%
tbl_summary(type = list(ps~"continuous2"),by=abcix,
statistic = all_continuous2() ~c(
"{median} ({p25}, {p75})",
"{min}, {max}"))| Characteristic | 0, N = 298 | 1, N = 698 |
|---|---|---|
| ps | ||
| Median (IQR) | 0.65 (0.55, 0.70) | 0.71 (0.65, 0.82) |
| Range | 0.25, 0.96 | 0.28, 0.98 |
Investigate overlap:
lindner %<>% mutate(overlap=ifelse(ps>=min(ps[abcix==1]) &
ps<=max(ps[abcix==0]),1,0))
# non-overlap: treatment group have higher ps than any non-abciximab user and
# control group have smaller ps than any abciximab user
with(lindner,table(overlap,abcix))## abcix
## overlap 0 1
## 0 1 8
## 1 297 690with(lindner,prop.table(table(overlap,abcix)),2)## abcix
## overlap 0 1
## 0 0.001004016 0.008032129
## 1 0.298192771 0.692771084In the successive steps, we remove subjects that does not overlap. This step reduce the original sample size, but we should respect the assumption of positivity in order to estimate a reasonable causal effect. It makes no sense including subjects that have “near-zero” probability to receive the treatment or to have a “match” in the successive analyses.
First option: “adjusting” for the propensity score
We can use the estimated PS as a covariate in a logistic regression model for the outcome:
# Model 1: Linear relationship between ps and outcome
fit.out<- glm(sixMonthDeath~abcix+ps,
data=lindner,
family = binomial,
subset=overlap==1)
# Model summary
summary(fit.out)##
## Call:
## glm(formula = sixMonthDeath ~ abcix + ps, family = binomial,
## data = lindner, subset = overlap == 1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4229 -0.2821 -0.1833 -0.1596 3.1184
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.211 1.165 -3.616 0.000300 ***
## abcix -1.445 0.439 -3.293 0.000992 ***
## ps 1.947 1.698 1.147 0.251522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 233.15 on 986 degrees of freedom
## Residual deviance: 222.01 on 984 degrees of freedom
## AIC: 228.01
##
## Number of Fisher Scoring iterations: 7The second step is to save the predicted probabilities for the treated and the untreated and estimate the causal effects of interest in the population:
# fitted values (probabilities)
lindner$predY0<-fit.out$family$linkinv(coef(fit.out)[1]+coef(fit.out)[3]*lindner$ps) # PCI subjects
lindner$predY1<-fit.out$family$linkinv(coef(fit.out)[1]+coef(fit.out)[2]+coef(fit.out)[3]*lindner$ps) # PCI+abciximab subjects
# ATE effect
Y1<-mean(lindner$predY1)
Y0<-mean(lindner$predY0)
# ATT effect
Y1_1<-mean(lindner$predY1[lindner$abcix==1])
Y0_1<-mean(lindner$predY0[lindner$abcix==1])
# ATE effect
Y1-Y0## [1] -0.04246731# ATT effect
Y1_1-Y0_1## [1] -0.04436082# Estimate odds ratios related to the "ATE" and the "ATT"
(Y1/(1-Y1))/(Y0/(1-Y0))## [1] 0.2363187(Y1_1/(1-Y1_1))/(Y0_1/(1-Y0_1))## [1] 0.236305The ATE effect is quite similar to the ATT effect, indicating that there is a protective effect of the PCI+abciximab vs PCI alone.
To obtain the corresponding confidence intervals we can use the bootstrap approach.
We do not here outline this procedure, see at the end of this practical the supplementary code.
This model relies on two additional assumptions: no interaction between propensity score and treatment, and a linear relationship between the propensity score and treatment.
Do these assumptions appear reasonable here?
We can try to fit different models, and then compare the AIC:
# Model 2: Non-linear relationship between ps and outcome
fit.out2<- glm(sixMonthDeath~abcix+ps+I(ps^2),
data=lindner,
family = binomial,
subset=overlap==1)
summary(fit.out2)##
## Call:
## glm(formula = sixMonthDeath ~ abcix + ps + I(ps^2), family = binomial,
## data = lindner, subset = overlap == 1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5751 -0.2881 -0.1730 -0.1407 3.0416
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.3787 3.8937 0.611 0.541255
## abcix -1.5088 0.4493 -3.358 0.000784 ***
## ps -17.7559 11.5272 -1.540 0.123475
## I(ps^2) 14.1968 8.3131 1.708 0.087681 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 233.15 on 986 degrees of freedom
## Residual deviance: 219.46 on 983 degrees of freedom
## AIC: 227.46
##
## Number of Fisher Scoring iterations: 7# Model 3: Non-linear relationship between ps and outcome and interaction between ps and treatment
fit.out3<- glm(sixMonthDeath~abcix*ps+I(ps^2),
data=lindner,
family = binomial,
subset=overlap==1)
summary(fit.out3)##
## Call:
## glm(formula = sixMonthDeath ~ abcix * ps + I(ps^2), family = binomial,
## data = lindner, subset = overlap == 1)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.6659 -0.2798 -0.1697 -0.1458 3.0262
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.361196 3.786408 0.624 0.5329
## abcix 0.009567 2.075020 0.005 0.9963
## ps -18.641792 11.183815 -1.667 0.0955 .
## I(ps^2) 15.513326 8.167099 1.899 0.0575 .
## abcix:ps -2.130967 2.862505 -0.744 0.4566
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 233.15 on 986 degrees of freedom
## Residual deviance: 218.93 on 982 degrees of freedom
## AIC: 228.93
##
## Number of Fisher Scoring iterations: 7It seems that the relationship could be partially non-linear, but there is no a statistical significance very strong, as well as for the interaction. So probably the best parsimonious model to keep is Model 1.
Second option: Stratification
Create propensity score strata: this could be an iterative process, since we should verify if we have enough subjects/event in each stratum.
lindner %<>% mutate(strata=cut(ps,quantile(ps,c(0,0.25,0.5,0.75,1)),include.lowest=T,labels=c(1:4))) #Check they have been created correctly
summary(lindner$strata)## 1 2 3 4
## 252 246 257 241tapply(lindner$ps,lindner$strata,summary)## $`1`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.2499 0.5168 0.5463 0.5321 0.5713 0.6051
##
## $`2`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6078 0.6505 0.6674 0.6605 0.6807 0.6846
##
## $`3`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.6877 0.7036 0.7307 0.7406 0.7728 0.8106
##
## $`4`
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.8106 0.8305 0.8706 0.8759 0.9102 0.9765#Look at numbers of events and patients in each strata/exposure group
table(lindner$sixMonthDeath,lindner$strata,lindner$abcix)## , , = 0
##
##
## 1 2 3 4
## 0 105 89 64 25
## 1 7 1 4 3
##
## , , = 1
##
##
## 1 2 3 4
## 0 138 156 185 208
## 1 2 0 4 5These strata seem quite “sparse” as number of events. Another possibility is :
#Create propensity score strata
lindner %<>% mutate(strata=cut(ps,quantile(ps,c(0,0.33,0.66,1)),include.lowest=T,labels=c(1:3)))
#Check they have been created correctly
summary(lindner$strata)## 1 2 3
## 356 302 338#Look at numbers of events and patients in each strata/exposure group
table(lindner$sixMonthDeath,lindner$strata,lindner$abcix)## , , = 0
##
##
## 1 2 3
## 0 146 96 41
## 1 7 4 4
##
## , , = 1
##
##
## 1 2 3
## 0 201 201 285
## 1 2 1 8Remind : we should also check for the balance of the confounders in each strata ! See the supplementary material for that. For now, let’s just estimate the OR in each stratum:
beta.treat<-numeric(3)
nstrata<-table(lindner$strata)
treated.strata<-table(lindner$strata,lindner$abcix)[,2]
for (i in 1:3){
ms<-glm(sixMonthDeath~abcix,data=lindner,subset = strata==i,family="binomial")
beta.treat[i]<-coef(ms)[2]
print(summary(ms))
}##
## Call:
## glm(formula = sixMonthDeath ~ abcix, family = "binomial", data = lindner,
## subset = strata == i)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.3060 -0.3060 -0.1407 -0.1407 3.0398
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0377 0.3869 -7.851 4.13e-15 ***
## abcix -1.5725 0.8091 -1.943 0.052 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 83.969 on 355 degrees of freedom
## Residual deviance: 79.319 on 354 degrees of freedom
## AIC: 83.319
##
## Number of Fisher Scoring iterations: 7
##
##
## Call:
## glm(formula = sixMonthDeath ~ abcix, family = "binomial", data = lindner,
## subset = strata == i)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.2857 -0.2857 -0.0996 -0.0996 3.2583
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.1781 0.5103 -6.228 4.73e-10 ***
## abcix -2.1253 1.1249 -1.889 0.0589 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 50.927 on 301 degrees of freedom
## Residual deviance: 46.200 on 300 degrees of freedom
## AIC: 50.2
##
## Number of Fisher Scoring iterations: 8
##
##
## Call:
## glm(formula = sixMonthDeath ~ abcix, family = "binomial", data = lindner,
## subset = strata == i)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4315 -0.2353 -0.2353 -0.2353 2.6835
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3273 0.5238 -4.443 8.88e-06 ***
## abcix -1.2458 0.6347 -1.963 0.0497 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 103.68 on 337 degrees of freedom
## Residual deviance: 100.39 on 336 degrees of freedom
## AIC: 104.39
##
## Number of Fisher Scoring iterations: 6And, finally, let’s estimate the weighted OR related to the ATE and the ATT as a weighted average of the ORs in the various strata:
exp(sum(beta.treat*nstrata)/nrow(lindner))## [1] 0.1960846exp(sum(beta.treat*treated.strata)/sum(treated.strata))## [1] 0.2028472Also here, we should use a bootstrap approach to estimate the corresponding confidence intervals.
Third option: Matching
Here we should create a reduced dataset retaining only patients with the overlap:
lindner.overlap <- lindner %>% filter(overlap==1)Now we proceed with the matching algorithm: there is plenty of different algorithms in R that produce matching, here we use one from the library(Matching).
match <- Match(Y=lindner.overlap$sixMonthDeath,
Tr=lindner.overlap$abcix,
X=lindner.overlap$ps,
caliper=0.2,# all matches not equal to or within 0.2 standard deviations of ps are dropped
M=1,
ties=F,
replace=TRUE # 1:1
)# Number of pairs
nn <- length(match$index.treated)
# Create matched dataset
lindnerMatched <- cbind(rbind(lindner.overlap[match$index.treated,],
lindner.overlap[match$index.control,]),
pair=c(1:nn,1:nn))
table(lindner.overlap$abcix)##
## 0 1
## 297 690#Check number of treated patients
table(lindnerMatched$abcix)##
## 0 1
## 689 689#Look at people being used multiple times in the matched sample
summary(as.factor(table(match$index.treated)))## 1
## 689summary(as.factor(table(match$index.control)))## 1 2 3 4 5 6 7 8 9 10 12 13 15 17 19
## 69 45 30 20 11 9 9 5 2 2 1 1 2 1 2#Look at the propensity score distribution in the matched dataset
ggplot(lindnerMatched) +
geom_boxplot(aes(y = ps,group = as.factor(abcix),col = as.factor(abcix))) +
scale_y_continuous("Estimated PS") +
scale_color_discrete("Treatment",labels = c("PCI", "PCI+abciximab")) +
theme_classic()
ggplot(lindnerMatched) +
geom_histogram(aes(x = ps,group = as.factor(abcix),fill = as.factor(abcix))) +
scale_y_continuous("Estimated PS") +
facet_grid(cols=vars(abcix))+
scale_fill_discrete("Treatment",labels = c("PCI", "PCI+abciximab")) +
theme_classic()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Let’s check now the balance:
# Balance Diagnostics before and after matching
bal <- MatchBalance(abcix~
stent+
female+
diabetic+
acutemi+
ejecfrac+
ves1proc,
data=lindner.overlap,
match.out = match)##
## ***** (V1) stent *****
## Before Matching After Matching
## mean treatment........ 0.70435 0.70537
## mean control.......... 0.58586 0.65602
## std mean diff......... 25.947 10.817
##
## mean raw eQQ diff..... 0.11785 0.049347
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.059245 0.024673
## med eCDF diff........ 0.059245 0.024673
## max eCDF diff........ 0.11849 0.049347
##
## var ratio (Tr/Co)..... 0.85663 0.92097
## T-test p-value........ 0.0004398 0.0057027
##
##
## ***** (V2) female *****
## Before Matching After Matching
## mean treatment........ 0.33333 0.33382
## mean control.......... 0.38384 0.30479
## std mean diff......... -10.706 6.151
##
## mean raw eQQ diff..... 0.050505 0.029028
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.025253 0.014514
## med eCDF diff........ 0.025253 0.014514
## max eCDF diff........ 0.050505 0.029028
##
## var ratio (Tr/Co)..... 0.9378 1.0495
## T-test p-value........ 0.13211 0.09767
##
##
## ***** (V3) diabetic *****
## Before Matching After Matching
## mean treatment........ 0.20725 0.2061
## mean control.......... 0.26599 0.21771
## std mean diff......... -14.483 -2.8684
##
## mean raw eQQ diff..... 0.057239 0.011611
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.029373 0.0058055
## med eCDF diff........ 0.029373 0.0058055
## max eCDF diff........ 0.058747 0.011611
##
## var ratio (Tr/Co)..... 0.83988 0.96072
## T-test p-value........ 0.050489 0.47958
##
##
## ***** (V4) acutemi *****
## Before Matching After Matching
## mean treatment........ 0.17101 0.17126
## mean control.......... 0.060606 0.16255
## std mean diff......... 29.302 2.3098
##
## mean raw eQQ diff..... 0.11111 0.0087083
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.055204 0.0043541
## med eCDF diff........ 0.055204 0.0043541
## max eCDF diff........ 0.11041 0.0087083
##
## var ratio (Tr/Co)..... 2.4853 1.0426
## T-test p-value........ 4.1797e-08 0.5361
##
##
## ***** (V5) ejecfrac *****
## Before Matching After Matching
## mean treatment........ 50.559 50.553
## mean control.......... 52.279 51.209
## std mean diff......... -16.899 -6.4414
##
## mean raw eQQ diff..... 1.8687 0.80406
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 20 15
##
## mean eCDF diff........ 0.033253 0.013763
## med eCDF diff........ 0.0087396 0.0043541
## max eCDF diff........ 0.10771 0.046444
##
## var ratio (Tr/Co)..... 0.97403 1.0808
## T-test p-value........ 0.016162 0.13002
## KS Bootstrap p-value.. 0.002 0.216
## KS Naive p-value...... 0.016165 0.44722
## KS Statistic.......... 0.10771 0.046444
##
##
## ***** (V6) ves1proc *****
## Before Matching After Matching
## mean treatment........ 1.4435 1.4456
## mean control.......... 1.2088 1.5109
## std mean diff......... 34.534 -9.6339
##
## mean raw eQQ diff..... 0.24242 0.065312
## med raw eQQ diff..... 0 0
## max raw eQQ diff..... 1 1
##
## mean eCDF diff........ 0.048684 0.013062
## med eCDF diff........ 0.014024 0.0043541
## max eCDF diff........ 0.1805 0.049347
##
## var ratio (Tr/Co)..... 2.0392 0.93506
## T-test p-value........ 9.0068e-10 0.0050611
## KS Bootstrap p-value.. < 2.22e-16 0.062
## KS Naive p-value...... 2.6627e-06 0.37114
## KS Statistic.......... 0.1805 0.049347
##
##
## Before Matching Minimum p.value: < 2.22e-16
## Variable Name(s): ves1proc Number(s): 6
##
## After Matching Minimum p.value: 0.0050611
## Variable Name(s): ves1proc Number(s): 6lindnerMatched %>%
dplyr::select(stent,female,diabetic,acutemi,ejecfrac,ves1proc,abcix) %>%
tbl_summary(by=abcix) %>%
add_overall() %>%
add_p()| Characteristic | Overall, N = 1,3781 | 0, N = 6891 | 1, N = 6891 | p-value2 |
|---|---|---|---|---|
| stent | 938 (68%) | 452 (66%) | 486 (71%) | 0.049 |
| female | 440 (32%) | 210 (30%) | 230 (33%) | 0.2 |
| diabetic | 292 (21%) | 150 (22%) | 142 (21%) | 0.6 |
| acutemi | 230 (17%) | 112 (16%) | 118 (17%) | 0.7 |
| ejecfrac | 55 (45, 56) | 55 (45, 60) | 55 (45, 56) | 0.3 |
| ves1proc | 0.3 | |||
| 0 | 2 (0.1%) | 0 (0%) | 2 (0.3%) | |
| 1 | 842 (61%) | 405 (59%) | 437 (63%) | |
| 2 | 434 (31%) | 231 (34%) | 203 (29%) | |
| 3 | 73 (5.3%) | 38 (5.5%) | 35 (5.1%) | |
| 4 | 27 (2.0%) | 15 (2.2%) | 12 (1.7%) | |
| 1 n (%); Median (IQR) | ||||
| 2 Pearson’s Chi-squared test; Wilcoxon rank sum test; Fisher’s exact test | ||||
Note that the number of stent has not been well balanced after the matching procedure.
For this reason, we use this covariate in the regression model for the outcome.
Now we estimate the causal effect:
fit.out4<-glm(sixMonthDeath~abcix+stent, data=lindnerMatched,family=binomial)
summary(fit.out4)##
## Call:
## glm(formula = sixMonthDeath ~ abcix + stent, family = binomial,
## data = lindnerMatched)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4039 -0.4039 -0.1888 -0.1888 3.1524
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.4089 0.3429 -9.941 < 2e-16 ***
## abcix -1.5529 0.3562 -4.359 1.3e-05 ***
## stent 0.9438 0.3727 2.532 0.0113 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 449.29 on 1377 degrees of freedom
## Residual deviance: 418.68 on 1375 degrees of freedom
## AIC: 424.68
##
## Number of Fisher Scoring iterations: 7We can see that the number of stent is statistically significant in the model,so it has been a good idea to control for it, since it was not well balanced in the matching procedure. As we already have discussed, sometimes also covariates that are not confounders for the effect of the treatment on the outcome could be included in the final model in order to obtain more accurate estimates of the effect.
Fourth option: IPTW
Definition of the weights:
# Definition of weights for ATE
lindner.overlap %<>%mutate(w_ATE=case_when(abcix==1~1/ps,
abcix==0~1/(1-ps)))
#Definition of weights for ATT
lindner.overlap %<>%mutate(w_ATT=case_when(abcix==1~1,
abcix==0~ps/(1-ps))) Check the extreme weights: sometimes it is useful to use truncated or stabilized weights, in order to reduce the variance of the final estimates, but we do not cover here this aspect.
#Check extremes
quantile(lindner.overlap$w_ATE[lindner.overlap$abcix==1],c(0,0.01,0.05,0.95,0.99,1))## 0% 1% 5% 95% 99% 100%
## 1.041578 1.052785 1.076734 1.882227 2.339352 3.627077quantile(lindner.overlap$w_ATE[lindner.overlap$abcix==0],c(0,0.01,0.05,0.95,0.99,1))## 0% 1% 5% 95% 99% 100%
## 1.566665 1.658869 1.791224 6.137759 14.832595 25.687628quantile(lindner.overlap$w_ATT[lindner.overlap$abcix==1],c(0,0.01,0.05,0.95,0.99,1))## 0% 1% 5% 95% 99% 100%
## 1 1 1 1 1 1quantile(lindner.overlap$w_ATT[lindner.overlap$abcix==0],c(0,0.01,0.05,0.95,0.99,1))## 0% 1% 5% 95% 99% 100%
## 0.5666647 0.6588687 0.7912238 5.1377591 13.8325951 24.6876276# Balance diagnostics
#ATE
bal_IPTW_ATE <- dx.wts(x=lindner.overlap$w_ATE,
data=lindner.overlap,
vars=colnames(lindner.overlap)[4:10],
treat.var = colnames(lindner.overlap)[3],
estimand = "ATE")
#ATT
bal_IPTW_ATT <- dx.wts(x=lindner.overlap$w_ATT,
data=lindner.overlap,
x.as.weights = T,
vars=colnames(lindner.overlap)[4:10],
treat.var = colnames(lindner.overlap)[3],
estimand = "ATT")
bal.table(bal_IPTW_ATE) ## $unw
## tx.mn tx.sd ct.mn ct.sd std.eff.sz stat p ks ks.pval
## stent 0.704 0.457 0.586 0.493 0.252 3.541 0.000 0.118 0.005
## height 171.419 10.728 171.458 10.605 -0.004 -0.053 0.958 0.024 1.000
## female 0.333 0.472 0.384 0.487 -0.106 -1.509 0.132 0.051 0.641
## diabetic 0.207 0.406 0.266 0.443 -0.141 -1.962 0.050 0.059 0.450
## acutemi 0.171 0.377 0.061 0.239 0.320 5.537 0.000 0.110 0.012
## ejecfrac 50.559 10.179 52.279 10.313 -0.168 -2.415 0.016 0.108 0.015
## ves1proc 1.443 0.680 1.209 0.476 0.370 6.207 0.000 0.181 0.000
##
## [[2]]
## tx.mn tx.sd ct.mn ct.sd std.eff.sz stat p ks ks.pval
## stent 0.669 0.471 0.668 0.472 0.000 0.005 0.996 0.000 1.000
## height 171.123 10.734 172.398 10.965 -0.119 -1.351 0.177 0.064 0.532
## female 0.347 0.476 0.330 0.471 0.035 0.468 0.640 0.017 1.000
## diabetic 0.225 0.418 0.242 0.429 -0.040 -0.461 0.645 0.017 1.000
## acutemi 0.137 0.344 0.149 0.357 -0.035 -0.307 0.759 0.012 1.000
## ejecfrac 51.077 9.910 51.093 10.182 -0.002 -0.020 0.984 0.044 0.909
## ves1proc 1.368 0.641 1.427 0.690 -0.093 -0.783 0.434 0.032 0.996bal.table(bal_IPTW_ATT) ## $unw
## tx.mn tx.sd ct.mn ct.sd std.eff.sz stat p ks ks.pval
## stent 0.704 0.457 0.586 0.493 0.259 3.541 0.000 0.118 0.005
## height 171.419 10.728 171.458 10.605 -0.004 -0.053 0.958 0.024 1.000
## female 0.333 0.472 0.384 0.487 -0.107 -1.509 0.132 0.051 0.641
## diabetic 0.207 0.406 0.266 0.443 -0.145 -1.962 0.050 0.059 0.450
## acutemi 0.171 0.377 0.061 0.239 0.293 5.537 0.000 0.110 0.012
## ejecfrac 50.559 10.179 52.279 10.313 -0.169 -2.415 0.016 0.108 0.015
## ves1proc 1.443 0.680 1.209 0.476 0.345 6.207 0.000 0.181 0.000
##
## [[2]]
## tx.mn tx.sd ct.mn ct.sd std.eff.sz stat p ks ks.pval
## stent 0.704 0.457 0.703 0.458 0.003 0.034 0.973 0.001 1.000
## height 171.419 10.728 172.793 11.090 -0.128 -1.264 0.206 0.068 0.597
## female 0.333 0.472 0.307 0.462 0.055 0.683 0.495 0.026 1.000
## diabetic 0.207 0.406 0.232 0.423 -0.061 -0.598 0.550 0.025 1.000
## acutemi 0.171 0.377 0.186 0.390 -0.041 -0.313 0.754 0.015 1.000
## ejecfrac 50.559 10.179 50.594 10.085 -0.003 -0.041 0.968 0.035 0.997
## ves1proc 1.443 0.680 1.519 0.743 -0.111 -0.827 0.409 0.043 0.966Finally, let’s estimate the ATE causal effect on the weighted dataset:
# Estimate ATE
design.lindnerATE <- svydesign(ids=~1,
weights = ~w_ATE,
data=lindner.overlap)
fit_itpw_ATE <- svyglm(sixMonthDeath~abcix,
family=binomial,
design=design.lindnerATE)
tbl_regression(fit_itpw_ATE,exponentiate = T)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| abcix | 0.17 | 0.06, 0.49 | <0.001 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
And the ATT:
design.lindnerATT <- svydesign(ids=~1,
weights = ~w_ATT,
data=lindner.overlap)
fit_iptw_ATT <- svyglm(sixMonthDeath~abcix,
family=binomial,
design=design.lindnerATT)## Warning in eval(family$initialize): #successi non interi in un binomial glm!tbl_regression(fit_iptw_ATT,exponentiate = T)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| abcix | 0.15 | 0.05, 0.46 | <0.001 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
Also here it is possible to estimate the 95% CI using boostrap methods (that are in general more robust).
Boostrap confidence intervals for METHOD 1: covariate adjustement
results <- data.frame(ATE=rep(NA,4),ATT=rep(NA,4))f_PSadj <- function(data, indices,outcome.formula) {
d <- data[indices,] # allows boot to select sample
# estimation of ps
m1<-glm(abcix~
acutemi+
ejecfrac+
ves1proc+
stent+
diabetic+
female,
data=d,
family = binomial)
d$ps<-fitted.values(m1)
# overlap
d %<>% mutate(overlap=ifelse(ps>=min(ps[abcix==1]) &
ps<=max(ps[abcix==0]),1,0))
# outcome model
m2<-glm(outcome.formula,data=d,family="binomial",subset = overlap==1)
if(!m2$converged) print("Model did not converged")
d$predY0<-m2$family$linkinv(coef(m2)[1]+coef(m2)[3]*d$ps)
d$predY1<-m2$family$linkinv(coef(m2)[1]+coef(m2)[2]+coef(m2)[3]*d$ps)
Y1<-mean(d$predY1)
Y0<-mean(d$predY0)
Y1_1<-mean(d$predY1[lindner$abcix==1])
Y0_1<-mean(d$predY0[lindner$abcix==1])
ATE_PSadj<-(Y1/(1-Y1))/(Y0/(1-Y0))
ATT_PSadj<-(Y1_1/(1-Y1_1))/(Y0_1/(1-Y0_1))
return(c(ATE_PSadj,ATT_PSadj))
}
res_boot <- function(obj.boot,type="percent",digits=3){
suppressWarnings({
orig <- round(obj.boot$t0,digits)
ciATE <- paste(round(boot.ci(obj.boot,index=1)[[type]][4:5],digits),collapse = "-")
ciATT <- paste(round(boot.ci(obj.boot,index=2)[[type]][4:5],digits),collapse = "-")
})
res <- paste(orig,rbind(ciATE,ciATT),sep="(")
res.u <- paste(res,rep(")",2),sep = "")
return(res.u)
}boot.out <- boot(data=lindner, statistic=f_PSadj,R=1000,outcome.formula=fit.out$formula)
print(boot.out)##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = lindner, statistic = f_PSadj, R = 1000, outcome.formula = fit.out$formula)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 0.2363187 0.02393070 0.1318934
## t2* 0.2363050 0.02394347 0.1318932# Get 95% confidence interval
results$ATE[1]<-res_boot(boot.out,digits=4)[1]
results$ATT[1]<-res_boot(boot.out,digits=4)[2]
results## ATE ATT
## 1 0.2363(0.0929-0.5899) 0.2363(0.0929-0.5899)
## 2 <NA> <NA>
## 3 <NA> <NA>
## 4 <NA> <NA>Boostrap confidence intervals for METHOD 2: stratification
Very often with the stratification method there are many problems of convergence of the regression algorithm, since in the strata we have very few events !
f_PSstrat <- function(data, indices) {
d <- data[indices,] # allows boot to select sample
m1<-glm(abcix~
acutemi+
ejecfrac+
ves1proc+
stent+
diabetic+
female,
data=d,family = binomial)
d$ps<-fitted.values(m1)
quart_PS<-quantile(d$ps,c(0,0.33,0.66,1))
d$strata<-cut(d$ps, quart_PS, labels=c(1:3))
for (i in 1:3){
ms<-glm(sixMonthDeath~abcix,data=d,subset = strata==i,family="binomial")
beta.treat[i]<-coef(ms)[2]
}
ATE <- exp(sum(beta.treat*nstrata)/nrow(lindner))
ATT <- exp(sum(beta.treat*treated.strata)/sum(treated.strata))
return(c(ATE,ATT))
}boot.out4 <- boot(data=lindner, statistic=f_PSstrat,R=1000)
# Get 95% confidence interval
results$ATE[2]<-res_boot(boot.out4,digits=3)[1]
results$ATT[2]<-res_boot(boot.out4,digits=3)[2]Robust confidence intervals for METHOD 3: matching
We now estimate the robust standard errors related to the estimate on the matched dataset:
cov <- vcovHC(fit.out4, type = "HC0")
std.err <- sqrt(diag(cov))
q.val <- qnorm(0.975)
r <- cbind(
Estimate = coef(fit.out4)
, "Robust SE" = std.err
, z = (coef(fit.out4)/std.err)
, "Pr(>|z|) "= 2 * pnorm(abs(coef(fit.out4)/std.err), lower.tail = FALSE)
, LL = coef(fit.out4) - q.val * std.err
, UL = coef(fit.out4) + q.val * std.err
)
#Exponential to get the OR
results$ATT[3]<- paste0(round(exp(r[2,1]),4),"(",round(exp(r[2,5]),4),"-",round(exp(r[2,6]),4),")")Boostrap confidence intervals for method 4: IPTW
f_IPTW <- function(data, indices) {
d <- data[indices,] # allows boot to select sample
# estimation of ps
m1<-glm(abcix~
acutemi+
ejecfrac+
ves1proc+
stent+
diabetic+
female,
data=d,
family = binomial)
d$ps<-fitted.values(m1)
# overlap
d %<>% mutate(overlap=ifelse(ps>=min(ps[abcix==1]) &
ps<=max(ps[abcix==0]),1,0)) %>%
filter(overlap==1)
# Definition of weights for ATE
d %<>%mutate(w_ATE=case_when(abcix==1~1/ps,
abcix==0~1/(1-ps)))
#Definition of weights for ATT
d %<>%mutate(w_ATT=case_when(abcix==1~1,
abcix==0~ps/(1-ps)))
# Estimate ATE
design.lindnerATE <- svydesign(ids=~1,
weights = ~w_ATE,
data=d)
fit_itpw_ATE <- svyglm(sixMonthDeath~abcix,
family=binomial,
design=design.lindnerATE)
# Estimate ATT
design.lindnerATT <- svydesign(ids=~1,
weights = ~w_ATT,
data=d)
fit_iptw_ATT <- svyglm(sixMonthDeath~abcix,
family=binomial,
design=design.lindnerATT)
ATE <- exp(fit_itpw_ATE$coefficients[2])
ATT <- exp(fit_iptw_ATT$coefficients[2])
return(c(ATE,ATT))
}boot.out5 <- boot(data=lindner, statistic=f_IPTW,R=1000)
results$ATE[4]<-res_boot(boot.out5,digits=4)[1]
results$ATT[4]<-res_boot(boot.out5,digits=4)[2]Final Comparison across methods !
# recall the crude OR from the original dataset
tbl_regression(fit.crude,exponentiate=T)| Characteristic | OR1 | 95% CI1 | p-value |
|---|---|---|---|
| abcix | 0.30 | 0.13, 0.66 | 0.003 |
| 1 OR = Odds Ratio, CI = Confidence Interval | |||
# PS results
row.names(results) <- c("PS Adj Linear","Stratification","Matching","IPTW")
results## ATE ATT
## PS Adj Linear 0.2363(0.0929-0.5899) 0.2363(0.0929-0.5899)
## Stratification 0.196(0-43.913) 0.202(0.001-35.486)
## Matching <NA> 0.2116(0.1046-0.4283)
## IPTW 0.1749(0.0592-0.6322) 0.1525(0.0474-0.6406)We can observe that using the propensity score based methods we obtain a stronger estimate of the treatment effect with respect to the crude odds ratio. We can also observe that stratification produce very unstable results due to the low sample size in the different strata. In conclusion, adjusting for confounders was important in this context !
Appendix 1 : METHOD 2 “Visually” checking balance across strata
# Function for graphical diagnostic to check balance
balance.Strat <- function(x){
if(length(unique(lindner[,x]))<10){
cat.psa(categorical = lindner[,x],
treatment = lindner$abcix,
strata = lindner$strata,
ylab="Relative Frequency")
title(x)
}
else{
box.psa(continuous = lindner[,x],
treatment = lindner$abcix,
strata = lindner$strata)
title(x)
}
}var <- lindner %>%dplyr::select(stent,
female,
diabetic,
acutemi,
ejecfrac,
ves1proc) %>% colnames() balance.Strat(var[1])
balance.Strat(var[2])
balance.Strat(var[3])
balance.Strat(var[4])
balance.Strat(var[5])## Warning in xy.coords(x, y): NA introdotti per coercizione
## Warning in xy.coords(x, y): NA introdotti per coercizione
balance.Strat(var[6])
Appendix 2: Methods to estimate weights
There are many R libraries that can be used to estimate weights see for example: https://cran.r-project.org/web/packages/PSweight/index.html This library enables the estimation and inference of average causal effects with binary and multiple treatments using overlap weights (ATO), inverse probability of treatment weights (ATE), average treatment effect among the treated weights (ATT), matching weights (ATM) and entropy weights (ATEN), with and without propensity score trimming. Reference: https://arxiv.org/pdf/2010.08893