名義反応モデル

名義反応モデル#

順序のないカテゴリカル変数を名義尺度（nominal scale）という。

名義尺度の反応を扱えるのが 名義反応モデル（nominal response model: NRM） softmax関数の形で多値化する。

名義反応モデル（nominal response model: NRM）

個人 \(i\)、項目 \(j\)、カテゴリ \(k \in \{1, \dots, K_j\}\) に対して、反応確率が次のように定義されるモデル

\[ P(y_{ij} = k) = \operatorname{softmax}(a_{jk} \theta_i + b_{jk}) = \frac{ \exp (a_{jk} \theta_i + b_{jk}) }{ \sum_{l=1}^K \exp (a_{jl} \theta_i + b_{jl}) } \]

名義カテゴリモデル（nominal categories model） という呼ばれ方もする

同値変換#

NRMのモデル

\[ P(y_{ij} = k) = \operatorname{softmax}(a_{jk} \theta_i + b_{jk}) \]

は

  • \(a_{jk}' = -a_{jk}\)

  • \(\theta_i' = -\theta_i\)

と置き換えても

\[ a_{jk}' \theta_i' + b_{jk} = (-a_{jk}) (-\theta_i) + b_{jk} = a_{jk} \theta_i + b_{jk} \]

となるため、推定に使うライブラリの違いなどで推定結果が真逆になることがある

多次元バージョン#

Thissen et al. (2010) の General-Purpose Multidimensional Nominal Model (GPMNM) は名義反応モデルで多次元の特性パラメータ\(\boldsymbol{\theta}\)と識別力\(\mathbf a_{jk}\)を扱えるよう一般化したモデル

General-Purpose Multidimensional Nominal Model (GPMNM)

\[ P(Y_{ij} = k \mid \boldsymbol{\theta}_i) = \frac{ \exp\!\left( \mathbf a_{jk}^\top \boldsymbol{\theta}_i + b_{jk} \right) }{ \sum_{c=1}^{K_j} \exp\!\left( \mathbf a_{jc}^\top \boldsymbol{\theta}_i + b_{jc} \right) } \]

  • \(\boldsymbol{\theta}_i \in \mathbb{R}^D\)：個人 \(i\) の 多次元潜在特性（能力ベクトル）

  • \(\mathbf a_{jk} \in \mathbb{R}^D\)：項目 \(j\) のカテゴリ \(k\) に固有な識別力ベクトル

  • \(b_{jk} \in \mathbb{R}\)：項目 \(j\) のカテゴリ \(k\) に固有な困難度

  • \(K_j\)：項目 \(j\) のカテゴリ数

実装例#

# データを生成
import numpy as np
import pandas as pd

def simulate_nrm(N=1000, J=10, K=4, seed=42):
    rng = np.random.default_rng(seed)
    theta = rng.normal(0, 1, size=N)

    # カテゴリ0を参照カテゴリ (a=0, b=0) とする
    a_free = np.sort(rng.normal(0, 1, size=(J, K - 1)), axis=1)
    b_free = rng.normal(0, 1, size=(J, K - 1))

    a_full = np.zeros((J, K))
    a_full[:, 1:] = a_free
    b_full = np.zeros((J, K))
    b_full[:, 1:] = b_free

    U = np.zeros((N, J), dtype=int)
    for i in range(N):
        for j in range(J):
            logits = a_full[j] * theta[i] + b_full[j]
            logits -= logits.max()
            P = np.exp(logits)
            P /= P.sum()
            U[i, j] = rng.choice(K, p=P)
    return U, theta, a_free, b_free

num_users = 1000
num_items = 20
U, true_theta, true_a, true_b = simulate_nrm(N=num_users, J=num_items)
df = pd.DataFrame(
    U,
    index=[f"user_{i+1}" for i in range(num_users)],
    columns=[f"question_{j+1}" for j in range(num_items)],
)
df.head()
question_1 question_2 question_3 question_4 question_5 question_6 question_7 question_8 question_9 question_10 question_11 question_12 question_13 question_14 question_15 question_16 question_17 question_18 question_19 question_20
user_1 0 0 3 2 0 3 3 0 0 3 3 3 3 0 1 1 1 3 3 2
user_2 2 1 0 1 1 1 0 0 1 2 1 3 1 1 1 1 2 3 2 0
user_3 1 2 2 3 0 3 3 2 0 0 3 3 1 2 2 3 3 3 3 2
user_4 2 1 3 0 0 3 3 2 0 3 3 3 3 0 0 1 1 3 3 0
user_5 1 1 0 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 1 2 
# 縦持ちへ変換
df_long = pd.melt(
    df.reset_index(),
    id_vars="index",
    var_name="item",
    value_name="response",
).rename(columns={"index": "user"}).astype({"user": "category", "item": "category"})

df_long.head()
user item response
0 user_1 question_1 0
1 user_2 question_1 2
2 user_3 question_1 1
3 user_4 question_1 2
4 user_5 question_1 1 
import pymc as pm
import pytensor.tensor as pt

user_idx = df_long["user"].cat.codes.to_numpy()
users = df_long["user"].cat.categories.to_numpy()
item_idx = df_long["item"].cat.codes.to_numpy()
items = df_long["item"].cat.categories.to_numpy()
responses = df_long["response"].to_numpy().astype("int64")

K = int(responses.max() + 1)
n_free = K - 1
n_items = len(items)

coords = {
    "user": users,
    "item": items,
    "category_free": np.arange(n_free),
    "obs_id": np.arange(len(df_long)),
}

with pm.Model(coords=coords) as model:
    user_idx_ = pm.Data("user_idx", user_idx, dims="obs_id")
    item_idx_ = pm.Data("item_idx", item_idx, dims="obs_id")

    theta = pm.Normal("theta", 0.0, 1.0, dims="user")
    a_free = pm.Normal("a_free", 0.0, 1.0, dims=("item", "category_free"))
    b_free = pm.Normal("b_free", 0.0, 2.0, dims=("item", "category_free"))

    # 参照カテゴリ（0列目）にゼロを付与
    a_full = pt.concatenate([pt.zeros((n_items, 1)), a_free], axis=1)
    b_full = pt.concatenate([pt.zeros((n_items, 1)), b_free], axis=1)

    # logits[n, k] = a_{j,k} * theta_i + b_{j,k}
    logits = a_full[item_idx_] * theta[user_idx_][:, None] + b_full[item_idx_]
    p = pm.math.softmax(logits, axis=1)

    pm.Categorical("obs", p=p, observed=responses, dims="obs_id")

推定#

%%time
with model:
    idata = pm.sample(random_seed=0, draws=1000)

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (2 chains in 2 jobs)

NUTS: [theta, a_free, b_free]



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File <timed exec>:2

File ~/work/notes/notes/.venv/lib/python3.10/site-packages/pymc/sampling/mcmc.py:957, in sample(draws, tune, chains, cores, random_seed, progressbar, progressbar_theme, step, var_names, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, blas_cores, model, compile_kwargs, **kwargs)
    953 t_sampling = time.time() - t_start
    955 # Packaging, validating and returning the result was extracted
    956 # into a function to make it easier to test and refactor.
--> 957 return _sample_return(
    958     run=run,
    959     traces=trace if isinstance(trace, ZarrTrace) else traces,
    960     tune=tune,
    961     t_sampling=t_sampling,
    962     discard_tuned_samples=discard_tuned_samples,
    963     compute_convergence_checks=compute_convergence_checks,
    964     return_inferencedata=return_inferencedata,
    965     keep_warning_stat=keep_warning_stat,
    966     idata_kwargs=idata_kwargs or {},
    967     model=model,
    968 )

File ~/work/notes/notes/.venv/lib/python3.10/site-packages/pymc/sampling/mcmc.py:1042, in _sample_return(run, traces, tune, t_sampling, discard_tuned_samples, compute_convergence_checks, return_inferencedata, keep_warning_stat, idata_kwargs, model)
   1040 # Pick and slice chains to keep the maximum number of samples
   1041 if discard_tuned_samples:
-> 1042     traces, length = _choose_chains(traces, tune)
   1043 else:
   1044     traces, length = _choose_chains(traces, 0)

File ~/work/notes/notes/.venv/lib/python3.10/site-packages/pymc/backends/base.py:624, in _choose_chains(traces, tune)
    622 lengths = [max(0, len(trace) - tune) for trace in traces]
    623 if not sum(lengths):
--> 624     raise ValueError("Not enough samples to build a trace.")
    626 idxs = np.argsort(lengths)
    627 l_sort = np.array(lengths)[idxs]

ValueError: Not enough samples to build a trace.













EAP推定量#






post_mean = idata.posterior.mean(dim=["chain", "draw"])

import matplotlib.pyplot as plt

fig, axes = plt.subplots(figsize=[12, 4], ncols=3)

ax = axes[0]
ax.scatter(true_theta, post_mean["theta"])
ax.plot(true_theta, true_theta, color="gray")
_ = ax.set(xlabel="true_theta", ylabel="theta_hat")

ax = axes[1]
ax.scatter(true_a.flatten(), post_mean["a_free"].to_numpy().flatten())
ax.plot(true_a.flatten(), true_a.flatten(), color="gray")
_ = ax.set(xlabel="true_a", ylabel="a_hat")

ax = axes[2]
ax.scatter(true_b.flatten(), post_mean["b_free"].to_numpy().flatten())
ax.plot(true_b.flatten(), true_b.flatten(), color="gray")
_ = ax.set(xlabel="true_b", ylabel="b_hat")















参考文献#


    

  • Lipovetsky, S. (2021). Handbook of Item Response Theory, Volume 1, Models. Technometrics, 63(3), 428-431.
    • 







    
    

                

              

              
              
              
              
                
              
            

            
            
              
                



  

  

     Contents