完全情報最尤推定法(FIML)のお気持ちを掴んで、選抜効果を考慮した相関の推定を行う

この記事は何

この記事は、完全情報最尤推定法(FIML)を用いて、欠測データに対して簡単な仮説検証を行なった記事です。

欠測が起こっているデータを用いて解析を行う場合、単純に欠測が起こっているデータを削除して解析を行うと誤った結果を導いてしまう場合がある。例えば、

• 入学試験の結果(X)と学年末試験の結果(Y)の相関
• 新卒採用試験の結果(X)とその後のパフォーマンスの相関(Y)
  

など、XとYの関係を調べたいときに、単純な回帰分析

$$ Y = \alpha + \beta X $$
を行って係数$\beta$を見たり、相関係数

$$\rho_{X, Y} = \frac{cov(X, Y)}{\sigma_{X} \sigma_{Y}}$$

の値を用いて判断しようとすると、真の偏回帰係数と比べて$\beta$がおかしくなったり、真の相関係数よりも値が小さくなってしまう事がある。データ本来の分布の中の一部分の偏ったデータだけを用いて相関を算出することにより、相関係数の値に影響を及ぼす場合を選抜効果という。

選抜効果を考慮して相関関係を推定する方法に完全情報最尤推定法(FIML)というのがある。

FIMLの理論の説明や実装してみた系の記事は既にあるが、実際にどの程度使えるのかを検証した日本語の記事が見当たらなかったので、勉強を兼ねて実装した。

対象読者

• FIML #とは
• FIMLという名前を聞いたけど、何に使えるのかわからない人

実装のコード部分は無視してOK
FIMLのお気持ちと使い方を何となく理解して、興味が湧けば(村山、2011)を読みに行くといいと思う。
http://koumurayama.com/koujapanese/missing_data.pdf

• FIMLの理論は何となくわかったが、実装の仕方がわからない人

お気持ちはサラッと読み流して大丈夫、コード等が参考になれば嬉しい。

この記事を読んだらわかる事

• FIMLを使うと何が嬉しいのか
• FIMLのざっくりとした理解
• FIMLの実装例(Python + numpy + scipy)
• FIMLを用いた仮説検証の例

この記事を読んでも(絶対に)わからない事

• FIMLのちゃんとした理論
  

FIMLを使うと何が嬉しいのか

観測データを用いて複数の変数の間の関係を調べたい時、観測データが発生した状況を考慮しないと、誤った解析結果を導いてしまう場合がある。

例えば、ある高校の入学試験の成績($X$)と高校入学後の学年末試験の成績($Y$)の相関を調べたいとする。入学試験の点数$x$が合格点$c$を下回った生徒は、仮に入学して学年末試験を受けていたら何点を取ったか、というデータ$y$が観測出来ない。なので、手元のデータは、受験したものの不合格になってしまった生徒の$Y$のデータが欠測している。

$Y$の欠測は完全にランダムに発生しているのではなく、

• $X$が合格点を超えた(x >= c)→$Y$は観測
• $X$が合格点を下回った(x < c)→$Y$は欠測
  

となっている。このように、全てが観測出来ている変数($X$)によって欠測するかどうかが決まるタイプの欠測をMARと言う。

$Y$が欠測するかどうかは完全にランダム(MCAR)であれば、$Y$が観測出来ていないデータを捨てて解析しても問題ないが、MARを仮定したデータに対して$Y$が欠測したデータを捨てて解析した場合、真の解析結果とは違う解析結果になってしまう場合がある。例えば、本当は$X$と$Y$に相関があるのに、$Y$が欠測した$X$のデータを捨てると相関が弱まったり見えなくなってしまう事がある(選抜効果)。

FIMLは、欠測がMARだった場合に、$Y$が観測出来たデータのみを用いて解析をするのではなく、$Y$が観測出来なかったデータも用いて解析する。
これにより、全てのデータが観測出来ていた場合の解析結果により近い解析を行うことが出来る。

