クイックノート

ちょっとした発見・アイデアから知識の発掘を

【R】傾向スコアマッチング

ある薬の効き目を調べたいと思った時に、
薬が投与された人と薬が投与されていない人の状態を比べるというのは、
自然なことですよね。

ところが注意しないといけないのは、この場合だと、
薬を投与する人を偏りなく選ばなければならないということです。

もし、症状の重い人に優先して投与して、
症状の軽い人にはあまり投与しないで、
投与の有無と、健康状態のデータを取るとどうなるでしょうか?

症状の重い人の方が、薬の投与に関わらず健康状態が悪化する可能性が高いのは当然ですね。
つまり、症状の重い人を狙って薬を投与すると、
*(症状が重いから)→薬を投与した→健康状態が悪化した *
となり、まるで、薬を投与した結果、健康状態が悪化したかのようなデータが取れてしまいます。

「症状が重い人に優先して投薬を行う」というのは自然なことだと考えられるので、
むしろ、このような偏りを持っているデータの方が多いでしょう。

逆に、このような偏りのないデータを得ようと思うと、
ちゃんとコントロールされた状況下で新しくデータを取り直さないといけませんが、
それはかなり手間がかかりそうです。

では、偏りを持ったデータを使って、
ある処置の効果の有無を調べることはできるのでしょうか?

その方法の一つが傾向スコアマッチングと呼ばれる方法です。

傾向スコアマッチングとは

傾向スコアマッチングでは、
データの中である処置が行われやすい状況と、
行われにくい状況が混在している時、
その行われやすさによる偏りを除くようにして、
データの分析を行う方法です。

単純に、処置を行われた人と行われなかった人を比べると、
上で述べたように、処置の行われやすさの違いによって生じる偏りの影響を受けてしまいます。

では、この偏りをなくすにはどうすれば良いかというと、
処置の行われやすさが同じくらいの人で比べればいい
ということになります。

傾向スコアマッチングでは、
この処置の行われやすさ(傾向スコア)が同じくらいのもの同士を
マッチングして、比較することで、
処置の行われやすさの偏りに影響されない比較を行うという方法です。

処置の行われやすさを計算する

傾向スコアマッチングで、一つの重要なステップは、
処置の行われやすさ、つまり傾向スコアを求めるという部分です。

上の例だと、薬が投与されやすい人とされにくい人をどう見分けるかという話になります。

これには、データの中の他の様々な属性から
薬が投与される確率を回帰によって求める
という方法がとられます。

例えば、症状の重さを示す数値が高い人ほど、
薬が投与される確率は高いなどの関係を回帰分析で求めておき、
それぞれの人に、薬が投与され得るであろう確率を計算するのです。

マッチング

上で計算した確率(傾向スコア)を使って、
その処置が行われる確率によるデータの偏りを除くために、
確率を揃えたデータで比較を行います。

この処理はマッチングと呼ばれます。

その方法は、同じ傾向スコアを持つデータをサンプルしてくるという単純なものです。
ただし、処置が行われたものと行われなかったもので比較したいので、
同じ傾向スコアを持つけど、「処置が行われたデータ」と「処置が行われなかったデータ」を、
ペアにしてデータから取り出してきます。

この組を比べていくことで、
処置の行われやすさを揃えた上で、
その処置が行われたことによる効果を見る
ことができるのです。

ただし、必ずしも全てのデータについて、ペアの相手が見つかるとは限りません。
極端に処置が行われやすい傾向スコアを持つデータはほとんど全てが処置済みのデータとなるでしょうから、
相手となる処置が行われないデータが見つかりません。
このようにマッチングに失敗したデータは分析から除外されてしまうので、
集めたデータの中の一部のみが分析に使われることとなります。

R で実践

それでは、Rで実際に傾向スコアマッチングをおこなってみましょう。

準備

傾向スコアマッチングを行うために、
R のパッケージMatchingを利用します。

インストールは

install.packages("Matching")

で行います。

インストール後は、パッケージを読み込んでおきましょう。

library(Matching)

データ生成

傾向スコアマッチングを試すためにサンプルデータを生成します。

ここでは、症状の重い人と症状のない人の2つのグループがあって、
症状の重い人に投薬すると死亡率が軽減して、
症状のない人に投薬すると死亡率に変化がないようなデータを作成します。

データは、各行が1人分のデータを表していて、
各列が次のような属性を保持しているものとします。

treated(処置が行われたか否か) dead(死亡したか否か) value(症状の重さを表す数値)

