【学習動機】
費用対効果分析をする際の感度分析で, モデルのパラメータを確率分布で乱数生成する確率的感度分析(Probabilistic Sensitivity Analysis; PSA)を行うことがある. パラメータが何かしらの確率を示している場合, ベータ分布(もしくはディリクレ分布)が使われることが多い. しかし, ベータ分布のパラメータ\(\alpha\)と\(\beta\)の設定が雑に見える分析がときどきある. そのような設定では, ベータ分布の確率変数の期待値が\(E[X]=\alpha/(\alpha+\beta)\)に基本分析(Base Case)の値が合うように設定されているのだが, そのスケールが無視されている. たとえば, \(p=0.2\)がBase Caseの分析のとき, PSAでは\(\alpha=2, \beta=8\)だったり, \(\alpha=20, \beta=80\)だったり, \(\alpha=0.2, \beta=0.8\)だったりする. 平均値だけ合わせているのである.
\(X～Beta(\alpha,\beta)\), \(Y～Beta(k\alpha,k\beta)\)のとき, \(E[X]=E[Y]\)だが, \(Var[X]=\frac{\alpha \beta}{(\alpha+\beta)^2 (\alpha+\beta+1)}\)である一方, \(Var[Y]=\frac{\alpha \beta}{(\alpha+\beta)^2 (k\alpha+k\beta+1)}=V[X]*\frac{\alpha+\beta+1}{k\alpha+k\beta+1}\)である.
二つのパラメータを平均値だけ合わせる形になっているので, 分散が定まらないのは当然な結果である. このような設定方法をしているとき, PSAの結果はおかしなものになる可能性がないのか気になった.
というわけで, 今回は簡単な費用対効果分析のモデルを設定し, ベータ分布のパラメータを平均値は同じだが分散が変わるように変化させることで, PSAの結果がどう変化するかを確認する. コードはRで書く.
【学習内容】
簡単な決定木モデルを考える.
イメージとしては, ある疾患Dに対する予防治療Aを受けるか否かの意思決定について, 期待費用の観点から費用を抑制するか否かで状況である. ここでは, 簡単のため, 費用のみを扱う.
Strategy Aは, 予防治療Aを受ける. 受けると疾患Dになる確率\(p_B\)の\(RR_{AB}\)倍. 予防治療Aの費用は\(cost_A\)かかる. 疾患になった際は治療費用\(cost_D\)かかる.
Strategy Bは, 予防治療Aを受けない. 疾患Dになる確率は\(p_B\). 疾患になった際は治療費用\(cost_D\)かかる.
基本分析(Base Case)では, \(p_B\)は20%, \(RR_AB\)は0.8, \(cost_D\)は500, \(cost_A\)は10とする. その後, 確率的感度分析(PSA)では, 簡単のため, \(p_B\)のみベータ分布で動かして, 費用を抑制できる割合を確認する. 以下でRコードでプログラムを回す.
パッケージ読み込みと乱数固定.
# パッケージ読み込みと乱数固定
library(tidyverse)
library(ggplot2)
set.seed(2024)パラメータの設定(Base Case)
\(p_B\)は20%, \(RR_AB\)は0.8, \(cost_D\)は500, \(cost_A\)は10とする.
# BaseCase
p_B <- 0.2
RR_AB <- 0.8
p_A <- p_B*RR_AB
cost_D <- 500
cost_A <- 10期待費用の計算(Base Case)
# 期待費用の計算
Ecost_A <- p_A*(cost_D+cost_A) + (1-p_A)*cost_A
Ecost_A## [1] 90Ecost_B <- p_B*cost_D + (1-p_B)*0
Ecost_B## [1] 100Ecost_AvsB <- Ecost_A - Ecost_B
Ecost_AvsB## [1] -10Aの期待費用は90, Bの期待費用は100なので, 期待値的にはAの方が費用を抑えられる.
PSAの関数作成
nsimはシミュレーション回数. 今回はnsim=1000で回していく. par_Beta_kはベータ分布パラメータの2と8をそれぞれ何倍するかの係数. 下でk=1, 10, 0.1で試す.
# PSA(Probabilistic Sensitivity Analysis)の関数作成
PSA_beta <- function(nsim, par_Beta_k){
# Beta分布のパラメータ設定
par_Beta_alpha <- 2*par_Beta_k
par_Beta_beta <- 8*par_Beta_k
# 結果を保存するベクトル作成
v_Ecost_A <- numeric(nsim)
v_Ecost_B <- numeric(nsim)
v_Ecost_AvsB <- numeric(nsim)
cost_D <- 500
cost_A <- 10
for(i in 1:nsim){
p_B <- rbeta(n=1, shape1=par_Beta_alpha, shape2=par_Beta_beta)
p_A <- p_B*RR_AB
# 期待費用の計算
Ecost_A <- p_A*(cost_D+cost_A) + (1-p_A)*cost_A
Ecost_B <- p_B*cost_D + (1-p_B)*0
Ecost_AvsB <- Ecost_A - Ecost_B
# 結果の保存
v_Ecost_A[i] <- Ecost_A
v_Ecost_B[i] <- Ecost_B
v_Ecost_AvsB[i] <- Ecost_AvsB
}
df_res <- data.frame(
Ecost_A = v_Ecost_A,
Ecost_B = v_Ecost_B,
Ecost_AvsB = v_Ecost_AvsB
)
return(df_res)
}期待費用の計算(k=1でPSA)
k <- 1
a <- 2*k
b <- 8*kベータ分布の分散確認
# p_Bの分散
v_p_B <- a*b/(a+b)^2/(a+b+1)
v_p_B## [1] 0.01454545ベータ分布の密度関数の形状確認
# ベータ分布の密度関数を表示
df_beta <- data.frame(x = seq(0, 1, length.out = 100)) %>%
mutate(y = dbeta(x, a, b))
ggplot(df_beta, aes(x = x, y = y)) +
geom_line(color = "red") +
labs(title = "ベータ分布の密度関数",
x = "x",
y = "density") +
theme_minimal()
PSAを回して, 期待費用の差分のプロット.
# PSA
df_k <- PSA_beta(nsim=1000 ,par_Beta_k=1)
ggplot(df_k, aes(x=Ecost_AvsB)) +
geom_histogram(colour="grey", bins=30) +
geom_vline(xintercept=0, colour="red") +
theme_bw()
費用が抑制できる割合
# 費用が抑制できる割合
df_k %>% filter(Ecost_AvsB<0) %>% nrow() / df_k %>% nrow() *100## [1] 77.4期待費用の計算(k=10でPSA)
k <- 10
a <- 2*k
b <- 8*kベータ分布の分散確認
# p_Bの分散
v_p_B <- a*b/(a+b)^2/(a+b+1)
v_p_B## [1] 0.001584158ベータ分布の密度関数の形状確認
# ベータ分布の密度関数を表示
df_beta <- data.frame(x = seq(0, 1, length.out = 100)) %>%
mutate(y = dbeta(x, a, b))
ggplot(df_beta, aes(x = x, y = y)) +
geom_line(color = "red") +
labs(title = "ベータ分布の密度関数",
x = "x",
y = "density") +
theme_minimal()
PSAを回して, 期待費用の差分のプロット.
# PSA
df_k <- PSA_beta(nsim=1000 ,par_Beta_k=10)
ggplot(df_k, aes(x=Ecost_AvsB)) +
geom_histogram(colour="grey", bins=30) +
geom_vline(xintercept=0, colour="red") +
theme_bw()
費用が抑制できる割合
# 費用が抑制できる割合
df_k %>% filter(Ecost_AvsB<0) %>% nrow() / df_k %>% nrow() *100## [1] 99.9期待費用の計算(k=0.1でPSA)
k <- 0.1
a <- 2*k
b <- 8*kベータ分布の分散確認
# p_Bの分散
v_p_B <- a*b/(a+b)^2/(a+b+1)
v_p_B## [1] 0.08ベータ分布の密度関数の形状確認
# ベータ分布の密度関数を表示
df_beta <- data.frame(x = seq(0, 1, length.out = 100)) %>%
mutate(y = dbeta(x, a, b))
ggplot(df_beta, aes(x = x, y = y)) +
geom_line(color = "red") +
labs(title = "ベータ分布の密度関数",
x = "x",
y = "density") +
theme_minimal()
PSAを回して, 期待費用の差分のプロット.
# PSA
df_k <- PSA_beta(nsim=1000 ,par_Beta_k=0.1)
ggplot(df_k, aes(x=Ecost_AvsB)) +
geom_histogram(colour="grey", bins=30) +
geom_vline(xintercept=0, colour="red") +
theme_bw()
費用が抑制できる割合
# 費用が抑制できる割合
df_k %>% filter(Ecost_AvsB<0) %>% nrow() / df_k %>% nrow() *100## [1] 38.1結果とまとめ
k=1, 10, 0.1と変えたとき, 費用が抑制できている割合は77.4%, 99.9%, 38.1%と変わった. 結構結果が変わった. モデルに重要なパラメータがこのような設定をされていると, 普通PSAをするときはモデルのパラメータ全体を確率的に回すので気が付きにくいが, 実はPSAの結果を大きく変えてしまっているかもしれない. 結論としては, モデルの確率のパラメータのベータ分布は分散(または, 二次モーメント)も考慮して, 確率的感度分析を回すべき.
【学習予定】
費用対効果分析の話を学んで, 何か気付いたら今後も記録したいところ.
