Stan で多分割 多相 Rasch モデルを推定する

~人材採用の最適化を事例として ~

2018年12月20日
株式会社 知能情報システム 井上暁光

はじめに

この記事では Stan を用いた、「polytomous (多分割)」 で 「multi-facet (多相)」 な Rasch モデルの推定方法を解説します。

本稿は以下のように構成されます。

• 本分析によって解決したい課題を解説します。
• Rasch モデルを解説します。
• Stan による実装例を示します。
• 推定例を示します。

解決したい課題

人材管理関連の仮想的な例を考えます。ある企業では採用において、複数の面接担当者に応募者の潜在的な能力を評価してもらい、その評価の素点を合計することで応募者の採否を決めてきました。しかし「なぜあの人を不採用にしたのか」「どういう経緯であの人が採用になったのか」という批判が多くあり、応募者の潜在的な能力を正しく推定できていないと考えられます。この業界においては、応募者の潜在的な能力は応募者全体で標準正規分布に従っていると考えられ、採用後に 1 標準偏差・1人あたり年間 100 万円程度の売上高の違いにつながると考えられています。

調査の結果、現状の方法で潜在的な能力を正しく推定できていない原因は、主に以下の 3 つにあることが分かりました。

1. 採否の境界近辺の応募者の能力の違いが、素点合計の違いに表れていない (A さんと B さんは印象がかなり違うのになぜ素点合計が同じなのか ?)
2. 面接担当者により評価の厳しさが異なるがその違いが反映されていない (X さんに当たると全員不採用、Y さんに当たると全員採用してないか ?)
3. 素点の間隔の意味の違い ( 1 と 2 の違いは 2 と 3 の違いと意味合いが違う) が反映されていない (明らかに問題のある側面があるとき 1 にするのに、なぜこの人は採用されたのか ?)

そこで、上記の原因を解消し、素点合計に代わる指標を開発することになりました。目的は、これまでの素点合計に代わる指標を用いた採用により、潜在的な能力が高い応募者を優先して採用できるようにすることです。目標は、今期採用メンバーによる売上高が、前期採用メンバーによる売上高より 50 万円増とすることです。この企業は年間 10 人を採用しており、この目標を一人当たりに換算すると、50 (万円) / 100 (万円) / 10 (人) = 0.05 となり、潜在的な能力が平均で 0.05 標準偏差分高い人材を採用できるようにすることが、この分析の直接の目標となります。

Polytomous Multi-facet Rasch モデル

基本的なアプローチ

• 原因 1 を Rasch モデルの項目難易度を推定することで解決します。
• 原因 2 はモデルの多相化、つまり面接担当者の評価の厳しさを推定することで解決します。
• 原因 3 はモデルの多分割化、つまり水準間の違いを推定することで解決します。

モデルの定義

Polytomous Multi-facet Rasch モデルは以下の式で表されます。

• $N$ : 対象の数です。
• $I$ : 項目の数です。
• $J$ : 評価者の数です。
• $h$ : 評点の水準の数です。評点は 1 origin の整数で表されます。
• $X_{n,i,j}$ : データを表す整数です。$n$ 番目の対象の $i$ 番目の項目の $j$ 番目の評価者による評価値です。$N \times I \times J$ のサイズをもつ整数テンソルです。
• $P_{n,i,j,k}$ : データが観測される確率です。$n$ 番目の対象の $i$ 番目の項目の $j$ 番目の評価者による評価値が、$k$ である確率を表します。$N \times I \times J \times h$ のサイズをもつ実数値テンソルです。(全ての $n, i, j$ について $\sum_{k=1}^{h} P_{n,i,j,k} = 1$ が成り立ちます)
• $B_n$ : 対象の能力です。$n$ 番目の対象の能力を表します。要素数 $N$ の実数値ベクトルです。
• $D_{i}$ : 項目の難しさです。$i$ 番目の項目の難しさを表します。要素数 $I$ の実数値ベクトルです。
• $C_j$ : 評価者の厳しさです。$j$ 番目の評価者の採点の厳しさを表します。要素数 $J$ の実数値ベクトルです。
• $F_k$ : 評点が水準間でどの程度異なるかを表します。$k-1$ 番目の水準と $k$ 番目の水準の差を表します。要素数 $h$ の実数値ベクトルです。$F_1 = 0$ と仮定します。

ここで、対象、項目、評価者を所与として、水準 $k-1$ が選ばれる確率に対する水準 $k$ が推定される確率の比の対数は以下で表されるという仮定をおくモデルが Rasch モデルです。

尤度で推定するために $P$ を陽に表した式が以下です。

Stan による実装

Stan による実装では、モデルが冗長なパラメータを含むことに気を付けてください。冗長であると、尤度を最大化する解が 1 つでなく無数に存在するので、パラメータの事後分布が正しく推定できません。これを防ぐため、本稿では、$B$, $D$, $C$ それぞれの和を 0 とする制約を与えます。

