Zenn
Aidemy Tech BlogPublicationへの投稿
↔️

限界値のあるデータをベイズモデリングで効率的に扱う

Shingo Nakazawa
2023/12/15に公開
Julia
データ分析
Python 3
pymc
ベイズモデリング
tech

はじめに - 限界値のあるデータとは

こんにちは。株式会社アイデミーデータサイエンティストの中沢(@shnakazawa_ja)です。

データ分析の世界では、完全な情報を得られない状況に頻繁に直面します。その一つとして 「限界値のあるデータ」 が挙げられます。例えば水質調査や医療診断、食品安全性評価などで機器の検出限界を超えた測定値が得られたときにこうしたデータが出てきます。具体的には以下のようなものですね。

サンプルID 汚染物質濃度 (mg/L)
1 0.008
2 0.007
3 0.006
4 <0.005
5 0.010
6 0.008
7 <0.005
8 0.009
9 0.011
10 <0.005

ここではいくつかのサンプルで汚染物質の濃度が検出限界 "0.005 mg/L" 未満になっています。
このようなデータでは実際の数値が確定しておらず、分析の際にはこの不確定性をどのように扱うかが鍵となります。

多くの場合は値を置換するアプローチが取られます。目的達成に向けて「検出限界未満ならゼロと考えてOK」であれば全てゼロに置換してしまうのが簡単ですし、「検出限界未満なら検出限界と考える」が適切であれば検出限界と同値に揃えるのがよいでしょう。「検出限界未満であるという事実が重要」と考える場合は新しいラベルの作成も有効かもしれません。しかしながら、いずれの方法も閾値未満の数値を一定の値に固定してしまうことで情報量を落としてしまっています[1]。
そこで本稿では、ベイズモデリングを用いて検出限界未満であるという情報を効率的に使う方法を紹介します[2]。

仮説をモデル式で表現する

本稿では簡単のためにシンプルなモデルを考えます。
真の平均値 \mu にノイズ \sigma が乗り、潜在的な測定値 y が生成されている、と仮説を置きましょう。
そして、検出下限値を T とすると、y が T 以上であれば測定値にはそのまま y が記録されますが、y が T 未満であれば <T と記録されると考えます[3]。

この仮説をモデル式で表現すると以下のようになります。

  • 測定できているデータについて (observedのobs)：
y[n] \sim Normal(\mu, \sigma)
y[n] ≥ T
n = 1, 2, ..., N_{obs}　
  • 検出下限未満のデータについて (thresholdのthr)：
y[n] \sim Normal(\mu, \sigma)　　
y[n] < T
n = 1, 2, ..., N_{thr}　

データを2群に分けて扱うことがポイントになります。

測定ごとの尤度を考える

次に、それぞれのモデル式において尤度を求めます。

測定できているデータについて、尤度は単純に

L(\mu,\sigma) = Normal(y|\mu,\sigma)

です。

検出下限未満のデータについて、尤度はyがT未満になる確率ですので

L(\mu,\sigma) = Probability[y<T]

これを変形していくと

Probability[y<T] = \int^T_{-\infin}Normal(y|\mu, \sigma) \\[15pt] =　\int^T_{-\infin}\frac{1}{\sqrt{2\pi\sigma^2}}exp[-\frac{1}{2}(\frac{y-\mu}{\sigma})^2]dy \\[15pt] = \int^\frac{T-\mu}{\sigma}_{-\infin}\frac{1}{\sqrt{2\pi}}exp[-\frac{1}{2}z^{2}]dz \\[15pt] = \int^\frac{T-\mu}{\sigma}_{-\infin}\phi(z)dz \\[15pt] = \Phi(\frac{T-\mu}{\sigma})

となります（ここで \phi は標準正規分布の確率密度関数で、 \Phi はその累積分布関数）。
すなわち、検出下限未満のデータの尤度は標準正規分布の累積分布関数として表現できるということです。よって、上式の最後の値の対数を取って、測定できているデータの対数尤度に合算すれば、検出下限未満の値を人為的に置換せずとも、その情報をモデルに組み込むことができます。

PyMCでの実装と結果

PyMCでの実装を見てみましょう。

import pandas as pd
import pymc as pm
import arviz as az

