1. Preparation for Analysis
1.1 Platform Info
print(version)platform x86_64-pc-linux-gnu
arch x86_64
os linux-gnu
system x86_64, linux-gnu
status
major 4
minor 1.3
year 2022
month 03
day 10
svn rev 81868
language R
version.string R version 4.1.3 (2022-03-10)
nickname One Push-Up 1.2 Loading packages
library(data.table)
library(dplyr)
library(tidyr)
library(ggplot2)
library(fixest)
library(table1)
library(ggsci)
library(cowplot)
library(kableExtra)2. Process of data manipulation
2.1 Import data
xiyi_2020 <- fread("~/EMR from 2018 to 2020/Dataset/Tidy Data/Xiyi_2020_tidy_2022-04-28.csv")xiyi_2019 <- fread("~/EMR from 2018 to 2020/Dataset/Tidy Data/Xiyi_2019_tidy_2022-04-28.csv")xiyi_2018 <- fread("~/EMR from 2018 to 2020/Dataset/Tidy Data/Xiyi_2018_tidy_2022-04-28.csv")2.2 Sample selection
#---Hospital selection---#
keep_hos_shaanxi <- c("陕西省人民医院",
"西安交通大学医学院第一附属医院",
"西安交通大学医学院第二附属医院",
"延安大学附属医院",
"榆林市第一医院",
"安康市中心医院",
"宝鸡市中心医院",
"汉中市中心医院",
"商洛市中心医院",
"西安国际医学中心医院",
"咸阳市中心医院",
"延安大学咸阳医院",
"榆林市第二医院",
"咸阳市第一人民医院",
"三二〇一医院",
"西安大兴医院",
"西安市红会医院",
"陕西省核工业二一五医院",
"铜川市人民医院",
"渭南市中心医院",
"延安市人民医院",
"宝鸡高新医院",
"长安医院",
"西安市中心医院",
"陕西省康复医院",
"西安市第三医院",
"宝鸡市人民医院",
"富平县医院",
"宝鸡市陈仓医院",
"西乡县人民医院",
"西安市第四医院",
"西安凤城医院",
"城固县医院",
"宝鸡市第三医院",
"兵器工业五二一医院",
"西安高新医院",
"西安医学院第一附属医院",
"蒲城县医院",
"铜川矿务局中心医院",
"西安市长安区医院",
"渭南市第二医院",
"西安医学院第二附属医院"
)
#---Variable selection---#
keep_vars <-
c("insti_name", "insti_code", "cardnumber",
"inpatient_time", "record_num", "dob_std", "age",
"year_std", "season", "sex_std", "payment_std", "nationality_std",
"ethic_std", "marriage_std", "admission_way_std",
"age_infant", "infant_weight_at_born", "infant_weight_admission", "place_of_born",
"hometown", "occupation", "admission_date_std",
"discharge_date_std", "indays", "total_fee", "oop",
"lab_diag_fee", "radio_diag_fee", "west_drug_fee",
"antibat_drug_fee", "chn_drug_fee",
"disposable_treat_mat_fee", "disposable_surg_mat_fee",
"diag_inpatient_code", "surg_code")
xiyi_2020_analysis <- xiyi_2020[insti_name %in% keep_hos_shaanxi, ..keep_vars]
xiyi_2019_analysis <- xiyi_2019[insti_name %in% keep_hos_shaanxi, ..keep_vars]
xiyi_2018_analysis <- xiyi_2018[insti_name %in% keep_hos_shaanxi, ..keep_vars]
rm(xiyi_2020, xiyi_2019, xiyi_2018)2.3 Data clean
clean_obs <- function(df) {
total_fee_q1 <- quantile(df[, "total_fee"], probs = 0.01, na.rm = TRUE)[[1]]
total_fee_q99 <- quantile(df[, "total_fee"], probs = 0.99, na.rm = TRUE)[[1]]
df_analysis <- df[(sex_std %in% c(1, 2) &
nationality_std == "CHN" &
!is.na(admission_way_std) &
(discharge_date_std - admission_date_std + 1) <= 90 &
total_fee >= total_fee_q1 &
total_fee <= total_fee_q99 &
lab_diag_fee >= 0 &
radio_diag_fee >= 0 &
west_drug_fee >= 0 &
disposable_treat_mat_fee + disposable_surg_mat_fee >= 0
)]
}
xiyi_2020_sub <- clean_obs(xiyi_2020_analysis)
rm(xiyi_2020_analysis)
xiyi_2019_sub <- clean_obs(xiyi_2019_analysis)
rm(xiyi_2019_analysis)
xiyi_2018_sub <- clean_obs(xiyi_2018_analysis)
rm(xiyi_2018_analysis)2.4 Derived variables
#---
xiyi_2020_sub[, c("len_admission",
"fee_insure",
"surg_code_std",
"mat_fee",
"diag_inpatient_code_1",
"diag_inpatient_code_3") := list(discharge_date_std - admission_date_std + 1,
total_fee - oop,
ifelse(surg_code == "", 0, 1),
disposable_treat_mat_fee + disposable_surg_mat_fee,
toupper(substr(diag_inpatient_code, 1, 1)),
toupper(substr(diag_inpatient_code, 1, 3))
)]
#---
xiyi_2019_sub[, c("len_admission",
"fee_insure",
"surg_code_std",
"mat_fee",
"diag_inpatient_code_1",
"diag_inpatient_code_3") := list(discharge_date_std - admission_date_std + 1,
total_fee - oop,
ifelse(surg_code == "", 0, 1),
disposable_treat_mat_fee + disposable_surg_mat_fee,
toupper(substr(diag_inpatient_code, 1, 1)),
toupper(substr(diag_inpatient_code, 1, 3))
)]
#---
xiyi_2018_sub[, c("len_admission",
"fee_insure",
"surg_code_std",
"mat_fee",
"diag_inpatient_code_1",
"diag_inpatient_code_3") := list(discharge_date_std - admission_date_std + 1,
total_fee - oop,
ifelse(surg_code == "", 0, 1),
disposable_treat_mat_fee + disposable_surg_mat_fee,
toupper(substr(diag_inpatient_code, 1, 1)),
toupper(substr(diag_inpatient_code, 1, 3))
)]2.5 Sample description
Table. 1
df_tab1 <-
rbind(xiyi_2018_sub, xiyi_2019_sub, xiyi_2020_sub) %>%
.[, "yday" := yday(admission_date_std)] %>%
filter(., yday < 23 | (yday > 50 & yday <= 335))
#---#
lbs_vars <- c("age", "sex_std", "payment_std",
"admission_way_std", "len_admission", "surg_code_std",
"total_fee", "oop", "fee_insure", "lab_diag_fee",
"radio_diag_fee", "west_drug_fee", "mat_fee")
df_tab1$sex_std <- as.character(df_tab1$sex_std)
df_tab1$payment_std <- as.character(df_tab1$payment_std)
df_tab1$admission_way_std <- as.character(df_tab1$admission_way_std)
df_tab1$surg_code_std <- as.character(df_tab1$surg_code_std)
lbs <- c("Age (yr)", "Sex", "Insurance Payment",
"Admission Way", "Length of Admission", "Surgical",
"Total Hospitalization Expenses (CNY)",
"Out-of-Pocket (CNY)",
"Claim Expenses (CNY)",
"Laboratory Test Expenses (CNY)",
"Imaging Test Expenses (CNY)",
"Drug Expenses (CNY)",
"Medical Consumables Expenses (CNY)")
for (i in 1:length(lbs)) {
label(df_tab1[[lbs_vars[i]]]) <- lbs[i]
}
#---#
table1(~ age + sex_std + payment_std +
admission_way_std + surg_code_std +
len_admission + total_fee + oop + fee_insure +
lab_diag_fee + radio_diag_fee + west_drug_fee +
mat_fee | factor(year_std),
data = df_tab1,
footnote = "CNY: Chinese Yuan",
caption = "Table.1 Summary of sample")
rm(df_tab1)3. Results
3.1 Descriptive analysis
3.1.1 Comparing of number of hospitalized admission of 36 tertiary hospitals from 2018 to 2020
#---Admission number
admission_daliy_2020_num <- xiyi_2020_sub[, .N, by = c("year_std", "admission_date_std")]
admission_daliy_2019_num <- xiyi_2019_sub[, .N, by = c("year_std", "admission_date_std")]
admission_daliy_2018_num <- xiyi_2018_sub[, .N, by = c("year_std", "admission_date_std")]
admission_daliy_num <- rbind(admission_daliy_2018_num,
admission_daliy_2019_num,
admission_daliy_2020_num)
admission_daliy_num[, "adm_date_day" := yday(admission_date_std)]Figure 1A
p_admission <-
ggplot(filter(admission_daliy_num, adm_date_day <= 92)) +
geom_line(aes(x = adm_date_day, y = N, color = factor(year_std))) +
geom_smooth(aes(x = adm_date_day, y = N, color = factor(year_std)),
method = "loess") +
labs(x = "Day of year", y = "No. of Admission") +
scale_colour_jama(name = "Year") +
scale_x_continuous(breaks = seq(1, 92, 7)) +
geom_vline(xintercept = c(23, 50), linetype="dotdash") +
geom_rect(aes(xmin = 23, xmax = 50, ymin=-Inf, ymax=Inf), fill = "yellow", alpha = .005) +
geom_text(aes(x = 23, y = 6500, label = "2020-1-23 \n First case"),
size = 5, family = "Times") +
geom_text(aes(x = 50, y = 6500, label = "2020-2-19 \n Last case"),
size = 5, family = "Times") +
theme_bw() +
theme(axis.text.x = element_text(vjust = 1, size = 15, family = "Times"),
text = element_text(family = "Times", size = 15),
legend.position = "bottom",
plot.title = element_text(hjust = 0.5, size = 20, family = "Times"))3.1.3 Percentage of patients surgery among 2018, 2019 and 2020
#---Admission number
admission_daliy_2020_surg <- xiyi_2020_sub[, .N,
by = c("year_std",
"admission_date_std",
"surg_code_std")]
admission_daliy_2019_surg <- xiyi_2019_sub[, .N,
by = c("year_std",
"admission_date_std",
"surg_code_std")]
admission_daliy_2018_surg <- xiyi_2018_sub[, .N,
by = c("year_std",
"admission_date_std",
"surg_code_std")]
admission_daliy_surg <- rbind(admission_daliy_2018_surg,
admission_daliy_2019_surg,
admission_daliy_2020_surg)
admission_daliy_surg_join <- left_join(admission_daliy_surg[surg_code_std == 0],
admission_daliy_surg[surg_code_std == 1],
by = c("year_std",
"admission_date_std")) %>%
mutate("N" = N.y * 100 / (N.x + N.y))
admission_daliy_surg_join[, "adm_date_day" := yday(admission_date_std)]
rm(admission_daliy_2018_surg, admission_daliy_2019_surg, admission_daliy_2020_surg)Figure S1 in supplementary
p_adm_surg <-
ggplot(filter(admission_daliy_surg_join, adm_date_day < 335)) +
geom_line(aes(x = adm_date_day, y = N, color = factor(year_std))) +
geom_smooth(aes(x = adm_date_day, y = N, color = factor(year_std)),
method = "loess") +
labs(x = "Day of year", y = "% of surgery") +
scale_colour_jama(name = "Year") +
scale_x_continuous(breaks = seq(1, 337, 14)) +
geom_vline(xintercept = c(23, 50), linetype="dotdash") +
geom_rect(aes(xmin = 23, xmax = 50, ymin=-Inf, ymax=Inf), fill = "yellow", alpha = .005) +
geom_text(aes(x = 23, y = 55, label = "2020-1-23 \n First case"),
size = 5, family = "Times") +
geom_text(aes(x = 50, y = 58, label = "2020-2-19 \n Last case"),
size = 5, family = "Times") +
theme_bw() +
theme(axis.text.x = element_text(vjust = 1, size = 15, family = "Times"),
text = element_text(family = "Times", size = 15),
legend.position = "bottom",
plot.title = element_text(hjust = 0.5, size = 20, family = "Times"))3.1.4 Percentage of admission from emergency room visit among 2018, 2019 and 2020
admission_daliy_2020_em <- xiyi_2020_sub[, .N,
by = c("year_std",
"admission_date_std",
"admission_way_std")]
admission_daliy_2019_em <- xiyi_2019_sub[, .N,
by = c("year_std",
"admission_date_std",
"admission_way_std")]
admission_daliy_2018_em <- xiyi_2018_sub[, .N,
by = c("year_std",
"admission_date_std",
"admission_way_std")]
#---------#
admission_daliy_2020_total <- xiyi_2020_sub[, .N,
by = c("year_std",
"admission_date_std")]
admission_daliy_2019_total <- xiyi_2019_sub[, .N,
by = c("year_std",
"admission_date_std")]
admission_daliy_2018_total <- xiyi_2018_sub[, .N,
by = c("year_std",
"admission_date_std")]
#---------#
em_percent <- function(dt1, dt2) {
df <- left_join(dt1[admission_way_std == 1],
dt2,
by = c("year_std", "admission_date_std")) %>%
mutate("N" = N.x * 100 / N.y)
return(df)
}
adm_2020_em <- em_percent(admission_daliy_2020_em, admission_daliy_2020_total)
adm_2019_em <- em_percent(admission_daliy_2019_em, admission_daliy_2019_total)
adm_2018_em <- em_percent(admission_daliy_2018_em, admission_daliy_2018_total)
admission_daily_em <- rbind(adm_2020_em,
adm_2019_em,
adm_2018_em)
admission_daily_em[, "adm_date_day" := yday(admission_date_std)]
rm(admission_daliy_2020_em, admission_daliy_2020_total,
admission_daliy_2019_em, admission_daliy_2019_total,
admission_daliy_2018_em, admission_daliy_2018_total)Figure S3 in supplementary
p_adm_em <-
ggplot(filter(admission_daily_em, adm_date_day < 335)) +
geom_line(aes(x = adm_date_day, y = N, color = factor(year_std))) +
geom_smooth(aes(x = adm_date_day, y = N, color = factor(year_std)),
method = "loess") +
labs(x = "Day of year", y = "% of admission from emergency room visit") +
scale_colour_jama(name = "Year") +
scale_x_continuous(breaks = seq(1, 337, 14)) +
geom_vline(xintercept = c(23, 50), linetype="dotdash") +
geom_rect(aes(xmin = 23, xmax = 50, ymin=-Inf, ymax=Inf), fill = "yellow", alpha = .005) +
geom_text(aes(x = 23, y = 55, label = "2020-1-23 \n First case"),
size = 5, family = "Times") +
geom_text(aes(x = 50, y = 58, label = "2020-2-19 \n Last case"),
size = 5, family = "Times") +
theme_bw() +
theme(axis.text.x = element_text(vjust = 1, size = 15, family = "Times"),
text = element_text(family = "Times", size = 15),
legend.position = "bottom",
plot.title = element_text(hjust = 0.5, size = 20, family = "Times"))3.1.5 Comparing the tread of total expenses, out-of-pocket, and claiming expenses of hospitalized patients in 36 tertiary hospitals from 2018 to 2020
#---Expenditure per admission day
admission_daliy_2020_fee <- xiyi_2020_sub[, .(fee_mean = mean(total_fee),
oop_mean = mean(oop),
insur_mean = mean(fee_insure)
),
by = c("year_std",
"admission_date_std")]
admission_daliy_2019_fee <- xiyi_2019_sub[, .(fee_mean = mean(total_fee),
oop_mean = mean(oop),
insur_mean = mean(fee_insure)
),
by = c("year_std",
"admission_date_std")]
admission_daliy_2018_fee <- xiyi_2018_sub[, .(fee_mean = mean(total_fee),
oop_mean = mean(oop),
insur_mean = mean(fee_insure)
),
by = c("year_std",
"admission_date_std")]
admission_daliy_fee <- rbind(admission_daliy_2018_fee,
admission_daliy_2019_fee,
admission_daliy_2020_fee)
admission_daliy_fee[, "adm_date_day" := yday(admission_date_std)]
#---Plot
line_plot <- function(df, var, ytitle) {
p <-
filter(df, adm_date_day <= 335) %>%
rename("outcome" = all_of(var)) %>%
ggplot() +
geom_line(aes(x = adm_date_day, y = outcome, color = factor(year_std))) +
geom_smooth(aes(x = adm_date_day, y = outcome, color = factor(year_std))) +
labs(x = "Day of year", y = ytitle) +
scale_color_jama(name = "") +
scale_x_continuous(breaks = seq(1, 337, 14)) +
scale_y_continuous(labels = function(x) ifelse(x >= 10000,
paste0(x / 1000, ",", "000"), x)) +
geom_vline(xintercept = c(23, 50), linetype="dotdash", color = c("red", "black")) +
geom_rect(aes(xmin = 23, xmax = 50, ymin=-Inf, ymax=Inf),
fill = "yellow", alpha = .005) +
theme_bw() +
theme(axis.text.x = element_text(size = 12, family = "Times"),
text = element_text(family = "Times", size = 12),
legend.position = c(0.92, 0.08),
legend.direction = "horizontal",
legend.text = element_text(size = 6),
legend.background = element_blank(),
plot.title = element_text(hjust = 0.5, size = 15, family = "Times"))
return(p)
}1) Total Expense of Each Patient
Figure 3A
p_total <-
line_plot(admission_daliy_fee, "fee_mean",
"Total Expense (CNY)")2) Out of Pocket of Each Patient
Figure 3B
p_oop <-
line_plot(admission_daliy_fee, "oop_mean",
"Out of Pocekt (CNY)")3) Expense paid by insurance of Each Patient
Figure 3C
p_insur <-
line_plot(admission_daliy_fee, "insur_mean",
"Claiming Expenses (CNY)")4) Merge plots
Figure 3: Merge Figure 3A, 3B and 3C
plot_grid(p_total, p_oop, p_insur, nrow = 3,
label_fontfamily = "Times",
labels = c("A", "B", "C"))
# "(A) Total Expenses", "(B) Out-of-Pocket", "(C) Claiming Expenses"3.2 Causal inference
3.2.1 Prepare data for causal inference
xiyi_2020_sub[, c("yday", "week") := list(yday(admission_date_std),
week(admission_date_std))]
xiyi_2019_sub[, c("yday", "week") := list(yday(admission_date_std),
week(admission_date_std))]
#---以2月19日为最后一天
xiyi_causal <- rbind(xiyi_2019_sub[yday < 23 | (yday > 50 & yday <= 335)],
xiyi_2020_sub[yday < 23 | (yday > 50 & yday <= 335)])
xiyi_causal[, c("period",
"treat",
"day_to_normal",
"week_to_normal",
"day_treated",
"week_treated",
"yday_for_sunab",
"week_for_sunab") := list(ifelse(yday >= 51, 1, 0),
ifelse(year_std == 2020, 1, 0),
ifelse(yday >= 51, yday - 51, yday - 23),
ifelse(week >= 8, week - 8, week - 4),
51 * ifelse(year_std == 2020, 1, 0),
8 * ifelse(year_std == 2020, 1, 0),
ifelse(yday < 23, yday + 28, yday),
ifelse(week < 8, week + 4, week)
)]3.2.2 Main results of DID
Sun and Abraham (2020)
models_sunab <- function(df) {
#---Total Expenses---#
#---Week
mod_fee_week <- feols(log(total_fee + 1) ~ sunab(week_treated, week_for_sunab) +
age + I(age^2) + factor(sex_std) + diag_inpatient_code_1 +
factor(payment_std) + factor(admission_way_std) |
week + insti_code,
cluster = "insti_code",
data = df[week != 4, ])
#---Out of Pocket---#
#---Week
mod_oop_week <- feols(log(oop + 1) ~ sunab(week_treated, week_for_sunab) +
age + I(age^2) + factor(sex_std) + diag_inpatient_code_1 +
factor(payment_std) + factor(admission_way_std) |
week + insti_code,
cluster = "insti_code",
data = df[week != 4, ])
#---Claiming Expense---#
#---Week
mod_insur_week <- feols(log(fee_insure + 1) ~ sunab(week_treated, week_for_sunab) +
age + I(age^2) + factor(sex_std) + diag_inpatient_code_1 +
factor(payment_std) + factor(admission_way_std) |
week + insti_code,
cluster = "insti_code",
data = df[week != 4, ])
return(list(mod_fee_week,
mod_oop_week,
mod_insur_week))
}
mods <- models_sunab(xiyi_causal)
mod_fee_week <- mods[[1]]
mod_oop_week <- mods[[2]]
mod_insur_week <- mods[[3]]3.2.3 Coefficient plot
Customized coefplot function for fixest package
fixest_coefplot <- function(mod, treatname, ref, xlabel) {
coef_df <- data.frame(mod$coeftable)
coef_df$varname <- row.names(coef_df)
coef_df_select <- coef_df[stringr::str_detect(coef_df$varname, treatname), ]
coef_df_select <- rbind(coef_df_select,
data.frame("Estimate" = 0, "Std..Error" = 0, "t.value" = 0,
"Pr...t.." = 0, "varname" = paste0("ref::", ref)))
coef_df_select$time <- stringr::str_extract(coef_df_select$varname, "-?\\d+")
ggplot(coef_df_select, aes(x = as.numeric(time), y = Estimate)) +
geom_point(fill = "blue", color = "blue", size = 2, shape = 22) +
geom_errorbar(aes(ymin = Estimate - 1.96*Std..Error,
ymax = Estimate + 1.96*Std..Error), width = 0.5) +
geom_hline(aes(yintercept = 0)) +
geom_vline(aes(xintercept = -1), linetype = "dashed") +
scale_x_continuous(breaks = scales::breaks_extended(n = 10)) +
scale_y_continuous(breaks = scales::breaks_extended(n = 10)) +
labs(y = "Estimate and 95% Conf. Int.", x = xlabel) +
theme_bw() +
theme(text = element_text(size = 15, family = "Times"))
}Figure 4
p_mod_total <- fixest_coefplot(mod_fee_week, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
p_mod_oop <- fixest_coefplot(mod_oop_week, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
p_mod_insur <- fixest_coefplot(mod_insur_week, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
plot_grid(p_mod_total, p_mod_oop, p_mod_insur, nrow = 3, align = "v",
label_y = 0.95,
label_x = 0.35,
label_fontfamily = "Times",
labels = c("(A) Total Expenses", "(B) OOP Expenses", "(C) Claiming Expenses"))3.2.4 Coefficient tables
Customized table function for fixest package
fixest_tab <- function(mod) {
coef_df <- data.frame(mod$coeftable)
coef_df$varname <- row.names(coef_df)
coef_df_select <- coef_df[stringr::str_detect(coef_df$varname, "week_for_sunab"), ]
coef_df_select <- rbind(coef_df_select,
data.frame("Estimate" = 0, "Std..Error" = 0, "t.value" = 0,
"Pr...t.." = 0, "varname" = paste0("ref::", -1)))
coef_df_select <- rename(coef_df_select, "pvalue" = "Pr...t..") %>%
mutate("time" = as.numeric(stringr::str_extract(varname, "-?\\d+")),
"Coef" = as.character(round(Estimate, 3)),
"result" = case_when(pvalue < 0.05 & pvalue > 0.01 ~ paste(Coef, "**"),
pvalue <= 0.01 ~ paste(Coef, "***"),
pvalue >= 0.05 ~ Coef
),
"se" = round(Std..Error, 4),
"lower" = round(Estimate - 1.96 * Std..Error, 3),
"upper" = round(Estimate + 1.96 * Std..Error, 3)
)
coef_df_select <- coef_df_select[order(coef_df_select$time), ]
return(coef_df_select[c("varname", "result", "se", "lower", "upper")])
}Table S1
total_table <- fixest_tab(mod_fee_week)
oop_table <- fixest_tab(mod_oop_week)
insur_table <- fixest_tab(mod_insur_week)
mod_table <- cbind(total_table, oop_table, insur_table)
mod_table[, -c(6, 11)] %>%
kbl(caption = "Regreesion results of models", row.names = FALSE,
col.names = c("Variable", rep(c("Coef", "S.E.", "95% CI Lower", "95% CI Upper"), 3))) %>%
kable_classic(full_width = F, html_font = "Times") %>%
add_header_above(c(" " = 1, "Total Expenses" = 4,
"Out of Pocket" = 4, "Claiming Expenses" = 4))3.3 Placebo Analysis
xiyi_2018_sub[, c("yday", "week") := list(yday(admission_date_std),
week(admission_date_std))]
#---
xiyi_causal_placebo <- rbind(xiyi_2018_sub[yday < 23 | (yday > 50 & yday <= 335)],
xiyi_2019_sub[yday < 23 | (yday > 50 & yday <= 335)])
#---
xiyi_causal_placebo[, c("period",
"treat",
"day_to_normal",
"week_to_normal",
"day_treated",
"week_treated",
"yday_for_sunab",
"week_for_sunab") := list(ifelse(yday >= 51, 1, 0),
ifelse(year_std == 2019, 1, 0),
ifelse(yday >= 51, yday - 51, yday - 23),
ifelse(week >= 8, week - 8, week - 4),
51 * ifelse(year_std == 2019, 1, 0),
8 * ifelse(year_std == 2019, 1, 0),
ifelse(yday < 23, yday + 28, yday),
ifelse(week < 8, week + 4, week)
)]
rm(xiyi_2018_sub, xiyi_2019_sub, xiyi_2020_sub)
#----#
mods_pb <- models_sunab(xiyi_causal_placebo)
mod_fee_week_pb <- mods_pb[[1]]
mod_oop_week_pb <- mods_pb[[2]]
mod_insur_week_pb <- mods_pb[[3]]Figure 5
p_mod_total_pb <- fixest_coefplot(mod_fee_week_pb, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
p_mod_oop_pb <- fixest_coefplot(mod_oop_week_pb, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
p_mod_insur_pb <- fixest_coefplot(mod_insur_week_pb, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
plot_grid(p_mod_total_pb, p_mod_oop_pb, p_mod_insur_pb, nrow = 3, align = "v",
label_y = 0.95,
label_x = 0.4,
label_fontfamily = "Times",
labels = c("(A) Total Expenses", "(B) OOP Expenses", "(C) Claiming Expenses"))3.4 Heterogeneous and sensitivity analysis
3.4.1 Heterogeneous
private_hos <- c("西安大兴医院", "长安医院",
"西安凤城医院", "西安高新医院", "延安大学咸阳医院")
xiyi_causal[, "private_yn" := ifelse(insti_name %in% private_hos, 1, 0)]
mods_priv <- models_sunab(xiyi_causal[private_yn == 1,])
mods_pub <- models_sunab(xiyi_causal[private_yn == 0,])mode_data <- function(mod, treatname, ref, var) {
coef_df <- data.frame(mod$coeftable)
coef_df$varname <- row.names(coef_df)
coef_df_select <- coef_df[stringr::str_detect(coef_df$varname, treatname), ]
coef_df_select <- rbind(coef_df_select,
data.frame("Estimate" = 0, "Std..Error" = 0, "t.value" = 0,
"Pr...t.." = 0, "varname" = paste0("ref::", ref)))
coef_df_select$time <- stringr::str_extract(coef_df_select$varname, "-?\\d+")
coef_df_select$group <- var
return(coef_df_select)
}
priv_pub_df <- rbind(mode_data(mods_priv[[1]], "week_for_sunab", -1, "Private"),
mode_data(mods_pub[[1]], "week_for_sunab", -1, "Public"))Figure 6A
p_priv_pub <-
ggplot(priv_pub_df, aes(x = as.numeric(time), y = Estimate, fill = group)) +
geom_point(position = position_dodge(width = 0.8), size = 2, shape = 22) +
geom_errorbar(aes(ymin = Estimate - 1.96*Std..Error,
ymax = Estimate + 1.96*Std..Error), width = 0.5,
position = position_dodge(width = 0.8)) +
geom_hline(aes(yintercept = 0)) +
geom_vline(aes(xintercept = -1), linetype = "dashed") +
scale_fill_npg(name = "") +
scale_x_continuous(breaks = scales::breaks_extended(n = 10)) +
scale_y_continuous(breaks = scales::breaks_extended(n = 10)) +
labs(y = "Estimate and 95% Conf. Int.",
x = "Time to reopen (Week, 0 = Week 8)") +
theme_bw() +
theme(text = element_text(size = 15, family = "Times"),
legend.position = "bottom")
rm(mods_priv, mods_pub)3.5 Sensitivity
3.5.1 Laboratory test expenses
Figure 6B
#---Week
mod_lab_week <- feols(log(lab_diag_fee + 1) ~ sunab(week_treated, week_for_sunab) +
age + I(age^2) + factor(sex_std) + diag_inpatient_code_1 +
factor(payment_std) + factor(admission_way_std) |
week + insti_code,
cluster = "insti_code",
data = xiyi_causal[(week != 4) & (private_yn == 0), ])
p_lab <-
fixest_coefplot(mod_lab_week, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
rm(mod_lab_week)3.5.2 Medical consumables expenses
Figure 6C
#---Week
mod_mat_week <- feols(log(mat_fee + 1) ~ sunab(week_treated, week_for_sunab) +
age + I(age^2) + factor(sex_std) + diag_inpatient_code_1 +
factor(payment_std) + factor(admission_way_std) |
week + insti_code,
cluster = "insti_code",
data = xiyi_causal[(week != 4) & (private_yn == 0), ])
p_mat <-
fixest_coefplot(mod_mat_week, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
rm(mod_mat_week)3.5.3 Drug expenses
Figure 6D
#---Week
mod_westdrug_week <- feols(log(west_drug_fee + 1) ~ sunab(week_treated, week_for_sunab) +
age + I(age^2) + factor(sex_std) + diag_inpatient_code_1 +
factor(payment_std) + factor(admission_way_std) + surg_code_std |
week + insti_code,
cluster = "insti_code",
data = xiyi_causal[(week != 4) & (private_yn == 0), ])
p_west <-
fixest_coefplot(mod_westdrug_week, "week_for_sunab", -1,
"Time to reopen (Week, 0 = Week 8)")
rm(mod_westdrug_week)Figure 6: Merge Figure 6A, 6B, 6C and 6D
plot_grid(p_priv_pub, p_lab, p_mat, p_west, nrow = 2, align = "v",
label_fontfamily = "Times",
labels = c("A", "B",
"C", "D"))
# "(A) Private & Public", "(B) Laboratory Test",
# "(C) Medical Consumables Expenses", "(D) Drug Expenses"