DPCビッグデータ分析

# 必要なパッケージのインストール
install.packages(c("dplyr", "MatchIt", "survival", "survminer", 
                   "tableone", "ggplot2", "lubridate"))

# ライブラリの読み込み
library(dplyr)
library(MatchIt)      # 傾向スコアマッチング
library(survival)     # 生存分析
library(tableone)     # 記述統計表作成
library(ggplot2)      # 可視化

-- 急性期脳卒中患者の抽出 SELECT * FROM dpc_data WHERE -- 脳梗塞、脳出血、くも膜下出血 (SUBSTR(dpc_code,1,6) IN ('010060', '010040', '010030')) AND age >= 18 -- 成人のみ AND admission_type IN (1,2,3) -- 緊急・救急入院 AND los >= 1 -- 在院日数1日以上 AND discharge_status IS NOT NULL; -- 転帰情報あり

# 在院日数カテゴリ data$los_category <- cut(data$length_of_stay, breaks = c(0, 7, 14, 30, Inf), labels = c("短期", "中期", "長期", "超長期")) # 30日再入院フラグ data$readmission_30d <- ifelse(data$days_to_readmission <= 30 & !is.na(data$days_to_readmission), 1, 0) # 在宅復帰率 data$home_discharge <- ifelse(data$discharge_destination == "home", 1, 0) # リハビリテーション強度 data$rehab_intensity <- data$total_rehab_units / data$length_of_stay # 併存症指標の作成（Charlson Comorbidity Index） data$cci_score <- data$dm + data$copd + data$chf + data$cad + data$ckd*2 + data$cancer*3 # 重み付けスコア

# 傾向スコアモデルの構築 ps_model <- glm(early_rehab ~ age + gender + severity_score + cci_score + hospital_volume + stroke_type + admission_jcs, family = binomial(link = "logit"), data = stroke_data) # 傾向スコアの算出 stroke_data$propensity_score <- predict(ps_model, type = "response") # 傾向スコア分布の確認 library(ggplot2) ggplot(stroke_data, aes(x = propensity_score, fill = factor(early_rehab))) + geom_histogram(alpha = 0.7, position = "identity", bins = 30) + labs(title = "傾向スコアの分布", x = "傾向スコア", y = "度数")

library(MatchIt) # 基本的な1:1マッチング matched_data <- matchit(early_rehab ~ age + gender + severity_score + cci_score + hospital_volume + stroke_type + admission_jcs, data = stroke_data, method = "nearest", distance = "logit", caliper = 0.25 * sd(stroke_data$propensity_score)) # マッチング結果の確認 summary(matched_data) plot(matched_data) # マッチ後データの抽出 matched_df <- match.data(matched_data)

# マッチング前後の共変量バランス評価 library(tableone) # マッチング前 table1_before <- CreateTableOne(vars = covariates, strata = "early_rehab", data = stroke_data, test = FALSE) print(table1_before, smd = TRUE) # マッチング後 table1_after <- CreateTableOne(vars = covariates, strata = "early_rehab", data = matched_df, test = FALSE) print(table1_after, smd = TRUE) # Love plotによる可視化 library(cobalt) love.plot(matched_data, stats = c("mean.diffs", "variance.ratios"), thresholds = c(m = 0.1, v = 2), binary = "std", abs = TRUE, var.order = "unadjusted")

library(survival) library(survminer) # 生存データの作成 stroke_data$time_to_readmission <- ifelse(stroke_data$readmission_30d == 1, stroke_data$days_to_readmission, 30) stroke_data$readmission_event <- stroke_data$readmission_30d # Kaplan-Meier推定 km_fit <- survfit(Surv(time_to_readmission, readmission_event) ~ early_rehab, data = matched_df) # 生存曲線の描画 ggsurvplot(km_fit, data = matched_df, pval = TRUE, conf.int = TRUE, risk.table = TRUE, xlab = "退院後日数", ylab = "再入院回避率", title = "早期リハビリの再入院予防効果") # Log-rank検定 survdiff(Surv(time_to_readmission, readmission_event) ~ early_rehab, data = matched_df)

# 単変量Cox回帰 cox_simple <- coxph(Surv(time_to_readmission, readmission_event) ~ early_rehab, data = matched_df) summary(cox_simple) # 多変量Cox回帰（必要に応じて残存交絡の調整） cox_adjusted <- coxph(Surv(time_to_readmission, readmission_event) ~ early_rehab + age + severity_score, data = matched_df) summary(cox_adjusted) # 比例ハザード性の確認 cox.zph(cox_adjusted) # ハザード比の可視化 forest_model(cox_adjusted)

# 入院費用の算出 calculate_hospitalization_cost <- function(data) { # DPC包括点数 + 出来高部分 data$dpc_cost <- data$dpc_points * 10 # 1点=10円 data$fee_for_service_cost <- data$additional_procedures * 10 data$total_cost <- data$dpc_cost + data$fee_for_service_cost # 患者負担と保険者負担の分離 data$patient_cost <- data$total_cost * data$copayment_rate data$insurance_cost <- data$total_cost - data$patient_cost return(data) } # QALYの推定（EQ-5Dデータがない場合のマッピング） estimate_qaly <- function(data) { # 在院日数と重症度からQOL低下を推定 qol_baseline <- 0.8 # 健康時のQOL qol_acute <- 0.3 # 急性期のQOL qol_recovery_rate <- 0.1 # 日当たり回復率 data$qaly_loss <- (qol_baseline - qol_acute) * data$length_of_stay / 365.25 return(data) }

calculate_icer <- function(intervention_group, control_group) { # 平均費用の差 cost_diff <- mean(intervention_group$total_cost) - mean(control_group$total_cost) # 平均効果の差（QALY） effect_diff <- mean(intervention_group$qaly_gain) - mean(control_group$qaly_gain) # ICER計算 icer <- cost_diff / effect_diff # ブートストラップによる信頼区間 bootstrap_icer <- replicate(1000, { boot_int <- sample_n(intervention_group, nrow(intervention_group), replace = TRUE) boot_ctrl <- sample_n(control_group, nrow(control_group), replace = TRUE) boot_cost_diff <- mean(boot_int$total_cost) - mean(boot_ctrl$total_cost) boot_effect_diff <- mean(boot_int$qaly_gain) - mean(boot_ctrl$qaly_gain) boot_cost_diff / boot_effect_diff }) return(list( icer = icer, ci_lower = quantile(bootstrap_icer, 0.025, na.rm = TRUE), ci_upper = quantile(bootstrap_icer, 0.975, na.rm = TRUE) )) }