こんにちは、Ubie Discoveryのmatsu-ryuです。
「事例と実装で学ぶ因果推論〜ダイエットアプリの効果を検証するには?」シリーズ3回目の記事になります。
前回は、「ポテンシャルアウトカムフレームワーク」について、解説を行いました。
観察データの効果検証を行うための「因果推論」に必要な、3条件
- 交換可能性：Exchangebility
- 正値性：Positivity
- 一貫性：Consistency
について解説し、第1回記事の「問題設定」で生成したダイエットアプリの疑似データは、この3条件を満たしていることを述べました。
3回目の今回は、「回帰モデル」によって、ダイエットアプリの効果を測定する手法について解説していきます。
疑似データのおさらい
- age:年齢
- gender：性別
- weight：体重
- diet：食習慣（0:良い、1:普通、2:悪い）
- treatment：アプリ利用有無
- outcome_0：アプリ利用無しの場合の体重増減（ポテンシャルアウトカム）
- outcome_1：アプリ利用有りの場合の体重増減（ポテンシャルアウトカム）
- outcome：実際に観測されたの体重増減
| age | gender | weight | diet | treatment | outcome_0 | outcome_1 | outcome |
|---|---|---|---|---|---|---|---|
| 50 | Male | 77.55406 | 0 | 0 | 0.8668771 | -1.1331229 | 0.8668771 |
| 34 | Female | 47.39310 | 0 | 0 | 2.4679330 | 0.4679330 | 2.4679330 |
| 33 | Male | 58.68436 | 2 | 0 | -0.3344975 | -2.3344975 | -0.3344975 |
| 22 | Female | 33.45681 | 0 | 0 | 2.5217297 | 0.5217297 | 2.5217297 |
| 61 | Male | 83.06569 | 2 | 1 | -0.4964269 | -2.4964269 | -2.4964269 |
| 69 | Female | 55.62998 | 0 | 1 | 2.5858146 | 0.5858146 | 0.5858146 |
回帰モデルによる効果検証
効果検証したいことを、改めて整理していきましょう。
前回記事「DAGと交絡因子」で述べた通り、問題は、 L((性別,年代,体重,食習慣))がT((アプリ利用))にも、Y((体重増減))にも影響を及ぼしていることです。このような変数を、交絡因子とい言います。
下記のDAGのイメージのように、Lの影響を排除して、TによるYの効果が測定することを試みます。
回帰モデルによる因果推論
分析の手順を示すと
の回帰モデルにおいて、
ダイエットアプリの効果は、treatmentの係数、a[1]になります。
では、実際にコードを実行してみましょう。
- 回帰分析の実行
# 回帰分析
> model <- lm(outcome ~ age + gender + weight + diet + treatment, data = data)
- 結果
> summary(model)
Call:
lm(formula = outcome ~ age + gender + weight + diet + treatment,
data = data)
Residuals:
Min 1Q Median 3Q Max
-0.80175 -0.13676 -0.00118 0.13772 0.76131
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.5798537 0.0414984 86.27 <2e-16 ***
age 0.0212767 0.0005022 42.36 <2e-16 ***
genderMale -0.4329524 0.0146359 -29.58 <2e-16 ***
weight -0.0425844 0.0006541 -65.10 <2e-16 ***
diet -0.8615105 0.0114992 -74.92 <2e-16 ***
treatment -2.0082445 0.0200480 -100.17 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2075 on 994 degrees of freedom
Multiple R-squared: 0.987, Adjusted R-squared: 0.987
F-statistic: 1.512e+04 on 5 and 994 DF, p-value: < 2.2e-16
treatmentのEstimateが
Estimate Std. Error t value Pr(>|t|)
treatment -2.0082445 0.0200480 -100.17 <2e-16 ***
上記のとおり、「-2.0082445」となっており、真の効果の「-2」に近い値が推定できています。
この回帰モデルの係数の解釈を数式で説明します。
となります。言葉で説明すると
L（性別、年齢、体重、食習慣）が同一の集団における、ダイエットアプリの平均的な効果はa[1]
手順（アウトカム〜処置変数+交絡因子の回帰モデルを実行する）だけ示せば、非常に単純ですが、この推定が有効なのは、回帰モデルの仮定を満たしている場合であり、その仮定を満たさない場合、推定が上手くいきません。疑似データの生成方法を変更することで、回帰モデルの仮定を満たさない仮想的な状況を作り、回帰モデルでは、真の効果を推定できない状況を見ていきたいと思います。
回帰モデルが適用できない状況
回帰モデルの仮定は、いつくかあるのですが、
代表的な仮定は、以下の式の表現の通り
- 線形性（アウトカムと説明変数の関係は線形和）を満たす
- パラメータaは一定
- 例えば年齢が２０でも３０でも、アプリの処置効果:a[1]は一定
という仮定があります。
回帰モデルの仮定を満たさないとはどういう状況でしょうか？
具体的な例を見ていきましょう。
ダイエットアプリの例を取ると、同じアプリ利用をしていても、性別、年代、体重、食事習慣が異なる場合、その効果は、一定ではないはずです。また、体重が多いほど減量はしやすく、少なければ減量は厳しくなってきます。減量と体重の関係はの一定の比例関係ではないことが想像できます。以下のようなイメージで体重増減のデータを生成するとどうなるでしょうか？
| 体重 | 食事 | ダイエットアプリ利用 | 体重増減 |
|---|---|---|---|
| 50kg未満 | 良い | あり | -2kg |
| 50kg台 | 良い | あり | -3kg |
| 60kg台 | 良い | あり | -4kg |
| 70kg台 | 良い | あり | -8kg |
| 80kg台 | 良い | あり | -10kg |
| 90kg以上 | 良い | あり | -20kg |
| 50kg未満 | 良い | なし | 0kg |
| 50kg台 | 良い | なし | -1kg |
| 60kg台 | 良い | なし | -3kg |
| 70kg台 | 良い | なし | -4kg |
| 80kg台 | 良い | なし | -5kg |
| 90kg以上 | 良い | なし | -6kg |
※一部のイメージを抜粋
もともと体重がある人で、食習慣も良く、ダイエットアプリしている人がより減量幅が大きくなるように設定します。
ここからは、回帰モデルの仮定を満たさない状況で疑似データを作成し
回帰モデルの推定が、真の効果とずれる状況をコードで示していきます。
library(tibble)
# 年齢、性別、体重、食習慣を生成する関数
generate_covariate <- function(n) {
age <- sample(20:69, n, replace = TRUE)
gender <- sample(c("Male", "Female"), n, prob = c(0.5, 0.5), replace = TRUE)
weight <- rnorm(n, 50, 10) + 0.1 * age + 10 * (gender == "Male") + 0.3 * rnorm(n)
diet <- sample(0:2, n, prob = c(0.4, 0.3, 0.3), replace = TRUE)
tibble(age, gender, weight, diet)
}
#アプリ利用有無を生成する関数
generate_treatment <- function(covariate, a) {
logit <- a[1] * covariate$age + a[2] * (covariate$gender == "Male") + a[3] * covariate$weight + a[4] * covariate$diet
prob <- exp(logit) / (1 + exp(logit))
treatment <- rbinom(nrow(covariate), 1, prob)
tibble(treatment)
}
# アウトカムを生成する関数
generate_outcome <- function(covariate, treatment, outcome_name) {
outcome <- covariate %>%
mutate( outcom =
case_when(
treatment == 1 & weight < 50 & diet == 0 ~ -2,
treatment == 1 & weight < 60 & diet == 0 ~ -3,
treatment == 1 & weight < 70 & diet == 0 ~ -4,
treatment == 1 & weight < 80 & diet == 0 ~ -8,
treatment == 1 & weight < 90 & diet == 0 ~ -10,
treatment == 1 & weight >= 90 & diet == 0 ~ -15,
treatment == 1 & weight < 50 & diet == 1 ~ -2,
treatment == 1 & weight < 60 & diet == 1 ~ -3,
treatment == 1 & weight < 70 & diet == 1 ~ -4,
treatment == 1 & weight < 80 & diet == 1 ~ -5,
treatment == 1 & weight < 90 & diet == 1 ~ -8,
treatment == 1 & weight >= 90 & diet == 1 ~ -10,
treatment == 1 & weight < 50 & diet == 2 ~ 0,
treatment == 1 & weight < 60 & diet == 2 ~ -2,
treatment == 1 & weight < 70 & diet == 2 ~ -3,
treatment == 1 & weight < 80 & diet == 2 ~ -4,
treatment == 1 & weight < 90 & diet == 2 ~ -5,
treatment == 1 & weight >= 90 & diet == 2 ~ -6,
treatment == 0 & weight < 50 & diet == 0 ~ 0,
treatment == 0 & weight < 60 & diet == 0 ~ -1,
treatment == 0 & weight < 70 & diet == 0 ~ -3,
treatment == 0 & weight < 80 & diet == 0 ~ -2,
treatment == 0 & weight < 90 & diet == 0 ~ -4,
treatment == 0 & weight >= 90 & diet == 0 ~ -4,
treatment == 0 & weight < 50 & diet == 1 ~ 2,
treatment == 0 & weight < 60 & diet == 1 ~ 3,
treatment == 0 & weight < 70 & diet == 1 ~ 4,
treatment == 0 & weight < 80 & diet == 1 ~ 3,
treatment == 0 & weight < 90 & diet == 1 ~ 5,
treatment == 0 & weight >= 90 & diet == 1 ~ 8,
treatment == 0 & weight < 50 & diet == 2 ~ 2,
treatment == 0 & weight < 60 & diet == 2 ~ 3,
treatment == 0 & weight < 70 & diet == 2 ~ 6,
treatment == 0 & weight < 80 & diet == 2 ~ 6,
treatment == 0 & weight < 90 & diet == 2 ~ 7,
treatment == 0 & weight >= 90 & diet == 2 ~ 10,
T ~ 1
)
) %>% .$outcom
tibble(!!outcome_name := outcome)
}
# 疑似データを生成する関数
generate_data <- function(n, a, b) {
covariate <- generate_covariate(n)
treatment <- generate_treatment(covariate, a)
# ポテンシャルアウトカム
outcome_0 <- generate_outcome(covariate, rep(0, nrow(covariate)), 'outcome_0')
outcome_1 <- generate_outcome(covariate, rep(1, nrow(covariate)), 'outcome_1')
# 観測値
outcome <- generate_outcome(covariate, treatment$treatment, 'outcome')
cbind(covariate, treatment, outcome_0, outcome_1, outcome)
}
# 1000件の疑似データを生成する
n <- 1000
a <- c(-0.1, 1, 0.02, 1)
set.seed(456)
data <- generate_data(n, a, b)
# 真の結果
data <-data %>% mutate(
true_effect = outcome_1-outcome_0
)
第1回からの主な変更は、generate_outcome関数です。
このデータの真のダイエットアプリ利用効果は
> mean(data$true_effect)
[1] -4.96
です。
ナイーブ（単純）な比較の場合は、
> library(tableone)
> # tableoneを利用して生成されたデータの傾向を比較する
> vars <- c("age", "gender", "weight", "diet", "outcome")
> CreateTableOne(vars = vars, data = data, factorVars = c("gender", "diet"),
+ strata = c("treatment"))
Stratified by treatment
0 1 p test
n 772 228
age (mean (SD)) 47.63 (13.94) 33.61 (11.04) <0.001
gender = Male (%) 341 (44.2) 148 (64.9) <0.001
weight (mean (SD)) 59.21 (11.31) 60.36 (10.60) 0.173
diet (%) <0.001
0 347 (44.9) 50 (21.9)
1 216 (28.0) 65 (28.5)
2 209 (27.1) 113 (49.6)
outcome (mean (SD)) 1.30 (2.97) -3.20 (1.80) <0.001
-3.2 - (1.3)=-4.5
と真の値(-4.5)を捉えられていません。（第1回で説明した状況と同じ）
回帰モデルで推定した場合は、どうなるでしょうか？
> # 回帰分析
> model <- lm(outcome ~ age + gender + weight + diet + treatment, data = data)
> summary(model)
Call:
lm(formula = outcome ~ age + gender + weight + diet + treatment,
data = data)
Residuals:
Min 1Q Median 3Q Max
-3.7246 -1.3906 -0.0972 1.2837 6.4914
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.3108526 0.3423048 -0.908 0.364
age 0.0002678 0.0040952 0.065 0.948
genderMale 0.1833914 0.1214686 1.510 0.131
weight -0.0107978 0.0054042 -1.998 0.046 *
diet 2.6217523 0.0657083 39.900 <2e-16 ***
treatment -5.7090012 0.1471010 -38.810 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.695 on 994 degrees of freedom
Multiple R-squared: 0.742, Adjusted R-squared: 0.7407
F-statistic: 571.6 on 5 and 994 DF, p-value: < 2.2e-16
treatmentのEstimateは-5.7kgと、真の値（-4.96）から、ややズレている状況です。
このズレ幅の影響は、扱う事業や問題設定によりますが、ここで言いたいことは、回帰分析はモデルの仮定が多く、扱う変数が多いほどモデルを誤設定しやすい ということです。
このようなデータに対して、効果検証するにはどうしたら良いでしょうか？
次回
今回は、因果推論するための前提条件を満たし、線形モデルで生成した疑似データにおいて回帰モデルで効果検証できることを示しました。一方、回帰モデルの仮定を崩した疑似データにおいては、回帰モデルによる推定値が真の値からズレる状況を示しました。
次回は、このような状況にフレキシブルに対応できる「傾向スコアモデリング」について、解説していきます。
引き続き、因果推論の手法について学んでいきましょう！
一緒に働く仲間を探しています！
Ubie Discoveryは、「テクノロジーで適切な医療に導く」ことをミッションとしています。
このミッションをデータドリブンで実現することを目指しています。
そのためには、データアナリスト、アナリティクスエンジニア、データエンジニアの力が必要です。
そこで、私たちは一緒に働く仲間を募集しています!
興味のある方は下記JDから直接応募していただいても良いですし、一度カジュアルにお話したい方はそれぞれのメンバーと話せるカジュアル面談ページがあるのでぜひお気軽にご連絡ください。
最近soの制度もupdateされました。より魅力的なsoに生まれ変わってます！
Discussion