上の例では、合格者iのデータ$(x_i, y_i)$のみを用いて解析をするのではなく、不合格者jのデータ$(x_j,)$も用いて解析する事で、不合格者jが学年末試験を受けていた場合の解析結果により近い解析を行うことが出来る。

で、FIMLは何なのか

理論的な説明は
欠損データ分析（missing data analysis）-完全情報最尤推定法と多重代入法-
http://koumurayama.com/koujapanese/missing_data.pdf
等に任せるとして、この記事ではFIMLのお気持ちを掴む。

ちゃんと理論を理解したければ、得体のしれないブログを読むよりちゃんとした教授が書いたPDFとか本を読んだ方が絶対理解が速い。
ただ、初手でガチな本を読んだりするのはとても骨が折れるし、興味を維持するのも大変なので、過度な数式は使わずに「お気持ち」をサクッと説明する。

最尤推定のお気持ち

FIMLは最尤推定に基づいた考え方になっている。

最尤推定の説明は省くが、ざっくり言うと

• データ$x$が、パラメータ$\sigma$によって決められた確率分布から発生している
• 確率分布の式の形はわかっているが、パラメータ$\sigma$の値がわからない
• 手元に観測されたデータ$X = (x_1, x_2, ..., x_n)$がある

としたときに、データ$X$を発生させる確率が最も高くなるパラメータ$\sigma$の値を推定する方法。

未知のパラメータで形が決まる確率分布から手元のデータが発生する確率が、「最」も「尤」もらしくなる、高くなるようなパラメータを「推定」する。

具体的には、個々のデータの対数尤度を計算し、データ全体の合計値を最大化する。

FIMLのお気持ち

FIMLのお気持ちを掴む。

• データ$(x, y)$が、複数のパラメータ$\mathbf{\sigma} = (\sigma_1, \sigma_2, ...)$によって決められたある同時確率分布から発生している
• 確率分布の式の形はわかっているが、パラメータ$\mathbf{\sigma}$の値がわからない
• パラメータの中には、$y$と関係のあるパラメータ($y$の平均$\mu_y$など)と関係のないパラメータ($x$の平均$\mu_x$など)がある
• 手元に観測された$x$, $y$のデータがある。データの中には$(x, y)$が揃ったデータがあれば、$y$が欠測したデータ$(x,)$もある

としたときに、

• $(x, y)$が揃っている→全てのパラメータの推定に使用
• $(x,)$しかない、$y$が欠測している→$y$と関係のないパラメータの推定に使用

する事で、手元のデータ全てを用いてパラメータを推定する方法。

基本的には最尤推定と一緒。
全てが観測されたデータは通常の対数尤度、一部の変数が欠測した場合は欠測してない変数のみに関する確率分布を元の確率分布から抜き出してデータの対数尤度を計算し、データ全体の対数尤度の合計値が最大化するようなパラメータを探索する。
欠測データを前提にしたモデルの推定方法であって、欠測データの補完方法ではない。

※2019/1/26追記
推定した母集団の確率分布と欠測が起こっていない変数のデータを条件付ければ、欠測したデータを補完する事も出来そうです。

FIMLを実際に使ってみる

実装コードは、以下を参考にしました。
kamadak / fiml-py
https://github.com/kamadak/fiml-py/blob/master/fiml.py

Python+NumpyでFIML（完全情報最尤推定法）試してたら色々勘違いをしていたっぽい話
https://ensekitt.hatenablog.com/entry/2018/03/09/200000#fn-d1cba4eb

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from tqdm import tqdm_notebook as tqdm

np.random.seed(1234)
sns.set_style('darkgrid')

%matplotlib inline

import scipy as sp
from scipy import optimize

問題設定

(村山, 2011)に載っている例に従って、以下のような問題を設定する。
10人の被験者に対して、動機付け、IQテスト、適性検査を行なった。ただし、IQが100以下の被験者は適性検査を受けないとする、つまりIQが100以下の被験者の適性検査の結果は観測されない。
ただし、"test_full"は欠測が無かった場合の結果で、実際の解析でこの列を用いる事は出来ない。

ゴール

