frmgham-homework-kuroda
Hiroyuki Kuroda
2023-05-30
HDS2023授業「臨床予測モデル」最終課題
フラミンガムデータを用いた予測モデルの構築
プログラム作成者：黒田浩行
【課題】
使用するデータはフラミンガム研究のデータといたします。
各自でこのデータで構築できるような予測モデルについてのリサーチクエスチョンを考えて、生存時間分析（ロジスティック回帰分析でも可）で予測モデルを構築してください。アウトカムは、通常はデータセットにある変数、ANGINA,
HOSPMI, MI_FCHD, ANYCHD, STROKE, CVD, HYPERTEN,
DEATHのいずれかになると思いますが、独自にアウトカム定義を作成していただいてもかまいません。リサーチクエスチョン、用いる変数、変数選択の方法、連続変数の扱い、欠測値の扱い、統計解析、対象集団の背景、全体の当てはまり（Brier
score、予測値と観測値の平均の比較、R2等々、）、判別能、Calibration、内的妥当性（Bootstrap,もしくは交差検証）、作成したモデル式、結果の簡単な考察（1枚程度）、を含めるようにしてください。
0. 構築モデルの設定
アウトカムをCVDとDEATHのコンポジットアウトカムとする。
20年後のイベント（コンポジットアウトカム）発生リスクを予測する。
1. データクリーニング
# ライブラリのロード
pacman::p_load(tidyverse,
rio,
here,
rms,
mice,
naniar,
gtsummary,
corrr,
survival,
survminer,
ggcorrplot,
ggpubr,
broom,
glmnet)# シードの固定
set.seed(2023)# データのインポート
dat <- import(here("frmgham2.csv"))ピリオドごとにデータが格納されているが、まずは1行1参加者のデータに加工する必要あり。
# IDのみ抽出
dat_ID <- dat %>%
dplyr::select(RANDID) %>%
unique()ピリオドごとのデータを抽出する。3つと少ないので手作業で。
dat_IDとdat_1のN数が同じ（コードは省略）なので、ピリオド2以降に初めて参加した人はいない。
どのピリオドに参加したか分かるよう、ダミー変数を作成する。
後のコードで、例えばPERIOD_2=1ならピリオドに2参加あり、0なら参加なし、というようにする。
dat_1 <- dat %>%
filter(PERIOD == 1) %>%
mutate(PERIOD = 1) # あとで使うダミー変数のもと
dat_2 <- dat %>%
filter(PERIOD == 2) %>%
mutate(PERIOD = 1) # 同上
dat_3 <- dat %>%
filter(PERIOD == 3) %>%
mutate(PERIOD = 1) # 同上変数の性質及び定義から、ピリオドで変わりうるものを選ぶ。
Ageも変わるが、ベースラインだけ分かればよい。
これらの変数はどのピリオドに由来するか分かるよう、末尾に数字を付与する。重複するので、性別など不変の変数はdat_2,
3から落とす。
# ピリオドで変わりうる変数
variable_vars <- c("PERIOD", "SYSBP", "DIABP", "BPMEDS", "CURSMOKE", "CIGPDAY",
"TOTCHOL", "HDLC", "LDLC", "BMI", "GLUCOSE", "DIABETES",
"HEARTRTE", "PREVAP", "PREVCHD", "PREVMI", "PREVSTRK", "PREVHYP")
# ピリオドごとに末尾に_1, _2, _3を付与する
dat_1 <- dat_1 %>%
rename_at(vars(all_of(variable_vars)), ~paste(., "_1", sep = ""))
dat_2 <- dat_2 %>%
dplyr::select(RANDID, all_of(variable_vars)) %>%
rename_at(vars(all_of(variable_vars)), ~paste(., "_2", sep = ""))
dat_3 <- dat_3 %>%
dplyr::select(RANDID, all_of(variable_vars)) %>%
rename_at(vars(all_of(variable_vars)), ~paste(., "_3", sep = ""))くっつけて一つのデータセットにする。
結果的にはdat_IDは不要（ピリオド1のデータで十分）。
参加ピリオドを示すダミー変数を1/0にする（ここでのNAは欠測ではなくて、単に参加していないだけなので）
# replace_naはlistに収めれば複数指定できる
dat_master <- dat_ID %>%
full_join(dat_1, by = "RANDID") %>%
full_join(dat_2, by = "RANDID") %>%
full_join(dat_3, by = "RANDID") %>%
dplyr::select(RANDID, PERIOD_1, PERIOD_2, PERIOD_3, everything()) %>%
replace_na(list(PERIOD_1 = 0, PERIOD_2 = 0, PERIOD_3 = 0))コンポジットアウトカム(CVD or DEATH)の変数EVENTを作成する。
EVENTまでの時間TIMEは、TIMECVDとTIMEDEATHのうち、小さい方をとればいい。
# CVD, DEATH, TIMECVD, TIMEDTHは欠測なし
# 下記コードで確認可能
# sum(is.na(dat_master$CVD), is.na(dat_master$DEATH),
# is.na(dat_master$TIMECVD), is.na(dat_master$TIMEDTH))
dat_master <- dat_master %>%
rowwise() %>%
mutate(EVENT = if_else((CVD==1 | DEATH==1), 1, 0)) %>%
mutate(TIME = min(TIMECVD, TIMEDTH)) %>%
ungroup()ピリオド2以降でも測定している検査値などは、本来は時間依存性変数として扱うべきと思われる。
ここまで頑張ったが、今回はベースライン（ピリオド1）のデータのみを使用してモデルを作成することにする。使用する可能性のあるデータを、整形しつつ選ぶ。
【重要】ベースラインでCHD, strokeの既往（PREVAP, PREVCHD, PREVMI,
PREVSTRK）がある者は、本研究から除外する。
ピリオド1ではLDLとHDLは測っていない？ようなので除く。
dat_use <- dat_master %>%
filter(PREVAP_1==0 & PREVCHD_1==0 & PREVMI_1==0 & PREVSTRK_1==0) %>% # CHD, strokeの既往者は除外
dplyr::select(EVENT, TIME, SEX, AGE, SYSBP_1, DIABP_1, BPMEDS_1,
CURSMOKE_1, CIGPDAY_1, TOTCHOL_1, BMI_1, GLUCOSE_1, DIABETES_1, PREVHYP_1) %>%
mutate(across(everything(), as.numeric)) %>%
rename(SBP = SYSBP_1,
DBP = DIABP_1,
BPMED = BPMEDS_1,
SMOKE = CURSMOKE_1,
CIGPDAY = CIGPDAY_1,
TCHO = TOTCHOL_1,
BMI = BMI_1,
GLU = GLUCOSE_1,
DIABETES = DIABETES_1,
PREHT = PREVHYP_1) %>%
mutate(SEX = SEX - 1) %>% # 不要かもだがSEXは1/0にする, 0=male, 1=female
mutate(across(c(SEX, BPMED, SMOKE, DIABETES, PREHT), factor)) # fator型にする変数# rcs描画のためのデータ範囲を指定
ddist <- datadist(dat_use)
options(datadist="ddist") 2. 変数の特徴をつかむ
gtsummary::tbl_summaryで、baseline characteristicsをまとめる。
tbl_1 <- dat_use %>%
dplyr::select(-EVENT, -TIME) %>% # EVENT, TIMEはいらないので除く
tbl_summary(missing_text = "Missing") %>% # 欠測値の表示方法
add_n() %>%
modify_header(label ~ "**Variable**") %>%
modify_caption("**Table1. Baseline Characteristics**") %>%
bold_labels()
tbl_1| Variable | N | N = 4,2151 |
|---|---|---|
| SEX | 4,215 | |
| 0 | 1,810 (43%) | |
| 1 | 2,405 (57%) | |
| AGE | 4,215 | 49 (42, 56) |
| SBP | 4,215 | 128 (117, 144) |
| DBP | 4,215 | 82 (75, 90) |
| BPMED | 4,163 | |
| 0 | 4,046 (97%) | |
| 1 | 117 (2.8%) | |
| Missing | 52 | |
| SMOKE | 4,215 | |
| 0 | 2,127 (50%) | |
| 1 | 2,088 (50%) | |
| CIGPDAY | 4,186 | 0 (0, 20) |
| Missing | 29 | |
| TCHO | 4,165 | 234 (206, 263) |
| Missing | 50 | |
| BMI | 4,198 | 25.4 (23.1, 28.0) |
| Missing | 17 | |
| GLU | 3,828 | 78 (71, 87) |
| Missing | 387 | |
| DIABETES | 4,215 | |
| 0 | 4,107 (97%) | |
| 1 | 108 (2.6%) | |
| PREHT | 4,215 | |
| 0 | 2,917 (69%) | |
| 1 | 1,298 (31%) | |
| 1 n (%); Median (IQR) | ||
連続変数間の相関を確認する。
cor <- dat_use %>%
dplyr::select(where(is.numeric)) %>% # 連続変数だけを選ぶ
dplyr::select(-EVENT, -TIME) %>%
drop_na() %>% # ここでは欠測を有する行は考えないことにする
correlate() %>%
shave()## Correlation computed with
## • Method: 'pearson'
## • Missing treated using: 'pairwise.complete.obs'cor## # A tibble: 7 × 8
## term AGE SBP DBP CIGPDAY TCHO BMI GLU
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 AGE NA NA NA NA NA NA NA
## 2 SBP 0.388 NA NA NA NA NA NA
## 3 DBP 0.201 0.784 NA NA NA NA NA
## 4 CIGPDAY -0.191 -0.0929 -0.0581 NA NA NA NA
## 5 TCHO 0.263 0.215 0.169 -0.0305 NA NA NA
## 6 BMI 0.134 0.326 0.378 -0.0913 0.114 NA NA
## 7 GLU 0.118 0.133 0.0605 -0.0566 0.0447 0.0889 NA相関をプロットしてみる→相関高めなのは血圧同士くらい。
corresult <- dat_use %>%
dplyr::select(where(is.numeric)) %>%
dplyr::select(-EVENT, -TIME) %>%
drop_na() %>%
cor(method = "pearson")
corrplot <- ggcorrplot(corr = corresult, hc.order = FALSE, method = "square", title = "cor plot",
colors = c("#4b61ba", "white", "red"), lab = TRUE)
corrplot
欠測を確認する：GLU, BPMED, TCHO, CIGPDAY, BMIのみ。あとでmiceしていく。
naplot_1 <- gg_miss_var(dat_use, show_pct = TRUE)
naplot_1
naplot_2 <- vis_miss(dat_use)
naplot_2
3. 粗解析
観察時間のヒストグラムを描画する。
2500/4215 =
6割くらいはイベントフリーだったことがわかる。
time_hist <- ggplot(data = dat_use) +
geom_histogram(aes(x = TIME)) +
theme_classic()
time_hist## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
overallのKM曲線
survfit_1 <- survfit(Surv(time = TIME, event = EVENT) ~ 1,
data = dat_use)
KM_1 <- ggsurvplot(survfit_1, data = dat_use)
KM_1
KM曲線～SEXごと→差あり
survfit_2 <- survfit(Surv(time = TIME, event = EVENT) ~ SEX,
data = dat_use)
KM_2 <- ggsurvplot(survfit_2, data = dat_use, pval = TRUE)
KM_2
KM曲線～SMOKEごと→視覚的には意外と大差ない（log-rankは<0.05）
survfit_3 <- survfit(Surv(time = TIME, event = EVENT) ~ SMOKE,
data = dat_use)
KM_3 <- ggsurvplot(survfit_3, data = dat_use, pval = TRUE)
KM_3
KM曲線～PREHTごと→差あり
survfit_4 <- survfit(Surv(time = TIME, event = EVENT) ~ PREHT,
data = dat_use)
KM_4 <- ggsurvplot(survfit_4, data = dat_use, pval = TRUE)
KM_4
全部の変数を入れて、比例ハザード性を検討してみる。
ここではcomplete caseで考える。
global
p=0.076であり、比例ハザード性は保たれている（保たれていないとは言えない）と判断される。
formura_1 <- Surv(time = TIME, event = EVENT) ~
SEX + AGE + SBP + DBP + BPMED + SMOKE + CIGPDAY + TCHO + BMI + GLU + DIABETES + PREHT
cph_1 <- coxph(formura_1, data = dat_use)
zph_1 <- cox.zph(cph_1)
zph_1## chisq df p
## SEX 0.0401 1 0.841
## AGE 1.0620 1 0.303
## SBP 4.9742 1 0.026
## DBP 2.6929 1 0.101
## BPMED 0.1408 1 0.707
## SMOKE 0.0239 1 0.877
## CIGPDAY 2.0579 1 0.151
## TCHO 4.1324 1 0.042
## BMI 1.1030 1 0.294
## GLU 0.2012 1 0.654
## DIABETES 0.4666 1 0.495
## PREHT 3.2491 1 0.071
## GLOBAL 19.5534 12 0.076zph_plot_1 <- ggcoxzph(zph_1)
zph_plot_1
連続変数とハザードが直線的に関連しているのか、検討する。
BMIとGLUは、U-shapeと考えられる。
# 連続変数の変数名をまとめる
con_var <- c("AGE", "SBP", "DBP", "CIGPDAY", "TCHO", "BMI", "GLU")
# プロット結果をまとめるリストを用意する
plot <- list()
for (x in con_var){
# cphに投入するformulaを文字列で、for文で順番に指定していく
formula_tmp <- as.formula(paste("Surv(time = TIME, event = EVENT) ~ rcs(", x, ", 4)"))
fit_tmp <- cph(formula_tmp, data = dat_use)
# plot(Predict())はfor文内で実行できないようなので、eval(parse())で無理やり回した（非推奨）。
# Predict関数の呼び出しを文字列として作成し、それをパースして評価する
plot_cmd <- paste("plot(Predict(fit_tmp, ", x, "))")
plot_tmp <- eval(parse(text = plot_cmd))
plot[[x]] <- plot_tmp
}
# 結果をggarangeでまとめて表示する
logHR_plot <- ggarrange(plotlist = plot, ncol = 3, nrow = 3)
logHR_plot
4. 欠測補完・変数選択
MICEで欠測値の代入を行う。
→補完データセットそれぞれに対し、stepwiseで変数選択する。
→50%以上のデータセットで採用された変数を、最終モデルに採用する変数にする。
# 多重補完で作成するデータセットの数は20にする。
m <- 20
# MICEデータの作成
dat_mice <- mice(dat_use,
m = m,
maxit = 5,
print = FALSE,
seed = 2023)多重代入後の確認を行う。
CIGPDAYはやや分布がずれている。
plot_mice <- plot(dat_mice)
plot_mice
densityplot_mice <- densityplot(dat_mice)
densityplot_mice
変数選択に使用するフォーミュラについて。
リスク予測算出式を導出する際に複雑になってしまうため、logHR_plotの形状から、二乗項を投入することを考える。rcs項はリスク予測算出式を導出するには複雑と考えた。
formula_cph <- Surv(TIME, EVENT) ~
SEX + AGE + SBP + DBP + BPMED + SMOKE + CIGPDAY + TCHO + pol(BMI, 2) + pol(GLU, 2) + DIABETES + PREHT変数選択を行っていく。
# もともとの候補変数
vars_original <- c("SEX", "AGE", "SBP", "DBP", "BPMED", "SMOKE", "CIGPDAY", "TCHO", "BMI", "GLU", "DIABETES", "PREHT")
# miceの各々のcompleteデータで採用された変数名を貯めていく
# 全データで採用されたらm=20回名前が載る
vars_stack <- c()
for (id in 1:m){
dat_tmp <- complete(dat_mice, id)
cph_tmp <- cph(formula_cph, data = dat_tmp)
bw_tmp <- fastbw (cph_tmp, type = "individual" , rule = "aic")
vars_stack <- c(vars_stack, bw_tmp$names.kept)
}採用率50%以上で最終モデルに採用することにしている。 DBPが一度も採用されず。BMIは採用されないこともあったが、50%(=10回)以上なので採用する。
vars_counts <- sort(table(vars_stack), decreasing = TRUE)
vars_counts## vars_stack
## AGE BPMED CIGPDAY DIABETES GLU PREHT SBP SEX
## 20 20 20 20 20 20 20 20
## SMOKE TCHO BMI
## 20 19 15DBPのみ一度も採用されなかったことを念のため確認
drop_vars <- setdiff(vars_original, unique(vars_stack))
drop_vars## [1] "DBP"さらに、shrinkage methodであるLASSOでも変数選択してみる。
vars_stack_lasso <- c()
for (id_lasso in 1:m){
dat_tmp_lasso <- complete(dat_mice, id_lasso)
surv_tmp_lasso <- Surv(time = dat_tmp_lasso$TIME, event = dat_tmp_lasso$EVENT)
model_tmp_lasso <- model.matrix(~ SEX + AGE + SBP + DBP + BPMED + SMOKE + CIGPDAY + TCHO +
pol(BMI, 2) + pol(GLU, 2) + DIABETES + PREHT,
data = dat_tmp_lasso)
# ハイパーパラメータlamdaがベストな時の変数を選ぶことにする。
# 5-fold cvでbest lambdaを検索する。
# alpha=1でLASSO
cvfit_tmp <- cv.glmnet(model_tmp_lasso, surv_tmp_lasso,
nfolds = 5,
family = "cox",
alpha = 1)
# best lmbdaをとりだし
best_lambda_tmp <- cvfit_tmp$lambda.min
# best lambdaを使用
best_model_tmp <- glmnet(model_tmp_lasso, surv_tmp_lasso,
family = "cox",
lambda = best_lambda_tmp)
# betaをとりだし（選択されて0になったものもある可能性あり）
coef_tmp <- coef(best_model_tmp, s = best_lambda_tmp)
# 切片項以外の、採用された変数（beta=0にならなかった変数）を格納する。
# [-1]は切片項を除いたことを、!=0はbetaが0でない（残った）変数をピックしている。
vars_stack_lasso <- c(vars_stack_lasso, rownames(coef_tmp)[-1][coef_tmp[-1] != 0])
}意外とDBPも採用された。
vars_counts_lasso <- sort(table(vars_stack_lasso), decreasing = TRUE)
vars_counts_lasso## vars_stack_lasso
## AGE BPMED1 CIGPDAY DIABETES1
## 20 20 20 20
## pol(BMI, 2)BMI^2 pol(GLU, 2)GLU PREHT1 SBP
## 20 20 20 20
## SEX1 SMOKE1 TCHO DBP
## 20 20 20 18LASSOではDBPも残ったが、臨床的意義やSBPとの相関の高さを考えて、DBPを落としたモデルを採用する。
formula_final <- Surv(TIME, EVENT) ~
SEX + AGE + SBP + BPMED + SMOKE + CIGPDAY + TCHO + pol(BMI, 2) + pol(GLU, 2) + DIABETES + PREHT5. 判別能の評価
5-fold CVとbootstrapで内的検証する。
miceで生成したcomplete
dataそれぞれに5-fold CV もしくは bootstrapを適用し、平均をとる(Rubin’s
ruleによる点推定)。なおbootstrapの実行時にcoefがinfでは、とエラーが出るが、実行内容を確認してもこれは問題ないと考えた。
# 結果を格納する
Cindex_original <- c()
Cindex_CV <- c()
Cindex_BS <- c()
options(warn = -1) # 警告を非表示
for (j in 1:m){
dat_tmp <- complete(dat_mice, j)
cph_tmp <- cph(formula_final, data = dat_tmp, x = TRUE, y = TRUE)
# 5-fold CV
CV_tmp <-validate(cph_tmp, method = "crossvalidation", B = 5, bw = FALSE, rule = "aic")
# bootstrap(100リサンプリング)
BS_tmp <-validate(cph_tmp, method = "boot", B = 100, bw = FALSE, rule = "aic")
# original dataのDxyを、Cindexに換算しつつとりだす（apparent Cindex）
Cindex_original <- c(Cindex_original, (CV_tmp[1, 1]*0.5 + 0.5))
# CVのDxyを、Cindexに換算しつつとりだす
Cindex_CV <- c(Cindex_CV, (CV_tmp[1, 3]*0.5 + 0.5))
# BSのDxyを、Cindexに換算しつつとりだす
Cindex_BS <- c(Cindex_BS, (BS_tmp[1, 5]*0.5 + 0.5))
}
options(warn = 0) # 警告を表示# apparent C indexを表示する
print_Cindex_original <- paste("apparent C-index is", round(mean(Cindex_original), 3))
print_Cindex_original## [1] "apparent C-index is 0.725"# 5-fold CVで検証したC incexを表示する
print_Cindex_CV <- paste("C-index corrected by 5-fold CV is", round(mean(Cindex_CV), 3))
print_Cindex_CV## [1] "C-index corrected by 5-fold CV is 0.723"# bootstrapで検証したC incexを表示する
print_Cindex_BS <- paste("C-index corrected by 100-resampled bootstrap is", round(mean(Cindex_BS), 3))
print_Cindex_BS## [1] "C-index corrected by 100-resampled bootstrap is 0.723"apparent/CV/BS C-indexに大差はなく、重大な過学習はないものと思われる。 miceによる結果を統合する（with→poolはcph使えないので、coxphで）。
cox_mi <- with(data = dat_mice,
coxph(Surv(TIME, EVENT) ~
SEX + AGE + SBP + BPMED + SMOKE + CIGPDAY + TCHO + pol(BMI, 2) + pol(GLU, 2) + DIABETES + PREHT))
pool_cox_mi <- pool(cox_mi)
# 結果をまとめる
result_cox_mi <- as.data.frame(summary(pool_cox_mi)) %>%
dplyr::select(term, estimate, std.error, p.value) %>%
mutate(est_exp = estimate,
CIlow = estimate - 1.96*std.error,
CIhigh = estimate + 1.96*std.error) %>%
mutate(across(c(est_exp, CIlow, CIhigh), exp)) %>% # expしてHRとして解釈できるようにする
dplyr::select(term, estimate, est_exp, CIlow, CIhigh)
result_cox_mi## term estimate est_exp CIlow CIhigh
## 1 SEX1 -6.295780e-01 0.5328166 0.4795286 0.5920262
## 2 AGE 6.707503e-02 1.0693757 1.0626812 1.0761124
## 3 SBP 1.040788e-02 1.0104622 1.0075290 1.0134040
## 4 BPMED1 4.344539e-01 1.5441196 1.2326348 1.9343162
## 5 SMOKE1 1.958333e-01 1.2163241 1.0509407 1.4077334
## 6 CIGPDAY 1.006390e-02 1.0101147 1.0042321 1.0160318
## 7 TCHO 9.167567e-04 1.0009172 0.9998255 1.0020100
## 8 pol(BMI, 2)BMI -3.821756e-02 0.9625035 0.8975699 1.0321347
## 9 pol(BMI, 2)BMI^2 8.085762e-04 1.0008089 0.9996341 1.0019851
## 10 pol(GLU, 2)GLU 6.527683e-03 1.0065490 1.0010994 1.0120284
## 11 pol(GLU, 2)GLU^2 -8.471995e-06 0.9999915 0.9999774 1.0000057
## 12 DIABETES1 4.326627e-01 1.5413562 1.1551049 2.0567648
## 13 PREHT1 2.424008e-01 1.2743049 1.1165199 1.45438786. 較正能の評価
どの時点を評価するかを指定する（今回は20年）。
u <- 365.25*20# stacked datasetsを作る
dat_mice_stack <- data.frame()
for (k in 1:m){
dat_mice_stack <- rbind(dat_mice_stack, complete(dat_mice, k))
}
# 作成したデータセット分重みづけする
dat_mice_stack$weight <- 1/m
# miceデータ（をスタックしたデータ）でcph
# クラスターを考慮した分散ではないので、信頼区間は狭く推定される
# ここでは点推定値だけ検討している
cph_stack <- cph(formula_final, data = dat_mice_stack,
weights = weight,
x = TRUE,
y = TRUE)
# 予測
estim_surv <- survest(cph_stack,
newdata = dat_mice_stack,
times = u)$surv
ext_val <- val.surv(cph_stack,
newdata = dat_mice_stack,
u = u,
est.surv = estim_surv,
S = Surv(dat_mice_stack$TIME, dat_mice_stack$EVENT))
# プロット
plot(ext_val,
xlab = "Predicted probability at 20 years",
ylab = "Observed probability at 20 years",
las = 1)
7. リスク算出式の導出
6.に引き続き、20年後のリスク算出とする。
まず、ベースラインに対するbetaを取り出す。
beta <- result_cox_mi$estimate %>%
round(., 5)
names(beta) <- result_cox_mi$term
beta## SEX1 AGE SBP BPMED1
## -0.62958 0.06708 0.01041 0.43445
## SMOKE1 CIGPDAY TCHO pol(BMI, 2)BMI
## 0.19583 0.01006 0.00092 -0.03822
## pol(BMI, 2)BMI^2 pol(GLU, 2)GLU pol(GLU, 2)GLU^2 DIABETES1
## 0.00081 0.00653 -0.00001 0.43266
## PREHT1
## 0.24240ベースラインの生存関数S(t=20y)を推定する。
# ベースライン変数
baseline <- data.frame(SEX = 0,
AGE = 0,
SBP = 0,
BPMED = 0,
SMOKE = 0,
CIGPDAY = 0,
TCHO = 0,
BMI = 0,
GLU = 0,
DIABETES = 0,
PREHT = 0)
# miceデータをまとめる用
baseline_surv_20y_mice <- c()
# miceデータそれぞれに、formura_finalをフィットさせる。
# 得られたモデルをbaselineに適用する。
# そこから20年後のベースラインの生存関数S(t=20y)を推定する。
for (l in 1:m){
dat_base_tmp <- complete(dat_mice, l)
cph_base_tmp <- cph(formula_final, data = dat_base_tmp, x = TRUE, y = TRUE)
survfit_base_tmp <- survfit(cph_base_tmp, newdata = baseline)
survprob_20y <- summary(survfit_base_tmp, times = u)$surv
baseline_surv_20y_mice <- c(baseline_surv_20y_mice, survprob_20y)
}
# 最後に統合する（点推定値なので平均）
# AGE=0, BMI=0などなので、これ自体は臨床的な意義はない。
baseline_surv_20y <- round(mean(baseline_surv_20y_mice), 5)
print_baseline_surv_20y <- paste("baseline 20-year survival rate is", baseline_surv_20y*100, "%")
print_baseline_surv_20y## [1] "baseline 20-year survival rate is 99.742 %"S(t|X)=S(t)^exp(β*X)である。
ここから任意の変数セットに対応する、20年後のイベント発生リスク算出式を導出する。
# 小さい値を指数表記にしないようにする
options(scipen = 1000)
# betaを再掲
beta## SEX1 AGE SBP BPMED1
## -0.62958 0.06708 0.01041 0.43445
## SMOKE1 CIGPDAY TCHO pol(BMI, 2)BMI
## 0.19583 0.01006 0.00092 -0.03822
## pol(BMI, 2)BMI^2 pol(GLU, 2)GLU pol(GLU, 2)GLU^2 DIABETES1
## 0.00081 0.00653 -0.00001 0.43266
## PREHT1
## 0.24240よって、求めたい20年後のイベント発生リスクは、1-生存確率なので、下記となる。
20-year event risk (%) = (1- 0.99742^exp[A]) * 100%, where
A =
- 0.62958 [if female]
+ 0.06708*age in years
+ 0.01041*SBP
in mmHg
+ 0.43445 [if BP meds is used]
+ 0.19583 [if current
smoker]
+ 0.01006*number of cigarette per day
+
0.00092*total cholesterol in mmol/L
- 0.03822*BMI
+
0.00081*BMI^2
+ 0.00653*fasting blood glucose
-
0.00001*fasting blood glucose^2
+ 0.43266 [if diabetes]
+
0.24240 [if past history of HT]
例えば、
SEX=1, AGE=60, SBP=120, BPMED=1, SMOKE=1, CIGPDAY=10,
TCHO=150, BMI=20, GLU=90, DIABETES=0, PREHT=1
とすると、
test_1 <- c(1, 60 ,120, 1, 1, 10, 150, 20, 20^2, 90, 90^2, 0, 1)
test_2 <- beta*test_1
test_3 <- sum(test_2)
# 求めたいリスクは下記。
test_risk <- round((1 - 0.99742^exp(test_3))*100, 2)
print_test_risk <- paste("test 20-year event risk is", test_risk, "%")
print_test_risk## [1] "test 20-year event risk is 58.2 %"