tread,deadTRUE,FALSEの2値として、
value正規分布に従う乱数を生成します。
正規分布の平均値は、症状の重さに応じて設定していて、
症状が重い人は平均が1、
症状のない人は平均が-1となるようにデータを生成しています。

データを生成するコードは次の通りです。

make_data = function(p,N){
  x = sapply(1:N,function(i){
    r = runif(1)
    if(p > r){ # weak
      r = runif(1)
      if(r < 0.8){ # treated
        r = runif(1)
        if(r < 0.5){ # alive
          c(TRUE,FALSE,rnorm(1,mean=1))
        }else{       # death
          c(TRUE,TRUE,rnorm(1,mean=1))
        }
      }else{ # non treated
        r = runif(1)
        if(r < 0.3){
          c(FALSE,FALSE,rnorm(1,mean=1))
        }else{
          c(FALSE,TRUE,rnorm(1,mean=1))
        }
      }
    }else{ # normal
      r = runif(1)
      if(r < 0.2){ # treated
        r = runif(1)
        if(r < 0.8){ # alive
          c(TRUE,FALSE,rnorm(1,mean=-1))
        }else{
          c(TRUE,TRUE,rnorm(1,mean=-1))
        }
      }else{ # non treated
        if(r < 0.8){ # alive
          c(FALSE,FALSE,rnorm(1,mean=-1))
        }else{
          c(FALSE,TRUE,rnorm(1,mean=-1))
        }
      }
    }
  })
  rownames(x) = c("treated","death","value")
  data.frame(t(x))
}

x = make_data(0.1,10000)

長くなってしまいましたが、条件で場合わけをしているだけの簡単なデータです。
Nがデータの数で、pが症状の重い人の割合をさしていて、
症状の重い人には80%の確率で処置が行われ、
症状のない人には20%の確率で処置が行われます。

また、症状の重い人に処置が行われると死亡率が70%から50%に軽減され、
症状のない人の死亡率は処置の有無にかかわらず20%としています。
細かい値はコードを書き換えることで変更できます。

傾向スコアの計算

それでは、それぞれのデータに対して傾向スコアを求めましょう。

ここではロジスティック回帰を使って、傾向スコアを計算しています。

  prop = glm(treated~value+death,
          family=binomial(link="logit"),
          data=x)
  score = prop$fitted.values

処置がされたかどうかを、その他の属性値から回帰して、
処置のされた確率をscoreとして計算しています。

マッチング

上で求めた傾向スコアを使って、
データのマッチングを行い、
処置のされたものとされていないものでの死亡率の差を比較します。

これは、Matchingパッケージの関数Matchを呼び出すだけで実行できます。

  res = Match(Y = x$death,
          Tr = x$treated,
          X = score)

Y には処置による効果を見たい属性(今の場合は死亡か否か)、
Tr には処置の有無を表す属性、
Xには傾向スコアを指定します。

結果

それでは傾向スコアマッチングの結果をみてみましょう。

単純な群の比較

傾向スコアマッチングのありがたみを知るために、
まずは、傾向スコアマッチングを行わず、
単純に全データに対して、処置を行ったものと行わなかったものの死亡率を比較してみましょう。

tmp = table(x[,c("treated","death")])
non_treated_death = tmp[1,2] / sum(tmp[1,])
treated_death = tmp[2,2] / sum(tmp[2,])
print("non_tr_death   tr_death   diff")
print(sprintf("%f       %f   %f",non_treated_death,treadted_death,treated_death-non_treated_death)) 

結果は、

[1] "non_tr_death   tr_death   diff"
[1] "0.257603       0.380952   0.050609"

このように、処置を行った方が死亡率が高いという結果になりました。

実際には処置をすることで、死亡率が軽減するようにデータを生成しているので、
この分析結果は正しくないですよね。

傾向スコアマッチングの結果

それでは、傾向スコアマッチングの結果を見てみましょう。

summary(res)

結果は、

Estimate...  -0.017217 
AI SE......  0.012401 
T-stat.....  -1.3884 
p.val......  0.16502 

Original number of observations..............  10000 
Original number of treated obs...............  2667 
Matched number of observations...............  2667 
Matched number of observations  (unweighted).  35485 

Estimateが処置を行った場合と行わなかった場合の死亡率の差ですが、
これが負の値なので、処置を行うことで死亡率が軽減するという
データにあった分析結果が得られていることが分かります。

まとめ

データのある性質の起こりやすさによる影響を除外して、
その性質による特徴を調べる方法として傾向スコアマッチングを紹介しました。

主に医療研究の分野を中心に使われることの多い分析方法ですが、
他の分野でも偏りを除いた分析は重要になってくるので、
応用の広い分析方法ですね。

プライバシーポリシー