IQと適性検査の間に相関があるのかを明らかにしたい。

確認したい事

• FIMLを用いて計算したIQテストと適性検査の相関係数$\rho$が、欠測データを捨てた場合と比べて、どれだけ真の値に近づけるか
• FIMLを用いて計算した適性検査の標本平均が、欠測データを捨てた場合と比べて、どれだけ真の値に近づけるか
• (村山, 2011)と同じ解析結果になっているか

モデル

データが発生する構造を定義する。
今回は、3変数(motivation, IQ, test)の確率変数ベクトル$\mathbf{x}$が平均$\mathbf{\mu}$, 分散$\Sigma$の多変量正規分布から発生するとする。

データ

df = pd.DataFrame({
    "motivation": [3,4,5,2,5,3,2,6,3,3],
    "IQ": [83,85,95,96,103,104,109,112,115,116],
    "test_obs": [np.nan, np.nan, np.nan, np.nan,128, 102, 111,113,117,133],
    "test_full": [93,99,98,103,128,102,111,113,117,133]
})
df = df[["motivation", "IQ", "test_obs", "test_full"]]

print(df)

上から4人分の"test_obs"が欠測している。

欠測は、"IQ"の値が100を下回るかどうかで決まっていることがわかる。
散布図を描かなくても"IQ"と"test_full"は相関ありそうだとわかるが、"test_obs"に関してはIQのわりに良いスコアを出している上から5番目の被験者によって相関が弱められそうなことが読み取れる。

df.mean()

"test_obs"と"test_full"で平均値が異なるのがわかる。
(村山, 2011)と"test_full"の平均値が違うけど、多分こちらが正しい。

前処理

データの標準化

地味に一番重要なところ。
値が大きいと計算途中でどんどん誤差が膨らんでいくので、予め標準化しておく。標準化しないと(村山, 2011)の値を復元出来なかった。相関係数の値はそのままで良いが、平均は最後に元の値に戻す。

data = df[["motivation", "IQ", "test_obs"]].values
mean_array = np.nanmean(data, axis=0)
std_array = np.nanstd(data, axis=0)
# nanを無視して要素ごとに計算が出来る、便利
data_std = (data - mean_array) / std_array

関数にデータを投入する準備

size, dim = data_std.shape

mis_num = np.isnan(data_std).sum()
mis_array = data_std[:mis_num, :-1]
obs_array = data_std[mis_num:, :]

obs_col = np.array([True for _ in range(dim)])
mis_col = np.array([True for _ in range(dim-1)] + [False])

data_blocks = [(mis_col, mis_array), (obs_col, obs_array)]

FIMLの実行

FIMLは、観測データの尤度を最大化するパラメータを求める。
観測データの対数尤度関数と初期パラメータを定義し、scipy.optimizeで最適化問題を解く。

初期パラメータの設定