なお、実用上は MMLE 等の推定法を用いる方が高速となりますが、MCMC 法で推定する方法は、各パラメータの事前分布を与えたり、共変量を加えるなどのモデル拡張も容易であるといったメリットがあります。

以下に Stan による実装を示します。

​x
data {
    int<lower=1> N;
    int<lower=1> I;
    int<lower=1> J;
    int<lower=2> h;
    int<lower=1> X[N, I, J];
}
​
parameters {
    // B, D, C をそれぞれ合計を 1.0 に制約を掛けたいので、
    // それぞれの最後の要素を除いたベクトルを推定対象のパラメータとします。
    real B_without_last[N - 1];
    real D_without_last[I - 1];
    real C_without_last[J - 1];
​
    // F[1] は 1 で固定であるため F[2 から h] を推定対象とします。
    // transformed parameters セクションで F という
    // ベクトルを作成します。
    real F_without_first[h - 1];
}
​
transformed parameters {
    real B[N];
    real D[I];
    real C[J];
​
    // parameter セクションの F_without_first を扱いやすいように
    // 最初の要素に定数を入れた F を作成します。
    real F[h];
​
    // P 自体は推定対象のパラメータではないので lower, upper の設定や
    // simplex として扱うことはしません。
    real P[N, I, J, h];
​
    // F の 1 番目からの累積和です。
    // たとえば cumulative_sum_F[3] == F[1] + F[2] + F[3] です。
    real cumulative_sum_F[h];
    
    real denominator;
    
    // B, D, C の計算です。
    // sum(B) == 1 となるように正規化します。
    for (n in 1:(N-1)) {
        B[n] = B_without_last[n];
    }
    B[N] = 1.0 - sum(B[1:(N-1)]);
​
    // sum(D) == 1 となるように正規化します。    
    for (i in 1:(I-1)) {
        D[i] = D_without_last[i];
    }
    D[I] = 1.0 - sum(D[1:(I-1)]);
​
    // sum(C) == 1 となるように正規化します。    
    for (j in 1:(J-1)) {
        C[j] = C_without_last[j];
    }
    C[J] = 1.0 - sum(C[1:(J-1)]);
    
    // F の計算です。
    F[1] = 0.0;
    for (k in 2:h) {
        F[k] = F_without_first[k-1];
    }
    
    // cumulative_sum_F の計算です。
    for (k in 1:h) {
        cumulative_sum_F[k] = 0.0;
        for (m in 1:k) {
            cumulative_sum_F[k] += F[m];
        }
    }
    
    // P の計算です。
    for (n in 1:N) {
        for (i in 1:I) {
            for (j in 1:J) {
                denominator = 0.0;
                for (k in 1:h) {
                    denominator += exp(k * (B[n] - D[i] - C[j]) - 
                                       cumulative_sum_F[k]);
                }
                
                for (k in 1:h) {
                    P[n, i, j, k] = exp(k * (B[n] - D[i] - C[j]) - 
                                        cumulative_sum_F[k]) / denominator;
                }
            }
        }
    }
}
​
model {
    for (n in 1:N) {
        for (i in 1:I) {
            for (j in 1:J) {
               for (k in 1:h) {
                   if (X[n,i,j] == k) {
                       target += log(P[n,i,j,X[n,i,j]]);
                   }
               }
            }
        }
    }
}

分析例

データ セットの生成

課題の説明に適したデータを人工的に生成します。Rasch モデルに基づいてデータを生成しています。今回は理解を優先し、モデルの mis-specification については考えません。

xxxxxxxxxx
import numpy as np
import pystan
import pickle
from hashlib import md5
import matplotlib.pyplot as plt
import seaborn as sns
​
np.random.seed(1)
​
# 候補者の数
N = 100
​
# 項目の数
I = 20
​
# 評価者の数
J = 10
​
# 評点のカテゴリ数
h = 5
​
# データ
X: np.ndarray = np.zeros(shape=(N, I, J), dtype=int)
​
# 候補者の潜在的な能力です。
# 標準正規分布に従うものとします。
B: np.ndarray = np.random.normal(loc=0, scale=1., size=N)
​
# 項目の難しさ
# ほどほどに高い難易度の問題が存在しないようにします。
D: np.ndarray = np.hstack((
    np.random.normal(loc=-1., scale=1., size=15),
    np.random.normal(loc=2, scale=0.1, size=5),
))
​
# 評価者の厳しさです。
# 標準正規分布に従うものとします。
C: np.ndarray = np.random.normal(loc=0., scale=1., size=J)
​
# 評点が水準間でどの程度異なるかを表します。
# 0 と 1 の間の差はかなり大きいものとします。
F: np.ndarray = np.array([0.0, 2.0, 0.5, 0.5, 0.5])
​
# 確率を定義します。
P: np.ndarray = np.zeros((N, I, J, h))
for n in range(N):
    for i in range(I):
        for j in range(J):
            denominator = 0.0
            for k in range(h):
                denominator += np.exp(k * (B[n] - D[i] - C[j] - np.sum(F[:k]) ))
            for k in range(h):
                P[n, i, j, k] = np.exp(k * (B[n] - D[i] - C[j] - np.sum(F[:k]) )) / denominator
