👌

事例と実装で学ぶ因果推論〜ダイエットアプリの効果を検証するには?(第3回:回帰モデル)

2023/04/24に公開

こんにちは、Ubie Discoveryのmatsu-ryuです。

「事例と実装で学ぶ因果推論〜ダイエットアプリの効果を検証するには?」シリーズ3回目の記事になります。

前回は、「ポテンシャルアウトカムフレームワーク」について、解説を行いました。
https://zenn.dev/matsu_ryu/articles/359b16f145abf7

観察データの効果検証を行うための「因果推論」に必要な、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の効果が測定することを試みます。

回帰モデルによる因果推論

分析の手順を示すと

outcome=a[0]+ a[1]\cdot treatment+ a[2]\cdot age + a[3]\cdot gender + a[4]\cdot weight + a[5]\cdot diet+ε(誤差項)

の回帰モデルにおいて、a[0]〜a[5]について最適なパラメータを求めます。
ダイエットアプリの効果は、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がa[1]に対応します。

              Estimate Std. Error t value Pr(>|t|)    
treatment   -2.0082445  0.0200480 -100.17   <2e-16 ***

上記のとおり、「-2.0082445」となっており、真の効果の「-2」に近い値が推定できています。

この回帰モデルの係数の解釈を数式で説明します。

E[Y|T=1,L]=a[0]+ a[1]+ a[2]\cdot age + a[3]\cdot gender + a[4]\cdot weight + a[5]\cdot diet

E[Y|T=0,L]=a[0]+ a[2]\cdot age + a[3]\cdot gender + a[4]\cdot weight + a[5]\cdot diet

E[Y|T=1,L]-E[Y|T=0,L] = a[1]
となります。言葉で説明すると
L(性別、年齢、体重、食習慣)が同一の集団における、ダイエットアプリの平均的な効果はa[1]

手順(アウトカム〜処置変数+交絡因子の回帰モデルを実行する)だけ示せば、非常に単純ですが、この推定が有効なのは、回帰モデルの仮定を満たしている場合であり、その仮定を満たさない場合、推定が上手くいきません。疑似データの生成方法を変更することで、回帰モデルの仮定を満たさない仮想的な状況を作り、回帰モデルでは、真の効果を推定できない状況を見ていきたいと思います。

回帰モデルが適用できない状況

回帰モデルの仮定は、いつくかあるのですが、
代表的な仮定は、以下の式の表現の通り

outcome=a[0]+ a[1]\cdot treatment+ a[2]\cdot age + a[3]\cdot gender + a[4]\cdot weight + a[5]\cdot diet+ε(誤差項)
  • 線形性(アウトカムと説明変数の関係は線形和)を満たす
  • パラメータaは一定
    • 例えば年齢が20でも30でも、アプリの処置効果: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から直接応募していただいても良いですし、一度カジュアルにお話したい方はそれぞれのメンバーと話せるカジュアル面談ページがあるのでぜひお気軽にご連絡ください。

https://recruit.ubie.life/jd_dev/eng_dataanalyst
https://recruit.ubie.life/jd_dev/eng_analytics
https://recruit.ubie.life/jd_dev/eng_data
https://youtrust.jp/recruitment_posts/66409c1aa5e042b45e1e581b6357dcd9
https://youtrust.jp/recruitment_posts/c19393ba65b2921c764e35a760269869

最近soの制度もupdateされました。より魅力的なsoに生まれ変わってます!
https://prtimes.jp/main/html/rd/p/000000051.000048083.html

Ubie テックブログ

Discussion