mean0 = np.zeros(dim)
cov0 = np.eye(dim)
params0 = np.empty(dim + (dim * (dim+1)) // 2)
params0[:dim] = mean0
params0[dim:] = cov0[np.tril_indices(dim)]

観測データ全体の対数尤度関数

def _obj_func(params, dim, data_blocks):
    mean = params[0:dim] # paramsから平均に関するパラメータをslice
    cov = np.empty((dim, dim)) # 分散共分散行列を定義
    ii, jj = np.tril_indices(dim) # 下三角行列の0じゃない部分を取り出す
    cov[ii, jj] = params[dim:] # 分散共分散行列は対称行列なので、
    cov[jj, ii] = params[dim:] # 下三角と上三角に同じ分散共分散を代入する

    # 分散共分散行列が対称な半正定値行列になっているかを判別する例外処理
    # 半正定値行列じゃないと分散共分散行列に使えない
    if (np.linalg.eigvalsh(cov) < 0).any():
        return np.inf
    objval = 0.0 # 対数尤度をここに加算していく 
    for obs, obs_data in data_blocks:
        obs_mean = mean[obs]
        obs_cov = cov[obs][:, obs]
        objval += _log_likelihood_composed(obs_data, obs_mean, obs_cov)
    return -objval

def _log_likelihood_composed(x, mean, cov):
    # 多変量正規分布の対数尤度を書き下した関数
    # scipy.stats.multivariate_normal.logpdfを使っても良いが、こちらの方が2倍速かった
    xshift = x - mean
    size = x.shape[0]
    t1 = x.shape[-1] * np.log(2*np.pi) # t1は無くても問題なさそう
    sign, logdet = np.linalg.slogdet(cov)
    t2 = logdet
    t3 = -0.5 * xshift.dot(np.linalg.inv(cov)) * xshift
    return (-0.5 * size * t1) + (-0.5 * size * t2) + t3.sum()

result = optimize.fmin_slsqp(
    _obj_func, params0, args=(dim, data_blocks), disp=True, iter=1000)

解析結果

mean = result[0:dim]
cov = np.empty((dim, dim))
ii, jj = np.tril_indices(dim)
cov[ii, jj] = result[dim:]
cov[jj, ii] = result[dim:]

IQと適性検査の相関係数

def corr(i, j, cov):
    return cov[i, j] / (cov[i,i]**(1/2) * cov[j,j]**(1/2))

# FIMLを用いて欠測を考慮した場合
print("FIMLで欠測データを考慮した相関係数: {:.3f}".format(corr(1, 2, cov)))

# 欠測無しの適性検査から相関係数を計算した場合
data_full = df[["motivation", "IQ", "test_full"]].values
print("真の相関係数: {:.3f}".format(np.corrcoef(data_full.T)[1,2]))

# 観測データのみから相関係数を計算した場合
data_obs = df.loc[mis_num:, ["motivation", "IQ", "test_obs"]].values
print("欠測データを考慮しない相関係数: {:.3f}".format(np.corrcoef(data_obs.T)[1,2]))

欠測データを考慮しないで相関係数を計算した結果は0.342となり、真の相関係数0.773と比べて相関がかなり過小評価されてしまっている。
一方で、FIMLを用いて相関係数を計算した結果は0.669となり、真値0.773に近いことがわかる。
一般的に、FIMLを用いた相関係数が0.669は正の相関があると言えるが、欠測を考慮しないで計算した相関係数0.342では、正の相関は弱いとしか言えない。

観測データが発生した状況を考慮し、場合によっては欠測が起こっている前提でモデルを構築しなければ、誤った判断を下す危険性がある。

※2019/1/26 追記
FIMLを用いた場合の相関係数は推定した母集団の確率分布の分散共分散行列の値を用いて計算した(多分)母集団の相関係数ですが、それ以外の相関係数は標本から直接計算した標本相関係数です。
厳密には2種類の相関係数を直接比較する事は出来ませんが、ここでは差分が明らかに大きいので、簡便のためそのまま比較します。

平均

# 標準化したmeanを元に戻す
print(mean * std_array + mean_array)

正解値: [動機付け, IQ, 適性検査(欠測無し)] = [101.8, 117.3, 111,7]
Mplus(村山, 2017)の予測値: [3.6, 101.8, 110.9]
ちなみに、適性検査(欠測あり)の平均値は117.3

FIMLを用いることで、元々欠測の起こっていない「動機づけ」や「IQ」が真値と一致するのはもちろん、欠測が起きている変数「適性検査」の平均値も110.9となり、真の平均値111.7に近い値となった。(村山, 2017)によるMplusの解析結果とも値が一致した。

まとめ

• 完全情報最尤推定法は、MAR仮定の欠測データに対し、欠測によるバイアスを取り除いた推定が出来る
• 仮想データを用いて、欠測を考慮しない推定法に比べたFIMLの有用性を確認した

次やれたらやりたいこと

stanでベイズ推定したり、多分切断分布仮定しても同じような事出来そう
データのスケールにどこまで対応出来るかも気になる
多くのデータは多分Xがバイナリ(アプリをダウンロードしたか/してないか、等)だから、それらの対応策も勉強したい
多重代入法とかで埋めた場合との比較検証も面白そう