​
# 確率に従いデータを生成します。
for n in range(N):
    for i in range(I):
        for j in range(J):
            k = np.random.choice(a=h, size=1, p=P[n,i,j])[0] + 1
            X[n, i, j] = k
​
# 担当表を整理します。
X_part = np.zeros(shape=(N, I, J), dtype=int)
n_assingments = 2
for n in range(N):
    # 評価者をランダムに選びます。
    js = np.random.choice(a=J, size=n_assingments)
    for i in range(I):
        for j in js:
            X_part[n, i, j] = X[n, i, j]

素点合計の算出

既存の方法で素点を合計します。

xxxxxxxxxx
# 候補者別の素点合計を示します。
raw_score_sum = np.sum(X_part, axis=(1, 2))

Rasch モデル推定値の算出

Rasch モデルの最尤推定値を求めます。

xxxxxxxxxx
data_stan = {
    'N': N,
    'I': I,
    'J': J,
    'h': h,
    'X': X,
}
​
# モデルの再コンパイルを防ぐコードです。
# 以下のリンクから取得し、文字コードの指定を ascii から utf-8 に変更しました。
# https://pystan.readthedocs.io/en/latest/avoiding_recompilation.html
def StanModel_cache(model_code, model_name=None, **kwargs):
    """Use just as you would `stan`"""
    code_hash = md5(model_code.encode('utf-8')).hexdigest()
    if model_name is None:
        cache_fn = 'cached-model-{}.pkl'.format(code_hash)
    else:
        cache_fn = 'cached-{}-{}.pkl'.format(model_name, code_hash)
    try:
        sm = pickle.load(open(cache_fn, 'rb'))
    except:
        sm = pystan.StanModel(model_code=model_code)
        with open(cache_fn, 'wb') as f:
            pickle.dump(sm, f)
    else:
        print("Using cached StanModel")
    return sm
​
model = StanModel_cache(model_code=open('multi_facet_polytomous_rasch.stan', 'r').read())
fit = model.sampling(data=data_stan, chains=1, iter=2000, init=0)
print(fit)
​
rasch_estimates = np.mean(fit['B'], axis=0)
print(np.corrcoef(rasch_estimates, B))

採用シミュレーション

100 人の応募者から、素点合計の上位 10 人を抽出した場合と、Rasch モデルのパラメータ推定値をもとに上位 10 人を抽出した場合の、抽出された 10人の潜在的能力を比較します。

Rasch モデルの推定値では平均値が高い 10 人を抽出できていることが分かります。

潜在的能力を 0.05 以上高めるという目標をこのデータで実現できたことが分かります。(誤差は評価していません)

xxxxxxxxxx
# 10人採用することにします。
n_recruit = 10
​
# 単純な素点合計の上位の応募者を抽出します。
index_descending_raw_score_sum = np.argsort(raw_score_sum)[::-1][:n_recruit]
# 単純な素点合計の上位の応募者の潜在的な能力の平均です。
print(np.mean(B[index_descending_raw_score_sum]))
​
# 事後分布の平均が上位の応募者を抽出します。
index_descending_rasch_estimates = np.argsort(rasch_estimates)[::-1][:n_recruit]
# 事後分布の平均が上位の応募者の潜在的な能力の平均です。
print(np.mean(B[index_descending_rasch_estimates]))
​
plt.figure()
plt.boxplot(
    [B, B[index_descending_raw_score_sum], B[index_descending_rasch_estimates]],
    labels=['all candidate', f'top {n_recruit}\n(sum of raw score)', f'top {n_recruit}\n(rasch estimates)'])
plt.title("applicant's potential abilities")
plt.savefig('boxplot.png')
plt.close()

免責事項

本稿の情報については十分な注意を払っておりますが、その内容の正確性等に対して一切保証するものではありません。本稿の利用で生じたいかなる結果についても、当社は一切責任を追わないものとします。

本稿の事例は、評価軸の数、種々の分布、採用プロセス、尺度を構成する項目群の設定などにおいて、単純化するための仮定を置いています。実務への応用に当たっては、これらの仮定をより詳細に検討してください。

参考文献

1. Linacre, J. M. (1999). Understanding Rasch measurement: estimation methods for Rasch measures. Journal of outcome measurement, 3, 382-405.
2. Mair & Hatzinger (2010), Fitting the Rasch Model in R - The eRm Package, http://statmath.wu-wien.ac.at/people/hatz/psychometrics/10w/RM_handouts_4.pdf

© Copyright 2018, 株式会社 知能情報システム