df = pd.DataFrame({
    "サンプルID": [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
    "汚染物質濃度 (mg/L)": [0.008, 0.007, 0.006, "<0.005", 0.010, 0.008, "<0.005", 0.009, 0.011, "<0.005"]
})

# 測定できているデータのみのDataFrame
df_obs = df[df["汚染物質濃度 (mg/L)"] != "<0.005"]
# 検出下限を下回ったレコードの数
n_thr = len(df[df["汚染物質濃度 (mg/L)"] == "<0.005"])

with pm.Model() as model:
    # パラメータの事前分布
    mu = pm.Normal('mu', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=10)

    # 測定できているデータの対数尤度
    y_obs_rv = pm.Normal.dist(mu=mu, sigma=sigma)
    y_obs_logp = pm.logp(y_obs_rv, df_obs["汚染物質濃度 (mg/L)"].astype(float)).sum()
    

    # 検出下限未満のデータの対数尤度
    prob_T = pm.logcdf(pm.Normal.dist(mu=mu, sigma=sigma), 0.005)
    
    #　対数尤度を足す
    pm.Potential('y_obs_logp', y_obs_logp)
    pm.Potential('prob_T', n_thr * prob_T)

with model:
    idata = pm.sample(1000)
    
summary = az.summary(idata, round_to=4)
  • pm.logpで測定できているデータの対数尤度を計算し、sum()で総和を取得
  • pm.logcdfで正規分布の累積分布関数の対数、すなわち検出下限未満のデータ1サンプル分の対数尤度を取得
  • pm.Potentialで取得した対数尤度を合算

として本稿の内容を実現しています。

実行結果は以下のようになります。

mu (平均値) に着目してみると0.0066という値が得られており、

  • 検出下限未満のデータを下限に揃えときの平均値は0.0074
  • ゼロとみなすと0.0059

ですので、 限界値を超えたデータを置換したときと比較し、より "現実を反映した" 結果が得られていると考えられます[4]。

Turing.jlでの実装と結果

(2024-07-05追加)

# 限界値のあるデータを扱う
using Turing
using Distributions
using DataFrames

# データの準備
df = DataFrame(
    サンプルID = 1:10,
    汚染物質濃度 = [0.008, 0.007, 0.006, missing, 0.010, 0.008, missing, 0.009, 0.011, missing]
)

# 測定できているデータのみのDataFrame
df_obs = dropmissing(df, :汚染物質濃度)
# 検出下限を下回ったレコードの数
n_thr = count(ismissing, df.汚染物質濃度)

# モデルの定義
@model function contamination_model(y_obs, n_thr)
    # パラメータの事前分布
    mu ~ Normal(0, 10)
    sigma ~ truncated(Cauchy(0, 10), 0, Inf)  # HalfNormalの代わりにtruncated Cauchyを使用

    # 測定できているデータの対数尤度
    for i in 1:length(y_obs)
        y_obs[i] ~ Normal(mu, sigma)
    end

    # 検出下限未満のデータの対数尤度
    for i in 1:n_thr
        Turing.@addlogprob! logcdf(Normal(mu, sigma), 0.005)
    end
end

# モデルのインスタンスを作成
model = contamination_model(df_obs.汚染物質濃度, n_thr)

# サンプリング
chain = sample(model, NUTS(), 1000)

# 結果の表示
println(chain)

さいごに

本稿では限界値のあるデータの不確定性を効果的に扱う方法を紹介しました。
ベイズモデリングではデータの背後にある構造を分布として捉えるからこそ、このような形でデータの不確定性を表現できます。本手法の考え方は他の分析にも応用ができるかもしれません。

本稿が皆様のデータ分析のお役に立てば幸いです。

参考文献

更新履歴

  • 2024-07-05
    • Julia での実装例を追加
脚注

  1. 限界値のあるデータに限らず、分析で最も大切なのは目的を明確にすること。目的を達成できるのであれば情報の損失があってもシンプルな方が良いでしょう。 ↩︎

  2. 上限を超えるケースでも考え方は同じです。 ↩︎

  3. 限界値 (Limit) には大文字Lを使うのが慣例のようなのですが、尤度 (Likelihood) のLと被るため本稿ではT (Threshold) を用いました。 ↩︎

  4. sigma (標準偏差) に着目しても、下限に揃えると0.0021、ゼロとみなすと0.0041で、ベイズモデリングの結果は中間値を示しています。 ↩︎

Aidemy Tech Blog
Aidemy Tech BlogPublication

株式会社アイデミーの技術ブログです。

Discussion