43. 項目反応理論 (IRT) 応用

• このセクションで使っている packages

library(DT)
library(glpkAPI)
library(irt)
library(irtoys)
library(ltm)
library(plink)
library(plyr)
library(psych)
library(reactable)
library(tidyverse)

7. 実際のテストデータを使った分析

7.0 データの準備

• 220人が受験した試験結果 (irt_was.csv)
  
• 含まれている項目は 36 個 (Q1,…, Q36)
  
• 回答結果は 1 が正答、0 が誤答

変数名	詳細
ID	受験者のID
Q1 〜 Q36	1番目から36番目の回答結果 (0 or 1)

• データを読み込み df1 と名前を付ける

df1 <- read_csv("data/irt_was.csv")

DT::datatable(df1)

• Q1からQ36までの総合点 total を作成して df2 として保存

df2 <- df1 |> 
  dplyr::mutate(total = rowSums(dplyr::across(Q1:Q36), 
    na.rm = TRUE))     

DT::datatable(df2)

7.1 正答率の計算: colMeans()

各項目の難易度を平均正答率で把握

• Q1（＝2行目）から Q36（＝37行目）までのスコアの平均値を計算し、crr1 と名前を付ける

crr2 <- colMeans(x = df2[2:37],
  na.rm = TRUE)

• 計算された値 (crr2) では使い勝手が悪いのでデータフレームに変換して df_crr2 と名前を付ける

df_crr2 <- data.frame(      # データフレーム名を指定（ここでは df_crr2 と指定）
  item = names(crr2),       # 変数名を指定（ここでは item と指定）
  seikai = as.numeric(crr2) # 変数名を指定（ここでは seikai と指定）
)

• df_crr2 の中身を確認

df_crr2

   item    seikai
1    Q1 0.7442922
2    Q2 0.8818182
3    Q3 0.8954545
4    Q4 0.3744076
5    Q5 0.4139535
6    Q6 0.5849057
7    Q7 0.8394495
8    Q8 0.5616438
9    Q9 0.5934579
10  Q10 0.4657534
11  Q11 0.2442396
12  Q12 0.5642202
13  Q13 0.5990783
14  Q14 0.2227273
15  Q15 0.9090909
16  Q16 0.8181818
17  Q17 0.3594470
18  Q18 0.7671233
19  Q19 0.5909091
20  Q20 0.6318182
21  Q21 0.7808219
22  Q22 0.7181818
23  Q23 0.1651376
24  Q24 0.2139535
25  Q25 0.3087558
26  Q26 0.3301887
27  Q27 0.9363636
28  Q28 0.9727273
29  Q29 0.9409091
30  Q30 0.3802817
31  Q31 0.7488584
32  Q32 0.3607306
33  Q33 0.8447489
34  Q34 0.9589041
35  Q35 0.8479263
36  Q36 0.8669725

• 棒グラフで可視化する
  
• seikai の値が大きい順に並べ替えて表示させる
  
• seikai の値が大きい順に因子の順序を指定

df_crr2$item <- factor(df_crr2$item, 
  levels = df_crr2$item[order(df_crr2$seikai)])

ggplot(df_crr2, aes(x = seikai, y = item)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  geom_text(aes(label = round(seikai, 2)),  # 小数第2位で丸める
            hjust = -0.2, size = 3) +        # 棒の内側に表示
  labs(
    title = "各項目の正答率",
    x = "項目",
    y = "正答率"
  ) +
  theme_minimal() +
  theme_bw(base_family = "HiraKakuProN-W3") # 文字化け対策  

• Q28の正解率が最も高く(97%)、Q23の正解率が最も低い (17%)

正答率の計算のポイント ・極端に正答率の高い/低い項目があるかどうか
・極端に高い/低い項目がある場合　→　問題あり
・極端に高い/低い項目がない場合　→　問題なし
・ここでは次の項目をチェックする
→　極端に高い項目（90%以上の正答率）・・・Q28〜Q3
→　極端に低い項目（20%以下の正答率）・・・Q23

→　次の分析に移る

7.2 I-T相関の計算: cor()

項目が全体スコアと一致するか（弁別力の目安）を確認

• cor()関数を使って、素点 (Q1 〜 Q36) と合計点 total との相関を計算
  
• Q1 〜 Q36の位置（＝列番号）を確認

dim(df2)

[1] 220  38

DT::datatable(df2)

• 受験生は 220名
  

• Q1・・・2列目
  

• Q36・・・37列目
  

• total・・・38列目
  

• 次の二つの相関を計算する

• Q1（＝2番目の項目）〜 Q36（＝36番目の項目）
  

• total（＝37番目の項目）

it2 <- cor(x = df2[, 2:37], 
  y = df2[, 38],
  use = "pairwise.complete.obs")

it2

         total
Q1  0.29336063
Q2  0.33562403
Q3  0.28121742
Q4  0.51560591
Q5  0.19104289
Q6  0.26743612
Q7  0.28791771
Q8  0.15180953
Q9  0.42685694
Q10 0.37875091
Q11 0.14832979
Q12 0.34877935
Q13 0.22713900
Q14 0.22160121
Q15 0.42822127
Q16 0.29070948
Q17 0.34604474
Q18 0.42356302
Q19 0.40954602
Q20 0.57798397
Q21 0.40066530
Q22 0.41258493
Q23 0.01552610
Q24 0.09734145
Q25 0.02357517
Q26 0.22185973
Q27 0.29641284
Q28 0.29901267
Q29 0.17982935
Q30 0.26133634
Q31 0.42912820
Q32 0.27226490
Q33 0.35369579
Q34 0.18603293
Q35 0.28148011
Q36 0.32466814

• ここで得られた結果 it2 は「行名付きの1列行列（matrix）」
  →　使い勝手が悪いので、この matrix をデータフレームに変換して、行名を項目名の列として追加する

# 行列をデータフレームに変換
df_it2 <- as.data.frame(it2)

# 行名を項目名として列に追加
df_it2$item <- rownames(df_it2)

# 列名をわかりやすく変更（オプション）
colnames(df_it2) <- c("correlation", "item")

• 相関係数が低い順位並べ変えて表示させてみる
  
• correlation の値順に因子の順序を指定

df_it2$item <- factor(df_it2$item, 
  levels = df_it2$item[order(df_it2$correlation)])

ggplot(df_it2, aes(x = correlation, y = item)) +
  geom_bar(stat = "identity", fill = "orange") +
  geom_text(aes(label = round(correlation, 3)), 
    hjust = -0.2, size = 3) +
  xlim(0, 0.65) +
  labs(
    title = "項目－合計相関（item-total correlation）",
    x = "項目",
    y = "相関係数"
  ) +
  theme_minimal() +
  theme_bw(base_family = "HiraKakuProN-W3") # 文字化け対策  

• 素点 (Q1 〜 Q36) と各項目得点 (total) との相関は 0.016と0.578の間

・各項目への反応 (Q1 〜 Q36) と合計点 (total) との間にI-T 相関の度合いは妥当？

IT相関値	評価	項目の扱い	該当項目
〜 0.2	極めて低い（要注意）	除外を検討する	Q5, Q25, Q24, Q11, Q8, Q29, Q34, Q5
0.2〜0.3	やや低い	内容によって再検討
0.3〜0.4	妥当なレベル	保留・文脈による判断
0.4以上	良好（望ましい）	採用して問題なし	Q20, Q4, Q31, Q15, Q9, Q18, Q22, Q19

→　ここでは0.2 以下の相関の項目が8つ　→　除外を検討

7.3 一次元性の検討: psych::fa.parellel()

全項目が1つの潜在特性で測れているかを確認

• 因子分析で利用される方法を用いて検討することが多い
  
• ここでは psych::fa.parallel() による平行分析スクリー・プロットを描いてみる（Horn, 1965)

1番目の固有値が大きく、2番目以下が急激に下がっていれば.

→ 一次元性が強い

• パッケージ読み込み

library(psych)

• データの読み取り

df2 <- read.csv("data/irt_was.csv")

• データの加工

item_df2 <- df2[, -1, -38] # 1列目と38列目は分析に含めない  

# 一次元性の検討（主成分＋平行分析）
fa.parallel(item_df2, 
  fa = "pc", 
  n.iter = 100) # 100回のシミュレーションを指定

Parallel analysis suggests that the number of factors =  NA  and the number of components =  7 

x軸とy軸が意味すること

軸	示している内容	解釈のポイント
X軸	因子番号（主成分番号）	因子1〜因子36
Y軸	固有値（説明される分散量）	1より大きい → 有意味な因子の可能性あり

→　固有値（Eigenvalue）
=> 主成分分析や因子分析において「その因子が説明するデータ全体の分散の量」を示す数値
→　各因子がデータの分散をどれだけ説明しているかを表す指標
→　固有値の値が大きいほど、その因子は多くの情報（＝ばらつき、構造）をキャッチしており
→　その因子で、たくさんの項目のスコアをまとめて説明できるということ

ポイント・・・注目すべきは青線 →　実データの固有値

• いくつの因子で青線が赤線の上にあるかどうか
  
• 因子1〜因子7の青線（実データ）だけが、赤線（ランダムデータ）より遙かに上にある
  
• 因子8以降は、青線は赤線を下回っている

分析結果 ・実データの固有値が12項目でランダムを超える
・明確な多次元性あり（一次元性モデルは使えない）

3つのデータの適合性（正答率、I-T相関、1次元性）をクリアできず

対処法

対処法	説明
① 多次元IRT（MIRT）モデルの推定	mirt(…, 2) や …, 3 などで2次元・3次元モデルを検討
② 因子分析（EFA）で項目を分類	似た因子に属する項目をまとめて、部分テストとしてIRTを適用
③ 次元を1つに統一するよう項目を精選	異なる内容の項目を削除・再構成して一次元性を確保

→　ここでは「③次元を1つに統一するよう項目を精選」して推定してみる

• 多次元性が検出された場合に、IRTで信頼できる1次元モデルを使うための「項目の精選」は、実務でもよく使われる重要な対応

因子分析と因子負荷に基づく項目選抜

セクション	内容
1. パッケージの読み込み	psych, mirt を使用
2. データ保存(item_data)	irt_was.csv から1列と38列を除いた項目データ
3. 因子分析（EFA）	fa(…, nfactors=2, rotate=“promax”) で因子負荷を算出
4. 項目の精選	因子1に強く、因子2に弱く関係する項目のみを抽出（例：ML1 > 0.4, ML2 < 0.3）
5. 条件分岐処理	項目が1つも抽出されなかったときはメッセージを表示し、それ以上の分析を行わない
6. 一次元性の再確認	fa.parallel() を再度実行し、精選項目の構造を確認
7. IRT（2PL）モデルの再推定	mirt(…, 1) を使ってモデルパラメータを推定

パッケージ読み込み

library(psych)
library(mirt)

データ読み込みと保存

df1 <- read_csv("data/irt_was.csv")
item_data <- df1[, -1, -38]

精選項目によるデータ抽出と分析（項目が存在する場合のみ実行）

• 主因子への負荷(ML1 > 0.4) + 他因子への負荷 (ML2 < 0.3)でフィルターをかける

efa_result <- fa(item_data, nfactors = 2, rotate = "promax") # psychパッケージ
num_items <- ncol(item_data)
loadings_df <- as.data.frame(as.matrix(efa_result$loadings)[1:num_items, ])
colnames(loadings_df) <- paste0("ML", seq_len(ncol(loadings_df)))

# 安全に項目を選抜（条件に合う項目がない場合に対応）
selected_items <- tryCatch({
  if (ncol(loadings_df) >= 2) {
    rownames(loadings_df[abs(loadings_df$ML1) > 0.4 & abs(loadings_df$ML2) < 0.3, ])
  } else {
    rownames(loadings_df[abs(loadings_df$ML1) > 0.4, ])
  }
}, error = function(e) character(0))

# 結果を表示
selected_items

[1] "Q27" "Q31" "Q33" "Q35" "Q36"

• 当てはまる項目は 5 つだけ

• 結果を可視化

if (length(selected_items) > 0) {
  item_data_reduced <- item_data[, selected_items]
  
  # 一次元性の再確認
  fa.parallel(item_data_reduced, fa = "pc", n.iter = 100)　# psychパッケージ
  
  # IRTモデルの再推定
  mod_1d <- mirt(item_data_reduced, 1, itemtype = "2PL")　# mirtパッケージ
  summary(mod_1d)
  
} else {
  cat("条件に合致する項目が1つもありません。閾値を緩めて再検討してください。")
}

Parallel analysis suggests that the number of factors =  NA  and the number of components =  1 
Iteration: 1, Log-Lik: -423.603, Max-Change: 0.33466Iteration: 2, Log-Lik: -417.094, Max-Change: 0.24168Iteration: 3, Log-Lik: -414.564, Max-Change: 0.16075Iteration: 4, Log-Lik: -413.265, Max-Change: 0.07673Iteration: 5, Log-Lik: -413.140, Max-Change: 0.05575Iteration: 6, Log-Lik: -413.083, Max-Change: 0.04113Iteration: 7, Log-Lik: -413.036, Max-Change: 0.01829Iteration: 8, Log-Lik: -413.031, Max-Change: 0.01410Iteration: 9, Log-Lik: -413.029, Max-Change: 0.01056Iteration: 10, Log-Lik: -413.025, Max-Change: 0.00493Iteration: 11, Log-Lik: -413.025, Max-Change: 0.00413Iteration: 12, Log-Lik: -413.025, Max-Change: 0.00302Iteration: 13, Log-Lik: -413.024, Max-Change: 0.00126Iteration: 14, Log-Lik: -413.024, Max-Change: 0.00107Iteration: 15, Log-Lik: -413.024, Max-Change: 0.00087Iteration: 16, Log-Lik: -413.024, Max-Change: 0.00033Iteration: 17, Log-Lik: -413.024, Max-Change: 0.00093Iteration: 18, Log-Lik: -413.024, Max-Change: 0.00028Iteration: 19, Log-Lik: -413.024, Max-Change: 0.00020Iteration: 20, Log-Lik: -413.024, Max-Change: 0.00013Iteration: 21, Log-Lik: -413.024, Max-Change: 0.00012Iteration: 22, Log-Lik: -413.024, Max-Change: 0.00046Iteration: 23, Log-Lik: -413.024, Max-Change: 0.00010Iteration: 24, Log-Lik: -413.024, Max-Change: 0.00006
       F1    h2
Q27 0.595 0.354
Q31 0.748 0.560
Q33 0.675 0.456
Q35 0.654 0.428
Q36 0.741 0.549

SS loadings:  2.347 
Proportion Var:  0.469 

Factor correlations: 

   F1
F1  1

分析結果 ・実データ（青瀬）がランダムデータ（赤線）より大きいのは因子1だけ
・そのような因子は因子1だけ
→　一次元性が成立している
→　2PLやRaschなど1因子IRTモデルを使用する前提が妥当
→　安心してIRTモデル（例：ltm() や mirt()）の分析へ進めてOK

• 条件が厳しすぎて項目が0件になることは実際よくある
  
• ここでは標準的なレベルの条件でフィルターに書けたところ 5 項目だけ条件をクリア

条件を緩めて、選定範囲を広げてみる

項目	推奨値（目安）	解釈・理由
ML1（主因子への負荷）	0.40以上 → 最低0.30程度	0.30は「実用的に弱いが意味のある関連性」の下限
ML2（他因子への負荷）	0.30未満 → 最大0.40程度	0.40を超えると交絡項目の可能性 ↑

efa_result <- fa(item_data, nfactors = 2, rotate = "promax")
num_items <- ncol(item_data)
loadings_df <- as.data.frame(as.matrix(efa_result$loadings)[1:num_items, ])
colnames(loadings_df) <- paste0("ML", seq_len(ncol(loadings_df)))

# 安全に項目を選抜（条件に合う項目がない場合に対応）
selected_items <- tryCatch({
  if (ncol(loadings_df) >= 2) {
    rownames(loadings_df[abs(loadings_df$ML1) > 0.3 & abs(loadings_df$ML2) < 0.4, ])
  } else {
    rownames(loadings_df[abs(loadings_df$ML1) > 0.4, ])
  }
}, error = function(e) character(0))

# 結果を表示
selected_items

[1] "Q15" "Q20" "Q27" "Q29" "Q31" "Q33" "Q34" "Q35" "Q36"

• 当てはまる項目が9つに増えた

精選項目によるデータ抽出と分析

if (length(selected_items) > 0) {
  item_data_reduced <- item_data[, selected_items]
  
  # 一次元性の再確認
  fa.parallel(item_data_reduced, fa = "pc", n.iter = 100)
  
  # IRTモデルの再推定
  mod_1d <- mirt(item_data_reduced, 1, itemtype = "2PL")
  summary(mod_1d)
  
} else {
  cat("条件に合致する項目が1つもありません。閾値を緩めて再検討してください。")
}

Parallel analysis suggests that the number of factors =  NA  and the number of components =  1 
Iteration: 1, Log-Lik: -693.607, Max-Change: 0.50905Iteration: 2, Log-Lik: -683.236, Max-Change: 0.24403Iteration: 3, Log-Lik: -680.515, Max-Change: 0.14078Iteration: 4, Log-Lik: -679.575, Max-Change: 0.07502Iteration: 5, Log-Lik: -679.332, Max-Change: 0.04290Iteration: 6, Log-Lik: -679.236, Max-Change: 0.02429Iteration: 7, Log-Lik: -679.192, Max-Change: 0.01491Iteration: 8, Log-Lik: -679.176, Max-Change: 0.00975Iteration: 9, Log-Lik: -679.168, Max-Change: 0.00633Iteration: 10, Log-Lik: -679.162, Max-Change: 0.00247Iteration: 11, Log-Lik: -679.161, Max-Change: 0.00192Iteration: 12, Log-Lik: -679.161, Max-Change: 0.00134Iteration: 13, Log-Lik: -679.161, Max-Change: 0.00033Iteration: 14, Log-Lik: -679.161, Max-Change: 0.00069Iteration: 15, Log-Lik: -679.161, Max-Change: 0.00028Iteration: 16, Log-Lik: -679.161, Max-Change: 0.00012Iteration: 17, Log-Lik: -679.160, Max-Change: 0.00010
       F1    h2
Q15 0.549 0.302
Q20 0.556 0.309
Q27 0.733 0.537
Q29 0.737 0.543
Q31 0.685 0.469
Q33 0.709 0.503
Q34 0.636 0.405
Q35 0.617 0.380
Q36 0.694 0.482

SS loadings:  3.931 
Proportion Var:  0.437 

Factor correlations: 

   F1
F1  1

3つのデータの適合性（正答率、I-T相関、1次元性）をクリア

• 1次元性に関しては、条件を弱めてクリア
  →　項目パラメーター（「識別力a」と「困難度b」）と潜在特性値（学力θ）の推定へ

7.4 項目パラメーターの推定

項目ごとの「識別力a」「困難度b」を数値として推定

・ここでは 2 パラメタ・ロジスティックモデル (2PL: 一般化ロジスティックモデル)を使って分析する

2 パラメタ・ロジスティックモデルの目的 テストに出題した項目パラメーターをもとに、受験者の正誤パターンが最も生じやすい学力θを推定する

• 潜在特性値は平均が０、分散が１の正規分布（標準正規分布）に従うと仮定して推定を行うのが一般的

項目パラメータの推定法

• 学力θは受験者の正誤パターンだけを使って推定する
  
• 項目パラメータも受験者の正誤パターンから推定する

推定対象	推定に必要な情報
学力θ	推定したい受験者の正誤パターン
項目パラメータ	全ての受験者の正誤パターン

・項目パラメータは、その項目を含むテストを受験した受験生の正答・誤答情報をもとに、IRTを使って分析して得られる

• 項目パラメータ（識別力・困難度）と受験者の学力θを「同時に」推定する
  →　「同時」最尤推定法を使う

同時最尤推定法

項目パラメータ推定に関してできることとできないこと

私たちができないこと：

・「識別力が○○、困難度が△△の問題」をあらかじめ作ること
←　試験結果を使って IRT分析しないと項目パラメータを得られないから

私たちができること：

・多くの項目を作成し、それを受験生に受けてもらって IRT 分析する
→　ひとつひとつの項目の検証・検討を積み重ねる

1PLモデルと2PLモデルの違い

特徴	1PLモデル（Raschモデル）	2PLモデル（一般化ロジスティックモデル）
モデル式	\(P（正答）= \frac{1}{1+e^-(\theta-b)}\)	\(P（正答）= \frac{1}{1+e^{-a}(\theta-b)}\)
識別力 a	すべての項目で同じ（固定）	項目ごとに推定
困難度 b	項目ごとに推定	項目ごとに推定
パラメータ数	項目数（b のみ）＋1（a 固定）	項目数 × 2（a と b をそれぞれ推定）
分析対象	能力と困難度 b の関係	能力、困難度 b、識別力 a の関係（より柔軟なモデル）

• 1PLモデル（Raschモデル）では「すべての項目が等しい識別力を持つ」という強い仮定がある
  
• 2PLモデルはその仮定を緩め、実際の項目の識別力の違いを反映する 　　

→　1PLモデルでは「困難度（bパラメータ）」だけが推定される
→　2PLモデルでは「困難度（bパラメータ）」と「識別力（aパラメータ）」が推定され
→ 2PLモデルでは、項目ごとの「能力の区別のしやすさ」が明らかになる

2 パラメタ・ロジスティックモデル (2PL)

• 2PLモデルは、1PLのRaschモデルより柔軟
  →　パラメータ数が増える分だけ推定の安定性はやや下がる
  
• 2PLモデルの数式：

\[P_j(\theta) = \frac{1}{1 + exp[-1.7a_j(\theta-b_j)]}\]

表記	詳細
\(P(\theta)\)	能力θを持つ受験者が、ある項目に正答する確率
\(\theta\)	受験者の能力で、平均0・標準偏差1の正規分布に従うと仮定
\(a\)	識別力（discrimination）パラメータ
\(b\)	困難度（difficulty）パラメータ＝位置パラメータ
\(j\)	項目の番号

• df2から一次元性の条件を満たす9つの項目だけを選んだデータフレームを作る

df_nine <- df2 |> 
  select("Q15", "Q20", "Q27", "Q29", "Q31", "Q33", "Q34", "Q35", "Q36")

2 パラメタ・ロジスティックモデル (2PL)

• ここでは 9 項目の「識別力a」と「困難度b」が計算する

DT::datatable(df_nine)

• ltm パッケージの est()関数で項目パラメーターを推定する

→　ex2 として保存

ex2 <- est(resp = df_nine[, 1:9], # テストデータを指定する引数
  model = "2PL",        # 2PLMを仮定
  engine = "ltm")　# ltmパッケージを利用して項目パラメーターを推定すると指定

ex2

$est
        [,1]       [,2] [,3]
Q15 1.116352 -2.4774085    0
Q20 1.138009 -0.5955641    0
Q27 1.831933 -2.1289528    0
Q29 1.856599 -2.1700494    0
Q31 1.598930 -0.9810200    0
Q33 1.709964 -1.4333128    0
Q34 1.402326 -2.8348121    0
Q35 1.331616 -1.6712500    0
Q36 1.638347 -1.6180664    0

$se
           [,1]      [,2] [,3]
 [1,] 0.3674157 0.6197462    0
 [2,] 0.2972607 0.1851252    0
 [3,] 0.5913987 0.3908171    0
 [4,] 0.5998113 0.3966968    0
 [5,] 0.4023416 0.1880701    0
 [6,] 0.4552229 0.2457418    0
 [7,] 0.5235105 0.7222579    0
 [8,] 0.3608220 0.3304256    0
 [9,] 0.4387965 0.2797887    0

$vcm
$vcm[[1]]
          [,1]      [,2]
[1,] 0.1349943 0.2106149
[2,] 0.2106149 0.3840853

$vcm[[2]]
           [,1]       [,2]
[1,] 0.08836392 0.03069508
[2,] 0.03069508 0.03427134

$vcm[[3]]
          [,1]      [,2]
[1,] 0.3497524 0.2012777
[2,] 0.2012777 0.1527380

$vcm[[4]]
          [,1]      [,2]
[1,] 0.3597736 0.2069513
[2,] 0.2069513 0.1573683

$vcm[[5]]
           [,1]       [,2]
[1,] 0.16187880 0.05227342
[2,] 0.05227342 0.03537038

$vcm[[6]]
           [,1]       [,2]
[1,] 0.20722789 0.08898875
[2,] 0.08898875 0.06038903

$vcm[[7]]
          [,1]      [,2]
[1,] 0.2740633 0.3502412
[2,] 0.3502412 0.5216565

$vcm[[8]]
          [,1]      [,2]
[1,] 0.1301925 0.1012023
[2,] 0.1012023 0.1091810

$vcm[[9]]
          [,1]       [,2]
[1,] 0.1925424 0.10061139
[2,] 0.1006114 0.07828174

7.5 項目特性曲線 (ICC)

能力θに対する項目の正答確率の変化を可視化　

ex3 <- est(resp = df_nine[, 1:9], # テストデータを指定する引数
  model = "2PL",        # 2PLMを仮定
  engine = "ltm")　# ltmパッケージを利用して項目パラメーターを推定すると指定

ex3

$est
        [,1]       [,2] [,3]
Q15 1.116352 -2.4774085    0
Q20 1.138009 -0.5955641    0
Q27 1.831933 -2.1289528    0
Q29 1.856599 -2.1700494    0
Q31 1.598930 -0.9810200    0
Q33 1.709964 -1.4333128    0
Q34 1.402326 -2.8348121    0
Q35 1.331616 -1.6712500    0
Q36 1.638347 -1.6180664    0

$se
           [,1]      [,2] [,3]
 [1,] 0.3674157 0.6197462    0
 [2,] 0.2972607 0.1851252    0
 [3,] 0.5913987 0.3908171    0
 [4,] 0.5998113 0.3966968    0
 [5,] 0.4023416 0.1880701    0
 [6,] 0.4552229 0.2457418    0
 [7,] 0.5235105 0.7222579    0
 [8,] 0.3608220 0.3304256    0
 [9,] 0.4387965 0.2797887    0

$vcm
$vcm[[1]]
          [,1]      [,2]
[1,] 0.1349943 0.2106149
[2,] 0.2106149 0.3840853

$vcm[[2]]
           [,1]       [,2]
[1,] 0.08836392 0.03069508
[2,] 0.03069508 0.03427134

$vcm[[3]]
          [,1]      [,2]
[1,] 0.3497524 0.2012777
[2,] 0.2012777 0.1527380

$vcm[[4]]
          [,1]      [,2]
[1,] 0.3597736 0.2069513
[2,] 0.2069513 0.1573683

$vcm[[5]]
           [,1]       [,2]
[1,] 0.16187880 0.05227342
[2,] 0.05227342 0.03537038

$vcm[[6]]
           [,1]       [,2]
[1,] 0.20722789 0.08898875
[2,] 0.08898875 0.06038903

$vcm[[7]]
          [,1]      [,2]
[1,] 0.2740633 0.3502412
[2,] 0.3502412 0.5216565

$vcm[[8]]
          [,1]      [,2]
[1,] 0.1301925 0.1012023
[2,] 0.1012023 0.1091810

$vcm[[9]]
          [,1]       [,2]
[1,] 0.1925424 0.10061139
[2,] 0.1006114 0.07828174

• 9 項目の特性曲線を描いてみる

P3 <- irf(ip = ex3$est) # irf()関数を使って正答確率を計算  

plot(x = P3,     # xは引数、irf関数で推定した結果を指定する
  co = NA,      # ICCの色を指定/項目毎に異なる色でICCを描く
  label = TRUE) # 各ICCに項目の番号がつく

abline(v = 0, lty    = 2) # x = 0 の縦点線を引く

横軸・・・潜在特性値 \(θ\) (Ability)
縦軸・・・正答確率 (Probability of a correct response)
- 曲線上の数字は「項目番号」ではなく「列番号」

識別力 a の解釈：

・項目が能力をどれだけ区別できるかを示すパラメータ

適切な識別力 a の大きさ

0.3 〜 2.0

（出典：芝祐順編『項目反応理論』p.34）

• 識別力を可視化

# 必要なパッケージ
library(irt)
library(ggplot2)
library(dplyr)

# モデル推定（再掲）
ex3 <- est(resp = df_nine[, 1:9], 
           model = "2PL",        
           engine = "ltm")

# 識別力（1列目）を取り出し
disc <- ex3$est[, 1]

# データフレームに変換して並び替え（小さい順に変更！）
disc_df <- data.frame(
  Item = names(disc),
  Discrimination = disc
) %>%
  arrange(Discrimination) %>%  # ★ここを修正
  mutate(Item = factor(Item, levels = Item))  # 並び替えを反映

# グラフ描画
ggplot(disc_df, aes(x = Item, y = Discrimination)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(aes(label = round(Discrimination, 2)), vjust = -0.5, size = 4) +
  labs(title = "識別力が小さい順に並べた項目", x = "項目", y = "識別力") +
  theme_minimal() +
  theme_bw(base_family = "HiraKakuProN-W3") # 文字化け対策 

結論
・9つの項目の識別力は 0.3〜2.0の間に収まっている

困難度 b の解釈：

適切な困難力 b の大きさ

−3 〜 3

（出典：芝祐順編『項目反応理論』p.34）

• 困難度の値順に可視化する

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定
ex3 <- est(resp = df_nine[, 1:9], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出
difficulty_df <- data.frame(
  Item = rownames(ex3$est),
  Difficulty = ex3$est[, 2]
)

# 3. 困難度で並べ替え（昇順）
difficulty_df <- difficulty_df %>%
  arrange(Difficulty) %>%
  mutate(Item = factor(Item, levels = Item))  # 項目順を固定する

# 4. ラベル位置を調整（棒の先に表示するため）
difficulty_df <- difficulty_df %>%
  mutate(label_position = ifelse(Difficulty >= 0, Difficulty + 0.2, Difficulty - 0.2))

# 5. プロット
ggplot(difficulty_df, aes(x = Difficulty, y = Item)) +
  geom_bar(stat = "identity", aes(fill = Difficulty > 0), width = 0.6) +
  geom_text(aes(label = round(Difficulty, 2), x = label_position), 
            size = 3, hjust = ifelse(difficulty_df$Difficulty >= 0, 0, 1)) +
  # ★ ここで赤い点線を引く ★
  geom_vline(xintercept = c(-3, 3), color = "red", linetype = "dashed", size = 0.2) +
  scale_fill_manual(values = c("skyblue", "salmon"), guide = "none") +
  labs(
    title = "困難度（Difficulty）の分布（-3と3に赤い点線）",
    x = "困難度",
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans")

標準誤差（識別力）の解釈：

標準誤差（識別力）の範囲

0.1 〜 0.4（識別力 a の標準誤差）

（出典：芝祐順編『項目反応理論』p.34）

library(irt)
library(ggplot2)
library(dplyr)

# モデル推定
ex3 <- est(resp = df_nine[, 1:9], 
           model = "2PL",        
           engine = "ltm")

# 識別力の標準誤差（標準偏差）を取り出す
disc_se <- ex3$se[, 1]  # 2列目が識別力SE


# 項目名を取得
item_names <- rownames(ex3$est)

# データフレームに変換して、小さい順に並べ替え
disc_se_df <- data.frame(
  Item = item_names,
  Difficulty_SE = disc_se
) %>%
  arrange(Difficulty_SE) %>%
  mutate(Item = factor(Item, levels = Item))  # 並び順を反映

# ラベルを棒グラフより少し右にずらすためオフセット作成
disc_se_df <- disc_se_df %>%
  mutate(label_position = Difficulty_SE + 0.02)

# グラフ描画
ggplot(disc_se_df, aes(x = Difficulty_SE, y = Item)) +
  geom_bar(stat = "identity", fill = "darkgreen") +
  geom_text(aes(x = label_position, label = round(Difficulty_SE, 3)), 
            hjust = 0, size = 4) +  # hjust=0で左寄せ
  geom_vline(xintercept = c(0.2, 0.5), color = "red", linetype = "dashed", size = 0.4) +
  labs(
    title = "識別力の標準誤差（SE）の値順に並べた9つの項目", 
    x = "識別力の標準誤差", 
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme_bw(base_family = "HiraKakuProN-W3")

標準誤差（困難度）の結果

標準誤差の規模	項目例	評価
小さい（安定）	Q20〜Q33	安心して使える
大きい（不安定）	Q29, Q27, Q34	注意。改善または除外を検討

標準誤差（困難度）の解釈：

標準誤差（困難度）の範囲

0.2 〜 0.5（困難度 b の標準誤差）

（出典：芝祐順編『項目反応理論』p.34）

library(irt)
library(ggplot2)
library(dplyr)

# モデル推定
ex3 <- est(resp = df_nine[, 1:9], 
           model = "2PL",        
           engine = "ltm")

# 困難度の標準誤差（標準偏差）を取り出す
disc_se <- ex3$se[, 2]  # 1列目が困難度の標準誤差

# 項目名を取得
item_names <- rownames(ex3$est)

# データフレームに変換して、小さい順に並べ替え
disc_se_df <- data.frame(
  Item = item_names,
  Discrimination_SE = disc_se
) %>%
  arrange(Discrimination_SE) %>%
  mutate(Item = factor(Item, levels = Item))  # 並び順をグラフに反映

# ★ ラベルを棒グラフより少し右にずらすためにオフセットを作成
disc_se_df <- disc_se_df %>%
  mutate(label_position = Discrimination_SE + 0.02)  # 0.02だけ右にずらす

# グラフ描画（横向き棒グラフ：左小 → 右大）
ggplot(disc_se_df, aes(x = Discrimination_SE, y = Item)) +
  geom_bar(stat = "identity", fill = "purple") +
  geom_text(aes(x = label_position, label = round(Discrimination_SE, 3)), 
            hjust = 0, size = 4) +  # ★hjust=0で左寄せ
  # ★ ここで赤い点線を引く ★
  geom_vline(xintercept = c(0.2, 0.5), color = "red", linetype = "dashed", size = 0.4) +
  labs(
    title = "困難度の標準誤差（SE）の値順に並べた9項目", 
    x = "困難度の標準誤差", 
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme_bw(base_family = "HiraKakuProN-W3")  # 文字化け対策

標準誤差（困難度）の結果

標準誤差の規模	項目例	評価
小さい（安定）	Q20〜Q29など	安心して使える
中程度（不安定）	Q15, Q34	除外を検討
大きい（推定不可）		除外

• 困難度（難易度）の標準誤差（SE）が大きいということは・・・
  ▶ その項目の困難度の推定が不安定である

標準誤差が大きいと次のような問題が起こる

1. 困難度の推定値が信頼できない

• 例えば、困難度 = -1.2 と推定されても
  
• 標準誤差が大きいと「本当は-1.8かもしれないし、-0.6かもしれない」とブレる
  
• 困難度がはっきりしない項目は、受験者の能力を正確に測る力が弱い

2. テスト全体の信頼性が下がる

• ブレブレの項目がテストに多いと
  →　テスト全体で測っている学力θの精度も落ちる

3. IRTモデルのフィット（適合度）が悪くなる

• 困難度が不安定な項目が多いと、モデル全体にとっても悪影響
  
• 項目反応曲線（ICC）がきれいな S 字ではなく変な形になる

分析結果（理想的な基準）

• 項目パラメタにそれぞれ適切な値を適応して項目の採用・不採用を検討する

適切な識別力、困難度、標準誤差
0.3 〜 2.0（識別力）
−3 〜 3（困難度）
0.1 〜 0.4（識別力 a の標準誤差）
0.2 〜 0.5（困難度 b の標準誤差）
（出典：芝祐順編『項目反応理論』p.34）

library(dplyr)

# --- estとseを読み込み（推定済みex3を使用）---
est <- ex3$est
se <- ex3$se

# --- データフレーム作成（列番号付き）---
irt_result <- data.frame(
  Col_No = 1:length(item_names),  # 列番号
  Item = item_names,              # 項目番号
  a = est[, 1],                # 識別力
  b = est[, 2],                # 困難度
  SE_a = se[, 1],              # 識別力の標準誤差
  SE_b = se[, 2]               # 困難度の標準誤差
)

# --- 判断と理由を付与 ---
irt_result <- irt_result %>%
  mutate(
    # 判断
    Judgment = case_when(
      a < 0.3 | a > 2 |
      b < -3 | b > 3 |
      SE_a < 0.1 | SE_a > 0.4 |
      SE_b < 0.2 | SE_b > 0.5 ~ "検討or削除",
      TRUE ~ "問題ない"
    ),
    
    # 理由
    Reason = case_when(
      a < 0.3 ~ "aが低",
      a > 2 ~ "aが高",
      b < -3 ~ "bが低",
      b > 3 ~ "bが高",
      SE_a < 0.1 ~ "SE_aが小",
      SE_a > 0.4 ~ "SE_aが大",
      SE_b < 0.2 ~ "SE_bが小",
      SE_b > 0.5 ~ "SE_bが大",
      TRUE ~ ""
    )
  ) %>%
  mutate(
    a = round(a, 3),
    b = round(b, 3),
    SE_a = round(SE_a, 3),
    SE_b = round(SE_b, 3)
  ) %>%
  arrange(factor(Judgment, levels = c("検討or削除", "問題ない")))

DT::datatable(irt_result)

• 「問題ない」・・・1項目
  
• 「検討or削除」・・・8項目

7.6 潜在特性値の推定（項目側）: ltm::ltm()

「項目側」のθ推定（モデルのフィット、aとbの推定）

• IRTモデルでは、観測された回答データ（正答/誤答）を使って
  →　受験者の学力θ（潜在特性値）を推定する
  
• その代表的な方法の一つが 最尤推定法（MLE: Maximum Likelihood Estimation）

最尤推定法の基本的な考え方 ・学力 θ の候補値を −3 から 3 まで 0.1 ごとに調べる
・どの θ のときに「そのパターンの回答が起こる可能性（尤度）」が最も高くなるかを探す
・その「最も可能性が高いθ」が、学力（潜在特性）の推定値

学力θを推定する確率計算

• 学力θを推定する確率計算

• ltmパッケージの中に入っているデータを使う
  • The Law School Admission Testへの解答結果を採点

• 受験者数: 220人
  

• 項目数：36項目
  

• 正答なら1、誤答なら０

観測された回答データ（正答/誤答）

DT::datatable(df2)

• theta（学力θのこと）の値が −3 から 3 までの範囲で、0.1 ごとの項目ごとの正答率を計算してみる

library(ltm)

# 2PLモデルによるIRT分析（5項目に限定）
mod <- ltm(df2[, 2:37] ~ z1, IRT.param = TRUE)

# 能力値 θ の範囲指定
theta_vals <- seq(-3, 3, by = 0.1)

# モデルから項目パラメーター（識別力 a、困難度 b）を抽出
coefs <- coef(mod)  # 戻り値は列1 = Dffclt (b), 列2 = Discr (a)

# 正答確率を格納するデータフレームの作成
icc_df <- data.frame(theta = theta_vals)

# 各項目に対して正答確率を計算（2PLモデルの公式）
for (i in 1:nrow(coefs)) {
  b <- coefs[i, 1]  # 困難度 b
  a <- coefs[i, 2]  # 識別力 a
  P_theta <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC式
  icc_df[[paste0("Item", i)]] <- round(P_theta, 4)
}

• theta（学力θのこと）の値が −3 から 3 までの範囲で、0.1 ごとの項目ごとの正答率が計算できた

item1 から item5 までの正答確率

DT::datatable(icc_df)

• theta（学力θのこと）の値が 0.1 ごとの項目ごとの正答率をすべて確認できる

• この結果を項目特性曲線として可視化してみる

plot(mod, type = "ICC", items = 2:36)

# 縦の点線を追加（θ = -3）
abline(v = -3, col = "red", lty = 2, lwd = 1)

• theta（学力θ） = -3 に赤い点線を引いた
  
• theta = -3 の場合、item1 の数値 0.5737
  →　黒線 (item1) は Probability = 0.57 を示している
  
• theta = -3 の場合、item3 の数値 0.0815
  →　緑線 (item3) は Probability = 0.0815 を示している

7.7. テスト情報曲線 (TIC)

テスト情報曲線でわかること

このテストはどの能力レベルを正確に測れているか？

• IRTではテストの推定精度を計算できる
  = 指定した学力θがどの程度の「誤差」を含むかがわかる
  
• テストの測定精度のことを「テスト情報量」とも呼ぶ
  
• テスト情報量は学力θごとに計算できる

テスト情報関数 ・2パラメタ・ロジスティック・モデルの場合のテスト情報量は次の式で表せる

\[I(\theta) = 1.7^2\sum_{j=1}^na_j^2P_j(\theta)Q_j(\theta)\]

変数	内容
\(I(θ)\)	テスト情報量
1.7	定数
\(a_j\)	項目\(j\) の識別力
\(P_j(\theta)\)	学力θに応じた項目 \(j\) に正答する確率
\(Q_j(\theta)\)	学力θに応じた項目 \(j\) に誤答する確率

• IRTにおける措定精度は、個々の学力θに応じて計算できる

• 学力θごとに計算したテスト情報量
  →　テスト情報曲線 (Test Information Curve，TIC) で可視化する

• 「テストがどの能力レベル（θ）でどれくらい正確に測れているか」を可視化できる
  • 潜在特性値の値毎にテスト情報量を計算しプロットする
  

• tif関数を使ってテスト情報量を計算
  →　plot関数で TIC を作成

• 様々な潜在特性値\(θ\)（能力）におけるテスト情報量（Information）が I

• irtoysパッケージの tif()関数を使って計算する

II <- irtoys::tif(ip = ex3$est) # データに対し2PLMを仮定
                           # x: tif関数で推定した結果を指定する引数
plot(x = II)                # ip: テストに含まれる各項目の項目パラメーターを指定する引数 

・横軸・・・潜在特性値\(θ\)（学力）
・縦軸・・・テスト情報量（測定の精度＝標準誤差の逆数）
・実線・・・テスト情報曲線
→　当該潜在特性値におけるテスト情報量をつなぎ合わせたもの

テスト情報曲線 (TIC) でわかること 1. どの能力レベルを正確に測れているか？
・情報量が高いところ → その学力θレベルでテストが精度が高い
→ \(\theta = −2\)付近の情報量が最も高い
→ \(\theta = −2\)付近のレベルでテストが高精度

・情報量が低いところ
→ その学力θレベルではテストの精度が低い
→ \(\theta = 4\)の情報量が最も低い
→ \(\theta = 4\)のレベルでテストが低精度

例：TICがθ = 0の周辺で高ければ
→「平均的な人を測るのに最適なテスト」だといえる

2. テストの設計意図が見える
・TICがどこで山になるかを見ることで、そのテストがどんな対象向けかが分かる
・TICが\(\theta = −2\)付近で山になっている
→　このテストは比較的能力の低い人向け

・TICがどこで山になるかを見ることで
→ そのテストがどんな対象向けかが分かる

TICの山の位置	意味	対象
θ = 0 周辺	平均均的な受験者向け	一般的な学力テスト
θ > 0（右寄り）	高能力者向け	難関資格・上級試験
θ < 0（左寄り）	初級者・低能力者向け	基礎力診断など

・ここではTICの山の位置が左より
→　平均均的な受験者向けだとわかる

3. 信頼性の高さ（精度）もわかる
・情報量が高い＝その範囲の 標準誤差（SE）が小さい
・標準誤差との関係：

\[SE(\theta) = \frac{1}{\sqrt{{I(\theta)}}}\]

・つまり、情報量が大きいと \(\theta\) の推定がブレにくい（= 信頼できる）

4. テスト情報量の基準を設定する：4以上

・例えば、テスト情報量の基準を「4以上」と設定する
→　情報精度の基準を満たしている学力θの範囲がわかる

項目適合度の結果のポイント どの潜在特性値 \(\theta\) の値で、情報量が最大になるか
・潜在特性値が −2 の辺りでテスト情報量が最大
⇒ 潜在特性値\(θ\)（能力）が低い（ −2 の辺り ）受験者の推定精度が最も高い

7.8 局所独立性の検討: irf() & cor()

7.8.1 局所独立性とは

• IRTに基づきテストデータや質問紙への回答を分析する
  ⇒ 局所独立性が仮定されている
  • IRTにおける局所独立性の検討
  
• 潜在特性値 \(\theta\) の値を固定したとき、項目相互間で解答が独立に生じている仮定を満たすこと
  
• 「局所独立の仮定」とは「一次元性の仮定」のこと

一次元性の仮定 各項目の正誤が、潜在特性値 \(\theta\) の値の大小によってのみばらつく

• つまり「受験者が正解するかどうかは、受験者の能力だけで決まる」という仮定

• 局所独立性の検討は \(Q_3\)統計量に基づいて行われることが多い
• \(Q_3\)統計量とは、各項目への回答（観測値）からその期待値を引き
→　得られた残差得点間の相関を求めることで得られる
- \(Q_3\)統計量は、各項目への反応（= 観測値） からその期待値（= 項目反応モデルにより計算される正答確率）を引き
→　得られた残差得点間の相関を求めることで得られる統計量
- その絶対値が 0 に近いほど、項目反応間に局所独立性を仮定できる

• たとえば今の場合、item1 の残差得点\(d_1\)は次の式で表せる

\[d_1 = u_1 - \hat{P_1(\theta)}\]

• \(u_1\): Q1 への解答結果（正答なら１、誤答なら０）
  
• \(\hat{P_1(\theta)}\): 項目パラメーターと潜在特性の推定値から計算される正答確率

7.8.2 Rで \(Q_3\) 統計量を計算する方法

step 1: mlebme関数を利用して潜在特性値を推定

• mlebme関数を利用して潜在特性値を推定
⇒ mlebme関数はirtoysパッケージの中に入っている関数

head(mlebme(resp = df_nine[, 1:9],　# テストデータを指定
  ip = ex3$est,       # データに対し2PLMを仮定
  method = "BM"))           # 最尤推定法 (ML) による潜在特性値の推定を指定

             est       sem n
[1,] -1.62540218 0.4295256 9
[2,]  0.03210580 0.6495970 9
[3,] -1.26007561 0.4500825 9
[4,] -1.22422032 0.4530553 9
[5,]  0.02299343 0.6477295 9
[6,]  0.58531124 0.7615398 9

• resp: テストデータを指定する引数
  
• ip: テストに含まれる各項目の項目パラメーターを指定する引数
  
• ex3$est と指定
  →　データに対し2PLMを仮定したときの項目パラメーターの推定値が各項目の項目パラメーターとして指定
  
• 項目パラメーターを推定した後に、潜在特性値の推定を行う
  
• method: どの推定法を用いて潜在特性値を推定するか
  
• ML と指定
  →　最尤推定法による潜在特性値の推定を指定
  ⇒ 全問正答/誤答の受験者がいる場合、コれらの受験者の推定値が求まらない
  
• BM と指定
  →　潜在特性値が標準正規分布に従っていることを仮定しているという情報を加味して推定が行えるようになり（事前情報）
  →　全問正答/誤答の受験者に対しても推定値を得ることができる（ベイズ推定）

theta.est <- mlebme(resp = df_nine[,1:9],
  ip = ex3$est,
  method="BM")

DT::datatable(theta.est)

• １列目：推定値 (est)
  
• ２列目：標準誤差 (sem)
  
• ３列目：解答した項目 (n)

Step 2. irf関数を利用して正答確率 ($f) を推定

• irf関数では2PLMを仮定
• ex3$est と指定
→　データに対し2PLMを仮定したときの項目パラメーターの推定値を各項目の項目パラメーターとして指定
• theta.est[, 1] と指定
→　データに対し2PLMを仮定したときの潜在特性の推定値を指定

⇒ 結果はPとして保存

P <- irf(ip = ex3$est, # 項目パラメーターを指定
  x = theta.est[, 1])  # 各受験者の潜在特性値を指定

変数	内容
$x	各受験者の潜在特性値\(\theta\)（能力）
$f	正答確率の推定値
行	受験者 (220名)
列	項目 (Q15, Q20, Q27, Q29, Q31, Q33, Q34, Q35, Q36)

Replace this figure!

Step 3. 推定された正答確率とテストデータより、残差得点 \(d\) を算出

• \(d_j = u_j - \hat{P_j(\theta)}\) の値を計算して \(d\) として保存 (1 ≦ j ≦9 )
• P$fと指定することで，正答確率の推定値が抽出される
  ⇒ テストデータ (df_nine,1:9]) から正答確率を引いた残差得点を \(d\) に保存

d <- df_nine[, 1:9] - P$f # P$f と指定して正答確率の推定値を得る  

• 220人の受験生の最初の6人分の計算結果を表示させてみる

head(d)

          Q15        Q20         Q27         Q29         Q31         Q33
1  0.27865585  0.7635000 0.284454670 0.266746778 -0.26302167 -0.41861396
2 -0.94275687  0.3286508 0.018726442 0.016487654  0.16521628  0.07545168
3  0.20440746 -0.3194675 0.169141358 0.155847381 -0.39026691 -0.57352051
4 -0.80202499 -0.3284032 0.160110065 0.147288882  0.59600831 -0.58844497
5  0.05779458 -0.6690572 0.019035671 0.016764250  0.16723559  0.07654586
6  0.03170468  0.2068782 0.006879199 0.005966569  0.07554652  0.03071599
           Q34         Q35         Q36
1 -0.845009896 -0.51525818 -0.49699539
2  0.017630213  0.09378980  0.06276394
3  0.099008149 -0.63356147  0.35743531
4  0.094612377  0.35542619 -0.65594073
5  0.017852899  0.09482622  0.06364790
6  0.008193977  0.04720625  0.02634196

受験者1のQ15 に関する残差得点 \(d_{11}\)をチェック ・例えば、受験者1のQ15 に関する残差得点 \(d_1\) は 0.27865585
・受験者1の Q15 への回答 \(u_{ij} = u_{11}\) を確認してみる

df_nine[1,1]

[1] 1

• 受験者1の item1 への回答は 1（正答）
  
• 受験者1の Q15 に関する正答確率の推定値 \(\hat{P_{11}(\theta)}\) は0.27865585（上の図を参照）
  →　受験者1の残差得点 \(d_{11}\) は

df_nine[1,1] - 0.27865585

[1] 0.7213441

Step 4. cor関数を利用して \(Q_3\) 統計量の値を計算

Q3 <- cor(x = d, 
  y = d, 
  use = "pairwise.complete.obs")

Q3

            Q15         Q20         Q27         Q29         Q31         Q33
Q15  1.00000000  0.09118227  0.14351291 -0.18093023 -0.16191130 -0.08141794
Q20  0.09118227  1.00000000 -0.02772476 -0.04527975 -0.07544948 -0.11149725
Q27  0.14351291 -0.02772476  1.00000000  0.16986533 -0.18460036 -0.22923594
Q29 -0.18093023 -0.04527975  0.16986533  1.00000000 -0.21958995 -0.08977776
Q31 -0.16191130 -0.07544948 -0.18460036 -0.21958995  1.00000000 -0.10826242
Q33 -0.08141794 -0.11149725 -0.22923594 -0.08977776 -0.10826242  1.00000000
Q34 -0.06464632 -0.21220045 -0.10987949  0.01530725 -0.06380572  0.07519163
Q35 -0.13450831 -0.17414860 -0.05136856 -0.04889832 -0.12272164 -0.06661985
Q36 -0.07171384 -0.13488838 -0.20999653 -0.12779358 -0.01523171 -0.17267718
            Q34         Q35         Q36
Q15 -0.06464632 -0.13450831 -0.07171384
Q20 -0.21220045 -0.17414860 -0.13488838
Q27 -0.10987949 -0.05136856 -0.20999653
Q29  0.01530725 -0.04889832 -0.12779358
Q31 -0.06380572 -0.12272164 -0.01523171
Q33  0.07519163 -0.06661985 -0.17267718
Q34  1.00000000 -0.01832574 -0.01656062
Q35 -0.01832574  1.00000000 -0.04327039
Q36 -0.01656062 -0.04327039  1.00000000

局所独立性の検討のポイント

各項目の正誤が、潜在特性値 \(\theta\) の値の大小によってのみばらつくかどうか
。\(Q_3\) の絶対値が 0 に近いほど、項目反応間に局所独立性を仮定できる
・\(Q_3\) の値の絶対値が 0.2以上の場合　→　問題あり: 局所依存性の疑い（Chen & Thissen, 1997）
・\(Q_3\) の値の絶対値が 0.2以下の場合　→　問題なし

Step 5. \(Q_3\) 統計量の値に関するヒートマップを作成する

• Q3 の絶対値をとり
  
• 対角成分（自己相関1.0）をNAにして除外し
  
• 絶対値が0.2以上の個数を数え、全体の割合を計算してみる

# Q3の絶対値をとる
Q3_abs <- abs(Q3)

# 対角成分（自己相関1.0）をNAにして除外する
diag(Q3_abs) <- NA

# 絶対値が0.2以上の個数を数える
count_0.2_or_more <- sum(Q3_abs >= 0.2, na.rm = TRUE)

# 全体の要素数（対角除外後の数）
total_elements <- sum(!is.na(Q3_abs))

# 割合を計算
proportion_0.2_or_more <- count_0.2_or_more / total_elements

# 結果を表示
cat("絶対値が0.2以上の個数:", count_0.2_or_more, "\n")

絶対値が0.2以上の個数: 8 

cat("割合:", proportion_0.2_or_more, "\n")

割合: 0.1111111 

• 𝑄3の値の絶対値が 0.2以上のケースは66（全体の5%程度）
  
• 得られた結果をヒートマップで可視化してみる

library(ggplot2)
library(reshape2)

# Q3行列を絶対値に
Q3_abs <- abs(Q3)
diag(Q3_abs) <- NA  # 対角成分除外

# meltしてロング形式に
Q3_long <- melt(Q3_abs, varnames = c("Item1", "Item2"), value.name = "Q3_value")

# 0.2超えのフラグをつける
Q3_long$Violation <- Q3_long$Q3_value > 0.2

# ヒートマップ作成
ggplot(Q3_long, aes(x = Item1, y = Item2, fill = Violation)) +
  geom_tile(color = "white") +
  scale_fill_manual(values = c("white", "deeppink")) +
  labs(
    title = "局所独立性違反ヒートマップ（ピンク＝違反）",
    x = "項目",
    y = "項目"
  ) +
  theme_bw(base_family = "HiraKakuProN-W3") +   # 先に白背景＋日本語フォント
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1),  # ←ここで傾ける
    axis.text.y = element_text(hjust = 1)
  )

• どの項目同士のペアのQ3値が高いのかを一覧にしてみる

library(dplyr)
library(reshape2)

# --- Q3絶対値＋小数点第3位四捨五入
Q3_abs <- abs(Q3)
Q3_abs <- round(Q3_abs, 3)
diag(Q3_abs) <- NA  # 対角成分はNA

# --- ロング形式に変換
Q3_long <- melt(Q3_abs, varnames = c("Item1", "Item2"), value.name = "Q3_abs")

# --- ダブりを除去するため、Item1とItem2を並び替え
Q3_long <- Q3_long %>%
  filter(!is.na(Q3_abs)) %>%
  mutate(Item_min = pmin(as.character(Item1), as.character(Item2)),
         Item_max = pmax(as.character(Item1), as.character(Item2))) %>%
  select(Item1 = Item_min, Item2 = Item_max, Q3_abs)

# --- さらに重複を除去
Q3_long <- distinct(Q3_long)

# --- 0.2以上違反だけ取り出して、大きい順に並べる
violations <- Q3_long %>%
  filter(Q3_abs >= 0.2) %>%
  arrange(desc(Q3_abs))

DT::datatable(violations)

局所独立性の検討のポイント

各項目の正誤が、潜在特性値 \(\theta\) の値の大小によってのみばらつくかどうか
・Q_3 の絶対値が 0 に近いほど、項目反応間に局所独立性を仮定できる

・\(Q_3\) の値の絶対値が 0.2以上の場合　→　問題あり
→　局所依存性の疑い（Chen & Thissen, 1997）

・\(Q_3\) の値の絶対値が 0.2以下の場合　→　問題なし

・\(Q_3\) の値の絶対値の値を計算してみた
→　\(Q_3\) の値の絶対値が0.2以上の個数が 66 　→　問題あり
→　局所独立性が成立していないことを示唆
→　対策が必要

• Q3行列は「IRTモデルで推定された正答確率」と「実際の回答データ」 のズレ（= 残差）から作った残差相関行列
  
• 従って、Q3が0.2以上というのは「局所独立性の仮定が破れている」ということ →　次のいずれかの可能性が高いことを示唆している

🔵 学力θで説明しきれない強い依存関係が残っている
🔵 何か別の理由で項目同士が結びついている

→　ここでは\(Q_3\) の値の絶対値が全て 0.2以下　→　問題なし
→　局所独立性が成立していることを示唆
→　次の分析に移る

・局所依存性：局所独立性が満たされていない状態のこと
⇒ \(Q_3\)はあくまで一つの基準に過ぎないので注意が必要

7.9 項目適合度の検討: itf()

各項目がIRTモデルにどれだけ適合しているかを検証

・「各項目が理論モデル（例：2PLモデル）にちゃんと従っているかどうか」を評価するものが項目適合度 (Item Fit)
• IRTにおいては項目反応モデルへの適合度を検討することも重要

• ここでは、itf関数を利用して Q15 の項目適合度を検討してみる

DT::datatable(df_nine)

• resp はテストデータを指定する引数

irtoys::itf(resp = df_nine[, 1:9], # 応答データ [受験者, 項目]  
  item = 1,                     # 1番目の項目適合度の検討を指定する
  ip = ex2$est, 　　　　   # 2PLモデルによる推定された項目パラメータ
  theta = theta.est[, 1])       # 受験者ごとの能力推定値（θ）

Statistic        DF   P-value 
 6.007858  5.000000  0.305455 

• 図は、当てはめた項目反応モデルとデータとの乖離度を示している
  
• 横軸・・・潜在特性値 (Ability) 　　
• 縦軸・・・正答確率 (Proportion right)
  
• 実線・・・当てはめた項目反応モデルに基づく正答確率の予測値
  
• 円・・・潜在特性の推定値に基づき受験者を群分けした際の群ごとの実際の正答率
  • 潜在特性値\(θ\)（能力）は平均０、分散１

項目適合度でわかること
その項目がIRTモデルに「合っているか」どうか
・各項目の「実際のデータによる反応パターン」と、IRTモデルが「理論的に予測する反応パターン」を比較
→　もしズレていたら、モデルの仮定がその項目には合っていないということ

なぜ重要？

・モデルに合っていない項目を使うと、潜在特性値\(θ\)（能力）の推定が不正確になる可能性あり
・項目の品質をチェックし、不適切な項目を修正・削除する判断材料になる
・バイアスの検出（DIF：Differential Item Functioning）の手がかりにもなる

適合度が悪いときに考えられること

現象	可能性
実際の正答率がモデルより低い	問題文がわかりにくい／迷いやすい選択肢
特定の能力層だけ挙動がおかしい	バイアスがある、ミスリードされやすい項目
正答率がランダムに近い	推測が強く影響（cパラメータが不十分）
認知的に複雑すぎる	単一の「能力θ」では説明できない

適合度を判断する指標： S-X²統計量（Orlando & Thissenの項目適合度指標）
・より精度の高い適合度検定（特に2PLや3PLに使う）
・能力をグループ（通常は10分位など）に分け
→　各グループでのモデルによる期待正答率と、実際の正答率の差を使う
→　カイ二乗型の統計量として適合度を評価
＝　これは、カイ二乗分布に従う統計量だが、S-X²特有の方法
→　通常のカイ二乗適合度検定とは区別される

結果の解釈

\[ 帰無仮説: 「当てはめたモデルがデータに適合している」\]

得られた結果:

• p値が有意水準 (0.05) よりも大きい: p-value = 0.305455
  →　帰無仮説は棄却できない
  → 「当てはめた項目反応モデルがデータに適合している」と判断される

• itf 関数を使用した際に出力される図において
  

• 実線と円の間の乖離が大きいほど、モデルがデータに当てはまっていないと判断

7.11 潜在特性値の推定（受験者側）

「受験者側」のθ推定（どのくらいできる人か）.

• Rで潜在特性値を推定する
  • mlebme関数を利用して潜在特性値を推定
  ⇒ mlebme関数はirtoysパッケージの中に入っている関数

head(mlebme(resp = df_nine[, 1:9],　# テストデータを指定
  ip = ex2$est,       # データに対し2PLMを仮定
  method = "BM"))           # 最尤推定法 (ML) による潜在特性値の推定を指定

             est       sem n
[1,] -1.62540218 0.4295256 9
[2,]  0.03210580 0.6495970 9
[3,] -1.26007561 0.4500825 9
[4,] -1.22422032 0.4530553 9
[5,]  0.02299343 0.6477295 9
[6,]  0.58531124 0.7615398 9

• resp: テストデータを指定する引数
  
• ip: テストに含まれる各項目の項目パラメーターを指定する引数
  
• ex1$est と指定
  →　データに対し2PLMを仮定したときの項目パラメーターの推定値が各項目の項目パラメーターとして指定
  
• 項目パラメーターを推定した後に、潜在特性値の推定を行う
  
• method: どの推定法を用いて潜在特性値を推定するか
  
• ML と指定
  →　最尤推定法による潜在特性値の推定を指定
  ⇒ 全問正答/誤答の受験者がいる場合、コれらの受験者の推定値が求まらない
  
• BM と指定
  →　潜在特性値が標準正規分布に従っていることを仮定しているという情報を加味して推定が行えるようになり（事前情報）
  →　全問正答/誤答の受験者に対しても推定値を得ることができる（ベイズ推定）

theta.est <- mlebme(resp = df_nine[,1:9],
  ip = ex2$est,
  method="BM")

DT::datatable(theta.est)

• １列目：推定値 (est)
  
• ２列目：標準誤差 (sem)
  
• ３列目：解答した項目 (n)

この分析結果からわかること

• この表には、各受験者ごとの以下の情報が示されている

項目	意味
est	能力値（θ）の推定値
sem	標準誤差（Standard Error of Measurement）
n	推定に用いた項目数（ここでは全員5項目）

主な傾向：

θが小さい（マイナス）人が多い

→　多くの受験者のθは-1.5〜0の範囲にある
→　平均的またはやや低めの学力層が多い

θの推定精度は一様に低い

• 標準誤差（sem）が約0.80前後でほぼ一定
  →　推定精度は高くはない
  
• これは項目数が少ない（n = 5）ため

受験者側からθを推定する意味

• 受験者の「潜在的な能力」を数値化する
  → 正答率ではなく「どの難しさの問題を解けたか」を加味して学力を測る

8. 一次元性を満たさない場合

• 一元性が厳密に満たされていない場合でも
  →　IRT分析では「実用的な洞察」や「問題傾向の把握」が可能

一次元性を満たさなくても、IRT分析が有益な理由

1. おおまかな弁別力・困難度の把握ができる

• 識別力や困難度などのIRTパラメータは、完全な一次元性がなくても指標として参考になる
  
• たとえば「この項目は難しい」「この項目は能力の差をよく見分ける」などの傾向は読み取れる

2. 問題作成や改訂のヒントになる

• 複数の下位因子があることがわかれば、それに応じてテスト構成を再検討したり、項目をグループに分けて分析する足がかりになる
  
• 特定の項目群が別因子に偏っている場合、「この項目は別の能力を測っているかも」という改善の材料になる

3. 学習者の回答パターンの診断

• 次元性がやや低くても、IRTモデルを使って学力推定（θ）や項目ごとの反応確率を出せば
  →　どのレベルの受験者がどこでつまずいているかが見えてくる

一次元性を満たすかどうかの確認 (df2)

DT::datatable(df2)

library(psych)

df2 <- read.csv("data/irt_was.csv")
item_df2 <- df2[, -1, -38] # 1列目と38列目は分析に含めない  

# 一次元性の検討（主成分＋平行分析）
fa.parallel(item_df2, 
  fa = "pc", 
  n.iter = 100) # 100回のシミュレーションを指定

Parallel analysis suggests that the number of factors =  NA  and the number of components =  7 

DT::datatable(df2)

項目パラメーターを推定 (ex2)

ex2 <- est(resp = df2[, 2:37], # テストデータを指定する引数
  model = "2PL",        # 2PLMを仮定
  engine = "ltm")　# ltmパッケージを利用して項目パラメーターを推定すると指定

ex2

$est
           [,1]        [,2] [,3]
Q1   0.83101563  -1.4837419    0
Q2   1.38529706  -1.8868385    0
Q3   0.96083273  -2.5781569    0
Q4   1.17398639   0.5348538    0
Q5   0.30308003   1.1720583    0
Q6   0.40541452  -0.8908177    0
Q7   0.84944311  -2.2112025    0
Q8   0.21668701  -1.1654871    0
Q9   1.02539301  -0.4665746    0
Q10  0.84664746   0.1589867    0
Q11  0.19059685   5.9706787    0
Q12  0.75312216  -0.4032208    0
Q13  0.08718103  -4.6166578    0
Q14  0.56079979   2.3682460    0
Q15  2.66234874  -1.5583871    0
Q16  0.99818390  -1.7964914    0
Q17  0.84008267   0.7663600    0
Q18  1.40381366  -1.1669284    0
Q19  0.48693593  -0.8133296    0
Q20  1.98871836  -0.4996180    0
Q21  1.35340407  -1.2622575    0
Q22  0.96305339  -1.1754989    0
Q23 -0.08172114 -19.8650474    0
Q24  0.03501333  37.1694908    0
Q25 -0.30777475  -2.6872822    0
Q26  0.02950535  23.9739215    0
Q27  1.55000393  -2.2765012    0
Q28  3.98038339  -1.9469339    0
Q29  0.76696232  -3.9245301    0
Q30  0.62852045   0.8365656    0
Q31  1.34152251  -1.1021836    0
Q32  0.15013542   3.8197754    0
Q33  1.08891672  -1.8915249    0
Q34  0.84665151  -4.0757930    0
Q35  0.79165359  -2.4259284    0
Q36  1.06302106  -2.1052554    0

$se
           [,1]        [,2] [,3]
 [1,] 0.2338576   0.3801713    0
 [2,] 0.3987228   0.3810303    0
 [3,] 0.3385227   0.7379581    0
 [4,] 0.2338250   0.1738089    0
 [5,] 0.1607692   0.7590312    0
 [6,] 0.1691979   0.4897666    0
 [7,] 0.2761040   0.6118341    0
 [8,] 0.1581844   1.0400471    0
 [9,] 0.2280855   0.1751453    0
[10,] 0.1980058   0.1896641    0
[11,] 0.1786044   5.5714769    0
[12,] 0.1942374   0.2162077    0
[13,] 0.1614898   8.6640289    0
[14,] 0.1947492   0.7913817    0
[15,] 0.7785051   0.2066898    0
[16,] 0.2852876   0.4183951    0
[17,] 0.2002724   0.2464511    0
[18,] 0.3247702   0.2081203    0
[19,] 0.1829278   0.3981055    0
[20,] 0.4084806   0.1163234    0
[21,] 0.3245696   0.2318685    0
[22,] 0.2397270   0.2718255    0
[23,] 0.2078260  50.4074626    0
[24,] 0.1857672 197.1842928    0
[25,] 0.1780322   1.5581067    0
[26,] 0.1656170 134.6250212    0
[27,] 0.5180435   0.5006750    0
[28,] 1.8721125   0.2602029    0
[29,] 0.4029644   1.7995226    0
[30,] 0.1808207   0.3242469    0
[31,] 0.3129724   0.2088459    0
[32,] 0.1631607   4.2265221    0
[33,] 0.3145162   0.4288678    0
[34,] 0.4816135   1.9908538    0
[35,] 0.2748202   0.7268334    0
[36,] 0.3270586   0.5094791    0

$vcm
$vcm[[1]]
           [,1]       [,2]
[1,] 0.05468939 0.07445014
[2,] 0.07445014 0.14453019

$vcm[[2]]
          [,1]      [,2]
[1,] 0.1589799 0.1342853
[2,] 0.1342853 0.1451841

$vcm[[3]]
          [,1]      [,2]
[1,] 0.1145976 0.2352867
[2,] 0.2352867 0.5445822

$vcm[[4]]
            [,1]        [,2]
[1,]  0.05467412 -0.01872764
[2,] -0.01872764  0.03020952

$vcm[[5]]
            [,1]        [,2]
[1,]  0.02584675 -0.09621237
[2,] -0.09621237  0.57612836

$vcm[[6]]
           [,1]       [,2]
[1,] 0.02862794 0.05677026
[2,] 0.05677026 0.23987132

$vcm[[7]]
          [,1]      [,2]
[1,] 0.0762334 0.1558379
[2,] 0.1558379 0.3743410

$vcm[[8]]
          [,1]      [,2]
[1,] 0.0250223 0.1302215
[2,] 0.1302215 1.0816979

$vcm[[9]]
           [,1]       [,2]
[1,] 0.05202298 0.01415844
[2,] 0.01415844 0.03067586

$vcm[[10]]
             [,1]         [,2]
[1,]  0.039206302 -0.008202823
[2,] -0.008202823  0.035972456

$vcm[[11]]
            [,1]       [,2]
[1,]  0.03189952 -0.9838605
[2,] -0.98386055 31.0413548

$vcm[[12]]
           [,1]       [,2]
[1,] 0.03772815 0.01413777
[2,] 0.01413777 0.04674576

$vcm[[13]]
           [,1]      [,2]
[1,] 0.02607895  1.375339
[2,] 1.37533850 75.065397

$vcm[[14]]
            [,1]       [,2]
[1,]  0.03792726 -0.1421844
[2,] -0.14218437  0.6262850

$vcm[[15]]
          [,1]       [,2]
[1,] 0.6060702 0.12211463
[2,] 0.1221146 0.04272067

$vcm[[16]]
           [,1]      [,2]
[1,] 0.08138904 0.1051536
[2,] 0.10515363 0.1750545

$vcm[[17]]
            [,1]        [,2]
[1,]  0.04010902 -0.03053288
[2,] -0.03053288  0.06073814

$vcm[[18]]
           [,1]       [,2]
[1,] 0.10547568 0.04825235
[2,] 0.04825235 0.04331407

$vcm[[19]]
           [,1]       [,2]
[1,] 0.03346259 0.04856354
[2,] 0.04856354 0.15848797

$vcm[[20]]
           [,1]       [,2]
[1,] 0.16685640 0.01270551
[2,] 0.01270551 0.01353113

$vcm[[21]]
           [,1]       [,2]
[1,] 0.10534545 0.05702124
[2,] 0.05702124 0.05376299

$vcm[[22]]
           [,1]       [,2]
[1,] 0.05746904 0.04876296
[2,] 0.04876296 0.07388908

$vcm[[23]]
             [,1]       [,2]
[1,]   0.04319165  -10.46569
[2,] -10.46568747 2540.91229

$vcm[[24]]
             [,1]        [,2]
[1,]   0.03450945   -36.61974
[2,] -36.61973784 38881.64534

$vcm[[25]]
            [,1]       [,2]
[1,]  0.03169548 -0.2635148
[2,] -0.26351485  2.4276966

$vcm[[26]]
             [,1]        [,2]
[1,]   0.02742898   -22.28111
[2,] -22.28110624 18123.89632

$vcm[[27]]
          [,1]      [,2]
[1,] 0.2683691 0.2366886
[2,] 0.2366886 0.2506755

$vcm[[28]]
          [,1]       [,2]
[1,] 3.5048053 0.37991626
[2,] 0.3799163 0.06770557

$vcm[[29]]
          [,1]     [,2]
[1,] 0.1623803 0.708381
[2,] 0.7083810 3.238282

$vcm[[30]]
            [,1]        [,2]
[1,]  0.03269611 -0.03853284
[2,] -0.03853284  0.10513604

$vcm[[31]]
           [,1]       [,2]
[1,] 0.09795175 0.04634588
[2,] 0.04634588 0.04361663

$vcm[[32]]
            [,1]       [,2]
[1,]  0.02662141 -0.6722509
[2,] -0.67225093 17.8634891

$vcm[[33]]
           [,1]      [,2]
[1,] 0.09892042 0.1201392
[2,] 0.12013920 0.1839276

$vcm[[34]]
          [,1]      [,2]
[1,] 0.2319515 0.9379179
[2,] 0.9379179 3.9634987

$vcm[[35]]
           [,1]      [,2]
[1,] 0.07552614 0.1869027
[2,] 0.18690270 0.5282868

$vcm[[36]]
          [,1]      [,2]
[1,] 0.1069673 0.1518345
[2,] 0.1518345 0.2595690

\[8.1〜8.4の分析は7.1〜7.4と同様なので省略\]

8.5 項目特性曲線 (ICC)

• ここでは 36 項目の「識別力a」と「困難度b」が計算できた
  
• 36 項目の特性曲線を描いてみる

P2 <- irf(ip = ex2$est) # irf()関数を使って正答確率を計算  

plot(x = P2,     # xは引数、irf関数で推定した結果を指定する
  co = NA,      # ICCの色を指定/項目毎に異なる色でICCを描く
  label = TRUE) # 各ICCに項目の番号がつく

abline(v = 0, lty    = 2) # x = 0 の縦点線を引く

横軸・・・潜在特性値 \(θ\) (Ability)
縦軸・・・正答確率 (Probability of a correct response)

• 例えば、Q25の曲線は左上から右下に下がっている → 能力の高まるほど正答率が下がる
  
• Q25の識別力は -0.30777475

識別力 a の解釈：

・項目が能力をどれだけ区別できるかを示すパラメータ

適切な識別力 a の大きさ

0.3 〜 2.0

（出典：芝祐順編『項目反応理論』p.34）

• 識別力を可視化
  
• 適切な識別力 a の上限と下限に赤の点線を引く

library(irt)
library(ggplot2)
library(dplyr)

# モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 識別力（aパラメータ）を取り出す
discrimination <- ex2$est[, 1]

# データフレーム化
disc_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = discrimination
) %>%
  arrange(Discrimination) %>%
  mutate(Item = factor(Item, levels = Item))  # 並び順固定

# ラベルを左右に分けるため、hjust列を作成
disc_df <- disc_df %>%
  mutate(
    hjust_pos = ifelse(Item %in% c("Q23", "Q25"), 1.1, -0.1)  # Q24とQ26だけ左側に、それ以外は右側に
  )

# グラフ描画
ggplot(disc_df, aes(x = Discrimination, y = Item)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  geom_text(
    aes(label = round(Discrimination, 2), hjust = hjust_pos),
    size = 3
  ) +
  geom_vline(xintercept = 0.3, color = "red", linetype = "dashed", size = 0.5) +
  geom_vline(xintercept = 2.0, color = "red", linetype = "dashed", size = 0.5) +
  labs(
    title = "識別力の大きさ順に並べた項目",
    x = "識別力",
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  coord_cartesian(xlim = c(min(disc_df$Discrimination) - 0.5, max(disc_df$Discrimination) + 0.5))  # 余白確保

識別力が 0.3 以下の項目の ICC

• Q8〜Q25だけを ICC曲線として描く

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. 識別力が0.3以下の項目だけ抽出
low_disc_items <- param_df %>%
  filter(Discrimination <= 0.3)

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%  # θ = 4 の位置
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +  # ★ ここでy軸を0〜1に固定 ★
  labs(
    title = "識別力が0.3以下の項目のICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")  # 凡例を消してスッキリ

識別力0.3以下の項目の評価
・識別力は、受験者のθ（能力）による正答確率の変化の鋭さを表す
・識別力が高いと　→　能力の違いに敏感に反応する
→　θが少し上がるだけで正答率が大きく上がる
・識別力が低い（0.3以下）と
→　θが変わっても正答確率があまり変わらない＝「区別がつかない」
・Q25の識別力は -0.30777475
→　識別力がマイナスの値　→　能力が高いほど正答率が低いということ
・識別力が小さくても、正答率が高ければ「基礎知識の確認問題」と位置づけて残すことも可能

結論：

・Q23とQ25・・・識別力がマイナス（異常値）なので除外
・Q32, Q11, Q26, Q24・・・識別力が小さので除外
・Q8とQ13・・・識別力が小さいが、正答率が比較的高いので「基礎知識の確認問題」として残すことも可能だが、基本的に除外

識別力が高すぎる項目の ICC

• 識別力が3.98 (Q28) と 2.66 (Q15) の ICC を確認する

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. 識別力が 0.2 以上の項目（＝Q30 & Q15）だけ抽出
low_disc_items <- param_df %>%
  filter(Discrimination > 2)  

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%  # θ = 4 の位置
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +  # ★ ここでy軸を0〜1に固定 ★
  labs(
    title = "Q28 (識別力 = 3.98) とQ15（識別力 = 2.66）のICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")  # 凡例を消してスッキリ

高すぎる識別力の問題点：

1. 現実を正しく反映していない可能性がある

• わずかな学力差（θの違い）で正答確率が極端に変化する
  → 「現実の受験行動」ではあまり起こらない
  → データのノイズ、サンプル数不足、異常応答などの影響を受けている可能性あり
  →　「統計的には識別力が高いけど、現実的には怪しい」

2. 異常値としてモデルの安定性を損なうリスクがある

• 極端に高い識別力項目は、モデル推定の不安定要因になる
  
• IRTでは、項目全体のバランスが大事
  
• 突出した項目があると、テスト全体の尺度（θ）推定が歪む可能性あり
  
• θ推定がその項目に引っ張られ、他の項目の情報が正しく活かされない

3. 実務的な理由

• 例えば受験者がほんの少し間違えただけで
  → θ（能力値）が過剰に低く評価されるリスクがある
  
• 逆にたまたま正答しただけで
  → θが過剰に高く評価されるリスクがある
  
• テスト結果の公平性・安定性を守るためには
  
• 異常に高い識別力の項目は除外または別扱い（参考扱い）にするのが一般的

識別力が高すぎる項目の評価
・識別力は、受験者のθ（能力）による正答確率の変化の鋭さを表す
・識別力が高いと　→　能力の違いに敏感に反応する
→　θが少し上がるだけで正答率が大きく上がる

識別力が2以上だと次にようなリスクがある

1. 現実を正しく反映できず
   
2. モデルの安定性を損ない
   
3. 受験者の能力を正しく評価できない

結論：

・Q28とQ15・・・識別力が異常値なので除外

適度な識別力項目 (0.3〜2.0) の ICC

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. 識別力が0.3以上、2.0以下の項目だけ抽出
low_disc_items <- param_df %>%
  filter(Discrimination > 0.3 & Discrimination < 2.0)

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%  # θ = 4 の位置
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +  # ★ ここでy軸を0〜1に固定 ★
  labs(
    title = "適度な識別力項目 (0.3〜2.0) の ICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")  # 凡例を消してスッキリ

適度な識別力項目 (0.3〜2.0) の評価

・θ（学力）が上がるにつれて正答確率がきちんと上昇
・学力がやや低い付近（−2 〜 0）で、正答確率がぐっと上がる項目が多い
・θが低い受験者と高い受験者で正答率に差が出る
→ つまり、能力をよく区別できる優れた項目群

識別力ごとのICC

識別力	評価	詳細
1.0〜2.0	◎ 非常に良い	能力差を鋭く捉えられる（例：Q20〜Q16)

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. 16〜Q20だけ抽出
low_disc_items <- param_df %>%
  filter(Item %in% c("Q20", "Q27", "Q18", "Q2", "Q21", "Q31", "Q4", "Q33", "Q36", "Q9", "Q16"))

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 400)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成（θ=-1付近）
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(abs(theta - (-1)) == min(abs(theta - (-1)))) %>%  # θ = -1に最も近い点を取る
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    vjust = -0.8,  # 少し上にずらして「曲線の上」に配置
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    title = "指定項目 (θ=-1あたりに項目ラベル)", 
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")

識別力	評価	詳細
0.6〜1.0	○ 良い	標準的な能力測定ができる（Q22〜Q30）

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. Q1, Q3, Q7, Q10, Q12だけ抽出
low_disc_items <- param_df %>%
  filter(Item %in% c("Q22", "Q3", "Q7", "Q34", "Q10", "Q17", "Q1", "Q35", "Q29", "Q12", "Q30"))

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    title = "指定項目 (Q23〜Q32) の ICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")

識別力	評価	詳細
0.3〜0.6	△ 許容範囲	使えるが改善を検討（Q14〜Q5）

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. Q5, Q6, Q14, Q19, Q32だけ抽出
low_disc_items <- param_df %>%
  filter(Item %in% c("Q14", "Q19", "Q6", "Q5"))

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    title = "指定項目 (Q14〜Q5) の ICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")

識別力	評価	詳細
0.3以下	× 削除すべき	推定不可(Q8〜Q25)

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（a=識別力, b=困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # aパラメータ
  Difficulty = ex2$est[, 2]       # bパラメータ
)

# 3. Q5, Q6, Q14, Q19, Q32だけ抽出
low_disc_items <- param_df %>%
  filter(Item %in% c("Q8", "Q11", "Q32", "Q13", "Q24", "Q26", "Q23", "Q25"))

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(low_disc_items), function(i) {
  a <- low_disc_items$Discrimination[i]
  b <- low_disc_items$Difficulty[i]
  item_name <- low_disc_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データ結合
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%
  ungroup()

# 8. ICC描画
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    title = "指定項目 (Q8〜Q25) の ICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")

困難度 b の解釈：

適切な困難力 b の大きさ

−3 〜 3

（出典：芝祐順編『項目反応理論』p.34）

• 困難度の値順に可視化する

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出
difficulty_df <- data.frame(
  Item = rownames(ex2$est),
  Difficulty = ex2$est[, 2]
)

# 3. 困難度で並べ替え（昇順）
difficulty_df <- difficulty_df %>%
  arrange(Difficulty) %>%
  mutate(Item = factor(Item, levels = Item))  # 項目順を固定する

# 4. ラベル位置を調整（棒の先に表示するため）
difficulty_df <- difficulty_df %>%
  mutate(label_position = ifelse(Difficulty >= 0, Difficulty + 0.2, Difficulty - 0.2))

# 5. プロット
ggplot(difficulty_df, aes(x = Difficulty, y = Item)) +
  geom_bar(stat = "identity", aes(fill = Difficulty > 0), width = 0.6) +
  geom_text(aes(label = round(Difficulty, 2), x = label_position), 
            size = 3, hjust = ifelse(difficulty_df$Difficulty >= 0, 0, 1)) +
  # ★ ここで赤い点線を引く ★
  geom_vline(xintercept = c(-3, 3), color = "red", linetype = "dashed", size = 0.2) +
  scale_fill_manual(values = c("skyblue", "salmon"), guide = "none") +
  labs(
    title = "困難度（Difficulty）の分布（-3と3に赤い点線）",
    x = "困難度",
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans")

• Q24, Q26, Q11, Q32が難しすぎる問題（困難度 b が 3 以上)
  
• Q23, Q13, Q34, Q29が簡単すぎる問題（困難度 b が -3 以下)

難しすぎる問題 (Q24, Q26, Q11, Q32) の ICC

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（識別力と困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # 識別力 a
  Difficulty = ex2$est[, 2]       # 困難度 b
)

# 3. 必要な項目だけフィルター（Q25, Q28, Q11, Q34）
selected_items <- param_df %>%
  filter(Item %in% c("Q24", "Q26", "Q11", "Q32"))

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(selected_items), function(i) {
  a <- selected_items$Discrimination[i]
  b <- selected_items$Difficulty[i]
  item_name <- selected_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLモデルのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データをまとめる
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成（θ = 4のときに表示）
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%
  ungroup()

# 8. ICCプロット
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    title = "Q24, Q26, Q11, Q32 の ICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")

• 特にQ24とQ26は問題が難し過ぎて、能力の大小によって正答確率はどちらも低い
  
• 識別力も低い

簡単過ぎる問題 (Q13, Q34, Q29) と異常値 (Q23) の ICC

library(irt)
library(ggplot2)
library(dplyr)

# 1. モデル推定（再掲）
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 2. パラメータ抽出（識別力と困難度）
param_df <- data.frame(
  Item = rownames(ex2$est),
  Discrimination = ex2$est[, 1],  # 識別力 a
  Difficulty = ex2$est[, 2]       # 困難度 b
)

# 3. 必要な項目だけフィルター（Q25, Q28, Q11, Q34）
selected_items <- param_df %>%
  filter(Item %in% c("Q23", "Q13", "Q34", "Q29"))

# 4. θの範囲を設定
theta_vals <- seq(-4, 4, length.out = 200)

# 5. 各項目のICCを計算
icc_list <- lapply(1:nrow(selected_items), function(i) {
  a <- selected_items$Discrimination[i]
  b <- selected_items$Difficulty[i]
  item_name <- selected_items$Item[i]
  
  P <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLモデルのICC公式
  data.frame(
    theta = theta_vals,
    Probability = P,
    Item = item_name
  )
})

# 6. データをまとめる
icc_long <- bind_rows(icc_list)

# 7. ラベル位置データ作成（θ = 4のときに表示）
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(theta == max(theta)) %>%
  ungroup()

# 8. ICCプロット
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1.2) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    hjust = -0.1,
    vjust = 0,
    size = 3,
    show.legend = FALSE
  ) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(
    title = "Q23, Q13, Q34, Q29 の ICC",
    x = expression(theta),
    y = "正答確率"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme(legend.position = "none")

• Q23の困難度・・・-19.87（異常値）
  
• Q13は能力が最低の受験生(−4) でも約50％の正答率
  →　能力が平均の受験生（θ= 0）だと60%
  
• Q28の正答率・・・97%
  
• Q34の正答率・・・96%

標準誤差（識別力）の解釈：

標準誤差（識別力）の範囲

0.1 〜 0.4（識別力 a の標準誤差）

（出典：芝祐順編『項目反応理論』p.34）

library(irt)
library(ggplot2)
library(dplyr)

# モデル推定
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 識別力の標準誤差（標準偏差）を取り出す
disc_se <- ex2$se[, 1]  # 2列目が識別力SE


# 項目名を取得
item_names <- rownames(ex2$est)

# データフレームに変換して、小さい順に並べ替え
disc_se_df <- data.frame(
  Item = item_names,
  Difficulty_SE = disc_se
) %>%
  arrange(Difficulty_SE) %>%
  mutate(Item = factor(Item, levels = Item))  # 並び順を反映

# ラベルを棒グラフより少し右にずらすためオフセット作成
disc_se_df <- disc_se_df %>%
  mutate(label_position = Difficulty_SE + 0.02)

# グラフ描画
ggplot(disc_se_df, aes(x = Difficulty_SE, y = Item)) +
  geom_bar(stat = "identity", fill = "darkgreen") +
  geom_text(aes(x = label_position, label = round(Difficulty_SE, 3)), 
            hjust = 0, size = 4) +  # hjust=0で左寄せ
  geom_vline(xintercept = c(0.2, 0.5), color = "red", linetype = "dashed", size = 0.4) +
  labs(
    title = "識別力の標準誤差（SE）の値順に並べた項目", 
    x = "識別力の標準誤差", 
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme_bw(base_family = "HiraKakuProN-W3")

標準誤差（困難度）の結果

標準誤差の規模	項目例	評価
小さい（安定）	Q8〜Q34	安心して使える
大きい（不安定）	Q28, Q15, Q27	注意。改善または除外を検討

標準誤差（困難度）の解釈：

標準誤差（困難度）の範囲

0.2 〜 0.5（困難度 b の標準誤差）

（出典：芝祐順編『項目反応理論』p.34）

library(irt)
library(ggplot2)
library(dplyr)

# モデル推定
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 困難度の標準誤差（標準偏差）を取り出す
disc_se <- ex2$se[, 2]  # 1列目が困難度の標準誤差

# 項目名を取得
item_names <- rownames(ex2$est)

# データフレームに変換して、小さい順に並べ替え
disc_se_df <- data.frame(
  Item = item_names,
  Discrimination_SE = disc_se
) %>%
  arrange(Discrimination_SE) %>%
  mutate(Item = factor(Item, levels = Item))  # 並び順をグラフに反映

# ★ ラベルを棒グラフより少し右にずらすためにオフセットを作成
disc_se_df <- disc_se_df %>%
  mutate(label_position = Discrimination_SE + 0.02)  # 0.02だけ右にずらす

# グラフ描画（横向き棒グラフ：左小 → 右大）
ggplot(disc_se_df, aes(x = Discrimination_SE, y = Item)) +
  geom_bar(stat = "identity", fill = "purple") +
  geom_text(aes(x = label_position, label = round(Discrimination_SE, 3)), 
            hjust = 0, size = 4) +  # ★hjust=0で左寄せ
  # ★ ここで赤い点線を引く ★
  geom_vline(xintercept = c(0.2, 0.5), color = "red", linetype = "dashed", size = 0.4) +
  labs(
    title = "困難度の標準誤差（SE）の値順に並べた項目", 
    x = "困難度の標準誤差", 
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme_bw(base_family = "HiraKakuProN-W3")  # 文字化け対策

• Q24の標準誤差 = 197
  
• Q26の標準誤差 = 134.625

なぜこんな異常が出るか？

• その項目がほぼ全員正答 or 誤答
  

• サンプルサイズが少ない
  

• 受験者の能力分布と項目難易度がずれている

• Q24の正答率・・・21%
  

• Q26の正答率・・・33%

• Q24, Q26, Q23, Q13, Q11, Q32 を除いて、可視化してみる

library(irt)
library(ggplot2)
library(dplyr)

# モデル推定
ex2 <- est(resp = df2[, 2:37], 
           model = "2PL",        
           engine = "ltm")

# 困難度の標準誤差（標準偏差）を取り出す
disc_se <- ex2$se[, 2]  # 2列目が困難度の標準誤差

# 項目名を取得
item_names <- rownames(ex2$est)

# データフレームに変換して、小さい順に並べ替え
disc_se_df <- data.frame(
  Item = item_names,
  Discrimination_SE = disc_se
) %>%
  filter(!(Item %in% c("Q24", "Q26", "Q23", "Q13", "Q11", "Q32"))) %>%   # ★ここで除外
  arrange(Discrimination_SE) %>%
  mutate(Item = factor(Item, levels = Item))        # 並び順をグラフに反映

# ラベルを棒グラフより少し右にずらすためにオフセットを作成
disc_se_df <- disc_se_df %>%
  mutate(label_position = Discrimination_SE + 0.02)

# グラフ描画（横向き棒グラフ：左小 → 右大）
ggplot(disc_se_df, aes(x = Discrimination_SE, y = Item)) +
  geom_bar(stat = "identity", fill = "purple") +
  geom_text(aes(x = label_position, label = round(Discrimination_SE, 3)), 
            hjust = 0, size = 4) +
  geom_vline(xintercept = c(0.2, 0.5), color = "red", linetype = "dashed", size = 0.4) +
  labs(
    title = "困難度の標準誤差の値順に並べた項目（Q24, Q26, Q23, Q13, Q11, Q32を除外）", 
    x = "困難度の標準誤差", 
    y = "項目"
  ) +
  theme_minimal(base_family = "Hiragino Sans") +
  theme_bw(base_family = "HiraKakuProN-W3")

標準誤差（困難度）の結果

標準誤差の規模	項目例	評価
小さい（安定）	Q20〜Q6など	安心して使える
中程度（不安定）	Q27〜Q34	除外を検討
大きい（推定不可）	Q24, Q26, Q23, Q13, Q11, Q32	除外

• 困難度（難易度）の標準誤差（SE）が大きいということは・・・
  ▶ その項目の困難度の推定が不安定である

標準誤差が大きいと次のような問題が起こる

1. 困難度の推定値が信頼できない

• 例えば、困難度 = -1.2 と推定されても
  
• 標準誤差が大きいと「本当は-1.8かもしれないし、-0.6かもしれない」とブレる
  
• 困難度がはっきりしない項目は、受験者の能力を正確に測る力が弱い

2. テスト全体の信頼性が下がる

• ブレブレの項目がテストに多いと
  →　テスト全体で測っている学力θの精度も落ちる

3. IRTモデルのフィット（適合度）が悪くなる

• 困難度が不安定な項目が多いと、モデル全体にとっても悪影響
  
• 項目反応曲線（ICC）がきれいな S 字ではなく変な形になる

分析結果（理想的な基準）

• 項目パラメタにそれぞれ適切な値を適応して項目の採用・不採用を検討する

適切な識別力、困難度、標準誤差
0.3 〜 2.0（識別力）
−3 〜 3（困難度）
0.1 〜 0.4（識別力 a の標準誤差）
0.2 〜 0.5（困難度 b の標準誤差）
（出典：芝祐順編『項目反応理論』p.34）

library(dplyr)

# --- estとseを読み込み（推定済みex2を使用）---
est <- ex2$est
se <- ex2$se

# --- データフレーム作成（列番号付き）---
irt_result <- data.frame(
  Col_No = 1:length(item_names),  # 列番号
  Item = item_names,              # 項目番号
  a = est[, 1],                # 識別力
  b = est[, 2],                # 困難度
  SE_a = se[, 1],              # 識別力の標準誤差
  SE_b = se[, 2]               # 困難度の標準誤差
)

# --- 判断と理由を付与 ---
irt_result <- irt_result %>%
  mutate(
    # 判断
    Judgment = case_when(
      a < 0.3 | a > 2 |
      b < -3 | b > 3 |
      SE_a < 0.1 | SE_a > 0.4 |
      SE_b < 0.2 | SE_b > 0.5 ~ "検討or削除",
      TRUE ~ "問題ない"
    ),
    
    # 理由
    Reason = case_when(
      a < 0.3 ~ "aが低",
      a > 2 ~ "aが高",
      b < -3 ~ "bが低",
      b > 3 ~ "bが高",
      SE_a < 0.1 ~ "SE_aが小",
      SE_a > 0.4 ~ "SE_aが大",
      SE_b < 0.2 ~ "SE_bが小",
      SE_b > 0.5 ~ "SE_bが大",
      TRUE ~ ""
    )
  ) %>%
  mutate(
    a = round(a, 3),
    b = round(b, 3),
    SE_a = round(SE_a, 3),
    SE_b = round(SE_b, 3)
  ) %>%
  arrange(factor(Judgment, levels = c("検討or削除", "問題ない")))

DT::datatable(irt_result)

• 「問題ない」・・・13項目
  
• 「検討or削除」・・・23項目

理想的な基準評価が「問題ない」13項目だけの ICC

library(ltm)  # IRTモデル推定用パッケージ

# 1. 問題ない項目だけ選択
items_ok <- irt_result %>%
  filter(Judgment == "問題ない") %>%
  pull(Item)  # Item名だけ取り出す

# 2. 元データから該当項目だけ抜き出し
# ※ここで元データ（ex2推定時に使ったデータフレーム）が必要です
# 例）df2とする（あなたの環境に合わせて読み替えてね）
df2_selected <- df2[, items_ok]

# 3. 問題ない項目だけでIRT再推定（2PLモデル）
mod_ok <- ltm(df2_selected ~ z1, IRT.param = TRUE)

# 4. ICCを描画
plot(mod_ok, legend = TRUE, cex = 0.6)

• ここに表示された13項目は条件をすべてクリアーた「問題なし」の項目

理想的な基準評価が「検討or削除」23項目だけの ICC

library(ltm)  # IRTモデル推定用パッケージ

# 1. 問題ない項目だけ選択
items_ok <- irt_result %>%
  filter(Judgment == "検討or削除") %>%
  pull(Item)  # Item名だけ取り出す

# 2. 元データから該当項目だけ抜き出し
# ※ここで元データ（ex2推定時に使ったデータフレーム）が必要です
# 例）df2とする（あなたの環境に合わせて読み替えてね）
df2_selected <- df2[, items_ok]

# 3. 問題ない項目だけでIRT再推定（2PLモデル）
mod_ok <- ltm(df2_selected ~ z1, IRT.param = TRUE)

# 4. ICCを描画
plot(mod_ok, legend = TRUE, cex = 0.6)

• ここに表示された23項目は何からの問題を抱えた「検討or削除」の項目
  
• 「検討or削除」の項目をさらに「検討」「削除すべき」に再分類する

分析結果（識別力を重視した基準）

基準値を緩める

• もしこの試験で最も重視したいのが「識別力」なら
  →　識別力だけを基準にして、「検討」「削除すべき」に再分類する

・項目が能力をどれだけ区別できるかを示すパラメータ

適切な識別力 a の大きさ

0.3 〜 2.0

（出典：芝祐順編『項目反応理論』p.34）

判断	基準（識別力）	説明
削除すべき	a ≤ 0.3 or a ≥ 2.5	情報なし or 特定範囲だけ超鋭敏すぎ
検討	0.3 < a < 0.6 or2.0 ≤ a < 2.5	即削除ではない
問題なし	0.6 ≤ a < 2.0	十分な情報を提供する標準的な範囲

library(dplyr)


est <- ex2$est

# --- データフレーム作成（識別力だけ使う） ---
irt_result <- data.frame(
  Item = item_names,
  Disc = est[, 1]
)

# --- Judgment（識別力だけで判断） ---
irt_result <- irt_result %>%
  mutate(
    Judgment = case_when(
      Disc <= 0.3 | Disc >= 2.5 ~ "削除すべき",
      (Disc > 0.3 & Disc < 0.6) | (Disc >= 2.0 & Disc < 2.5) ~ "検討",
      Disc >= 0.6 & Disc < 2.0 ~ "問題なし",
      TRUE ~ NA_character_  # 念のため（どれにも該当しないケース）
    )
  ) %>%
  mutate(
    Disc = round(Disc, 3)  # 小数第3位に丸め
  ) %>%
  arrange(factor(Judgment, levels = c("削除すべき", "検討", "問題なし")))

DT::datatable(irt_result)

• 「問題なし」・・・22項目
  
• 「検討すべき」・・・4項目
  
• 「削除すべき」・・・10項目

「識別力」で「問題なし」22項目の ICC曲線

library(dplyr)
library(ggplot2)
library(tidyr)

# --- 2. 推定結果をまとめる ---
est <- ex2$est  # パラメータ推定結果
irt_result <- data.frame(
  Row_No = 1:length(item_names),
  Item = item_names,
  Disc = est[, 1],   # 識別力
  Diff = est[, 2]    # 困難度
)

# --- 3. Judgment付与 ---
irt_result <- irt_result %>%
  mutate(
    Judgment = case_when(
      Disc <= 0.3 | Disc >= 2.5 ~ "削除すべき",
      (Disc > 0.3 & Disc < 0.6) | (Disc >= 2.0 & Disc < 2.5) ~ "検討",
      Disc >= 0.6 & Disc < 2.0 ~ "問題なし",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    Disc = round(Disc, 3),
    Diff = round(Diff, 3)
  ) %>%
  arrange(factor(Judgment, levels = c("削除すべき", "検討", "問題なし")))

# --- 4. 問題なし項目だけ選ぶ ---
items_to_plot <- irt_result %>%
  filter(Judgment == "問題なし")

# --- 5. θの範囲設定 ---
theta_vals <- seq(-4, 4, length.out = 100)

# --- 6. 各項目ごとにICCを計算 ---
icc_list <- lapply(1:nrow(items_to_plot), function(i) {
  a <- items_to_plot$Disc[i]  # 識別力
  b <- items_to_plot$Diff[i]  # 困難度
  P_theta <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PL ICC式
  data.frame(
    theta = theta_vals,
    Probability = P_theta,
    Item = items_to_plot$Item[i]
  )
})

# --- 7. ICCデータ結合 ---
icc_long <- bind_rows(icc_list)

# --- 8. ラベル位置データ作成（θ=-1付近＋分散）---
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(abs(theta - (-1)) == min(abs(theta - (-1)))) %>%
  ungroup()

# θを分散させる（少しずつずらす）
n_labels <- nrow(label_positions)
theta_offsets <- seq(-0.5, 0.5, length.out = n_labels)  # -0.5〜0.5の範囲で均等にオフセット
label_positions <- label_positions %>%
  arrange(Item) %>%  # 項目番号順（必要に応じて変更）
  mutate(theta = theta + theta_offsets)  # 少しずつθをずらす

# --- 9. ggplotでICC曲線＋ラベル描画 ---
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    vjust = -0.8,     # 曲線の上に配置
    size = 3,
    show.legend = FALSE
  ) +
  labs(
    title = "問題なし項目のICC曲線",
    x = "θ（能力値）",
    y = "正答確率"
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "Hiragino Sans"),  # フォント（環境に合わせて変えてもOK）
    legend.position = "none"  # 凡例を非表示
  )

「識別力」で「検討」4項目の ICC曲線

library(dplyr)
library(ggplot2)
library(tidyr)

# --- 2. 推定結果をまとめる ---
est <- ex2$est  # パラメータ推定結果
irt_result <- data.frame(
  Row_No = 1:length(item_names),
  Item = item_names,
  Disc = est[, 1],   # 識別力
  Diff = est[, 2]    # 困難度
)

# --- 3. Judgment付与 ---
irt_result <- irt_result %>%
  mutate(
    Judgment = case_when(
      Disc <= 0.3 | Disc >= 2.5 ~ "削除すべき",
      (Disc > 0.3 & Disc < 0.6) | (Disc >= 2.0 & Disc < 2.5) ~ "検討",
      Disc >= 0.6 & Disc < 2.0 ~ "問題なし",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    Disc = round(Disc, 3),
    Diff = round(Diff, 3)
  ) %>%
  arrange(factor(Judgment, levels = c("削除すべき", "検討", "問題なし")))

# --- 4. 問題なし項目だけ選ぶ ---
items_to_plot <- irt_result %>%
  filter(Judgment == "検討")

# --- 5. θの範囲設定 ---
theta_vals <- seq(-4, 4, length.out = 100)

# --- 6. 各項目ごとにICCを計算 ---
icc_list <- lapply(1:nrow(items_to_plot), function(i) {
  a <- items_to_plot$Disc[i]  # 識別力
  b <- items_to_plot$Diff[i]  # 困難度
  P_theta <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PL ICC式
  data.frame(
    theta = theta_vals,
    Probability = P_theta,
    Item = items_to_plot$Item[i]
  )
})

# --- 7. ICCデータ結合 ---
icc_long <- bind_rows(icc_list)

# --- 8. ラベル位置データ作成（θ=-1付近＋分散）---
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(abs(theta - (-1)) == min(abs(theta - (-1)))) %>%
  ungroup()

# θを分散させる（少しずつずらす）
n_labels <- nrow(label_positions)
theta_offsets <- seq(-0.5, 0.5, length.out = n_labels)  # -0.5〜0.5の範囲で均等にオフセット
label_positions <- label_positions %>%
  arrange(Item) %>%  # 項目番号順（必要に応じて変更）
  mutate(theta = theta + theta_offsets)  # 少しずつθをずらす

# --- 9. ggplotでICC曲線＋ラベル描画 ---
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    vjust = -0.8,     # 曲線の上に配置
    size = 3,
    show.legend = FALSE
  ) +
  labs(
    title = "検討項目のICC曲線",
    x = "θ（能力値）",
    y = "正答確率"
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "Hiragino Sans"),  # フォント（環境に合わせて変えてもOK）
    legend.position = "none"  # 凡例を非表示
  )

「識別力」で「削除すべき」10項目の ICC曲線

library(dplyr)
library(ggplot2)
library(tidyr)

# --- 2. 推定結果をまとめる ---
est <- ex2$est  # パラメータ推定結果
irt_result <- data.frame(
  Row_No = 1:length(item_names),
  Item = item_names,
  Disc = est[, 1],   # 識別力
  Diff = est[, 2]    # 困難度
)

# --- 3. Judgment付与 ---
irt_result <- irt_result %>%
  mutate(
    Judgment = case_when(
      Disc <= 0.3 | Disc >= 2.5 ~ "削除すべき",
      (Disc > 0.3 & Disc < 0.6) | (Disc >= 2.0 & Disc < 2.5) ~ "検討",
      Disc >= 0.6 & Disc < 2.0 ~ "問題なし",
      TRUE ~ NA_character_
    )
  ) %>%
  mutate(
    Disc = round(Disc, 3),
    Diff = round(Diff, 3)
  ) %>%
  arrange(factor(Judgment, levels = c("削除すべき", "検討", "問題なし")))

# --- 4. 問題なし項目だけ選ぶ ---
items_to_plot <- irt_result %>%
  filter(Judgment == "削除すべき")

# --- 5. θの範囲設定 ---
theta_vals <- seq(-4, 4, length.out = 100)

# --- 6. 各項目ごとにICCを計算 ---
icc_list <- lapply(1:nrow(items_to_plot), function(i) {
  a <- items_to_plot$Disc[i]  # 識別力
  b <- items_to_plot$Diff[i]  # 困難度
  P_theta <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PL ICC式
  data.frame(
    theta = theta_vals,
    Probability = P_theta,
    Item = items_to_plot$Item[i]
  )
})

# --- 7. ICCデータ結合 ---
icc_long <- bind_rows(icc_list)

# --- 8. ラベル位置データ作成（θ=-1付近＋分散）---
label_positions <- icc_long %>%
  group_by(Item) %>%
  filter(abs(theta - (-1)) == min(abs(theta - (-1)))) %>%
  ungroup()

# θを分散させる（少しずつずらす）
n_labels <- nrow(label_positions)
theta_offsets <- seq(-0.5, 0.5, length.out = n_labels)  # -0.5〜0.5の範囲で均等にオフセット
label_positions <- label_positions %>%
  arrange(Item) %>%  # 項目番号順（必要に応じて変更）
  mutate(theta = theta + theta_offsets)  # 少しずつθをずらす

# --- 9. ggplotでICC曲線＋ラベル描画 ---
ggplot(icc_long, aes(x = theta, y = Probability, color = Item)) +
  geom_line(size = 1) +
  geom_text(
    data = label_positions,
    aes(label = Item),
    vjust = -0.8,     # 曲線の上に配置
    size = 3,
    show.legend = FALSE
  ) +
  labs(
    title = "削除すべき項目のICC曲線",
    x = "θ（能力値）",
    y = "正答確率"
  ) +
  theme_minimal() +
  theme(
    text = element_text(family = "Hiragino Sans"),  # フォント（環境に合わせて変えてもOK）
    legend.position = "none"  # 凡例を非表示
  )

8.6 潜在特性値（学力θ）の推定

• IRTモデルでは、観測された回答データ（正答/誤答）を使って
  →　受験者の学力θ（潜在特性値）を推定する
  
• その代表的な方法の一つが 最尤推定法（MLE: Maximum Likelihood Estimation）

最尤推定法の基本的な考え方 ・学力 θ の候補値を −3 から 3 まで 0.1 ごとに調べる
・どの θ のときに「そのパターンの回答が起こる可能性（尤度）」が最も高くなるかを探す
・その「最も可能性が高いθ」が、学力（潜在特性）の推定値

学力θを推定する確率計算

• 学力θを推定する確率計算

• ltmパッケージの中に入っているデータを使う
  • テスト結果を採点

• 受験者数: 220人
  

• 項目数：36項目
  

• 正答なら1、誤答なら０

観測された回答データ（正答/誤答）

• 観測された解答データ (irt_was.csv) を読み取った df2 を使う

DT::datatable(df2)

dim(df2)

[1] 220  37

• theta（学力θのこと）の値が −3 から 3 までの範囲で、0.1 ごとの項目ごとの正答率を計算してみる

library(ltm)

# 2PLモデルによるIRT分析（36項目に限定）
mod <- ltm(df2[, 2:37] ~ z1, IRT.param = TRUE)

# 能力値 θ の範囲指定
theta_vals <- seq(-3, 3, by = 0.1)

# モデルから項目パラメーター（識別力 a、困難度 b）を抽出
coefs <- coef(mod)  # 戻り値は列1 = Dffclt (b), 列2 = Discr (a)

# 正答確率を格納するデータフレームの作成
icc_df <- data.frame(theta = theta_vals)

# 各項目に対して正答確率を計算（2PLモデルの公式）
for (i in 1:nrow(coefs)) {
  b <- coefs[i, 1]  # 困難度 b
  a <- coefs[i, 2]  # 識別力 a
  P_theta <- 1 / (1 + exp(-a * (theta_vals - b)))  # 2PLのICC式
  icc_df[[paste0("Q", i)]] <- round(P_theta, 4)
}

• theta（学力θのこと）の値が −3 から 3 までの範囲で、0.1 ごとの項目ごとの正答率が計算できた

Q1 から Q36 までの正答確率

DT::datatable(icc_df)

• theta（学力θのこと）の値が 0.1 ごとの項目ごとの正答率をすべて確認できる

• この結果を項目特性曲線として可視化してみる

plot(mod, 
  type = "ICC", 
  items = 1:36,
  cex = 0.6) 

# 縦の点線を追加（θ = -3）
abline(v = -3, col = "red", lty = 2, lwd = 1)

• theta（学力θ） = -3 に赤い点線を引いた
  
• theta = -3 の場合、Q25 の数値 0.5306

8.7. テスト情報曲線 (TIC)

テスト情報曲線でわかること

このテストはどの能力レベルを正確に測れているか？

• IRTではテストの推定精度を計算できる
  = 指定した学力θがどの程度の「誤差」を含むかがわかる
  
• テストの測定精度のことを「テスト情報量」とも呼ぶ
  
• テスト情報量は学力θごとに計算できる

テスト情報関数 ・2パラメタ・ロジスティック・モデルの場合のテスト情報量は次の式で表せる

\[I(\theta) = 1.7^2\sum_{j=1}^na_j^2P_j(\theta)Q_j(\theta)\]

変数	内容
\(I(θ)\)	テスト情報量
1.7	定数
\(a_j\)	項目\(j\) の識別力
\(P_j(\theta)\)	学力θに応じた項目 \(j\) に正答する確率
\(Q_j(\theta)\)	学力θに応じた項目 \(j\) に誤答する確率

• IRTにおける措定精度は、個々の学力θに応じて計算できる

• 学力θごとに計算したテスト情報量
  →　テスト情報曲線 (Test Information Curve，TIC) で可視化する

• 「テストがどの能力レベル（θ）でどれくらい正確に測れているか」を可視化できる
  • 潜在特性値の値毎にテスト情報量を計算しプロットする
  

• tif関数を使ってテスト情報量を計算
  →　plot関数で TIC を作成

• 様々な潜在特性値\(θ\)（能力）におけるテスト情報量（Information）が I

• irtoysパッケージの tif()関数を使って計算する

II <- irtoys::tif(ip = ex2$est) # データに対し2PLMを仮定
                           # x: tif関数で推定した結果を指定する引数
plot(x = II)                # ip: テストに含まれる各項目の項目パラメーターを指定する引数 

・横軸・・・潜在特性値\(θ\)（学力）
・縦軸・・・テスト情報量（測定の精度＝標準誤差の逆数）
・実線・・・テスト情報曲線
→　当該潜在特性値におけるテスト情報量をつなぎ合わせたもの

テスト情報曲線 (TIC) でわかること 1. どの能力レベルを正確に測れているか？
・情報量が高いところ → その学力θレベルでテストが精度が高い
→ \(\theta = −2\)付近の情報量が最も高い
→ \(\theta = −2\)付近のレベルでテストが高精度

・情報量が低いところ
→ その学力θレベルではテストの精度が低い
→ \(\theta = 4\)の情報量が最も低い
→ \(\theta = 4\)のレベルでテストが低精度

例：TICがθ = 0の周辺で高ければ
→「平均的な人を測るのに最適なテスト」だといえる

2. テストの設計意図が見える
・TICがどこで山になるかを見ることで、そのテストがどんな対象向けかが分かる
・TICが\(\theta = −2\)付近で山になっている
→　このテストは比較的能力の低い人向け

・TICがどこで山になるかを見ることで
→ そのテストがどんな対象向けかが分かる

TICの山の位置	意味	対象
θ = 0 周辺	平均均的な受験者向け	一般的な学力テスト
θ > 0（右寄り）	高能力者向け	難関資格・上級試験
θ < 0（左寄り）	初級者・低能力者向け	基礎力診断など

・ここではTICの山の位置が左より
→　平均均的な受験者向けだとわかる

3. 信頼性の高さ（精度）もわかる
・情報量が高い＝その範囲の 標準誤差（SE）が小さい
・標準誤差との関係：

\[SE(\theta) = \frac{1}{\sqrt{{I(\theta)}}}\]

・つまり、情報量が大きいと \(\theta\) の推定がブレにくい（= 信頼できる）

4. テスト情報量の基準を設定する：4以上

・例えば、テスト情報量の基準を「4以上」と設定する
→　情報精度の基準を満たしている学力θの範囲がわかる

項目適合度の結果のポイント どの潜在特性値 \(\theta\) の値で、情報量が最大になるか
・潜在特性値が −2 の辺りでテスト情報量が最大
⇒ 潜在特性値\(θ\)（能力）が低い（ −2 の辺り ）受験者の推定精度が最も高い

潜在特性値（学力θ）の推定:mlebme()

• Rで潜在特性値を推定する
  • mlebme関数を利用して潜在特性値を推定
  ⇒ mlebme関数はirtoysパッケージの中に入っている関数

head(mlebme(resp = df2[, 2:37],　# テストデータを指定
  ip = ex2$est,       # データに対し2PLMを仮定
  method = "BM"))           # 最尤推定法 (ML) による潜在特性値の推定を指定

            est       sem  n
[1,] -0.8767052 0.3681238 36
[2,] -0.6930227 0.3818157 36
[3,] -0.8473560 0.3703335 36
[4,] -1.4528580 0.3200674 36
[5,] -0.4407062 0.4012669 36
[6,]  0.7551844 0.5320713 36

• resp: テストデータを指定する引数
  
• ip: テストに含まれる各項目の項目パラメーターを指定する引数
  
• ex2$est と指定
  →　データに対し2PLMを仮定したときの項目パラメーターの推定値が各項目の項目パラメーターとして指定
  
• 項目パラメーターを推定した後に、潜在特性値の推定を行う
  
• method: どの推定法を用いて潜在特性値を推定するか
  
• ML と指定
  →　最尤推定法による潜在特性値の推定を指定
  ⇒ 全問正答/誤答の受験者がいる場合、コれらの受験者の推定値が求まらない
  
• BM と指定
  →　潜在特性値が標準正規分布に従っていることを仮定しているという情報を加味して推定が行えるようになり（事前情報）
  →　全問正答/誤答の受験者に対しても推定値を得ることができる（ベイズ推定）

theta.est <- mlebme(resp = df2[,2:37],
  ip = ex2$est,
  method="BM")

DT::datatable(theta.est)

• １列目：推定値 (est)
  
• ２列目：標準誤差 (sem)
  
• ３列目：解答した項目 (n)

8.8 局所独立性の検討: irf() & cor()

8.8.1 局所独立性とは

• IRTに基づきテストデータや質問紙への回答を分析する
  ⇒ 局所独立性が仮定されている
  • IRTにおける局所独立性の検討
  
• 潜在特性値 \(\theta\) の値を固定したとき、項目相互間で解答が独立に生じている仮定を満たすこと
  
• 「局所独立の仮定」とは「一次元性の仮定」のこと

一次元性の仮定 各項目の正誤が、潜在特性値 \(\theta\) の値の大小によってのみばらつく

• つまり「受験者が正解するかどうかは、受験者の能力だけで決まる」という仮定

• 局所独立性の検討は \(Q_3\)統計量に基づいて行われることが多い
• \(Q_3\)統計量とは、各項目への回答（観測値）からその期待値を引き
→　得られた残差得点間の相関を求めることで得られる
- \(Q_3\)統計量は、各項目への反応（= 観測値） からその期待値（= 項目反応モデルにより計算される正答確率）を引き
→　得られた残差得点間の相関を求めることで得られる統計量
- その絶対値が 0 に近いほど、項目反応間に局所独立性を仮定できる

• たとえば今の場合、item1 の残差得点\(d_1\)は次の式で表せる

\[d_1 = u_1 - \hat{P_1(\theta)}\]

• \(u_1\): Q1 への解答結果（正答なら１、誤答なら０）
  
• \(\hat{P_1(\theta)}\): 項目パラメーターと潜在特性の推定値から計算される正答確率

8.8.2 Rで \(Q_3\) 統計量を計算する方法

Step 1. irf関数を利用して正答確率 ($f) を推定

• irf関数では2PLMを仮定
• ex1$est と指定
→　データに対し2PLMを仮定したときの項目パラメーターの推定値を各項目の項目パラメーターとして指定
• theta.est[, 1] と指定
→　データに対し2PLMを仮定したときの潜在特性の推定値を指定

⇒ 結果はP1として保存
- P1の中はこんな感じ（かなり長い）

P1 <- irtoys::irf(ip = ex2$est, # 項目パラメーターを指定
  x = theta.est[, 1])  # 各受験者の潜在特性値を指定

変数	内容
$x	各受験者の潜在特性値\(\theta\)（能力）
$f	正答確率の推定値
行	受験者 (220名)
列	項目 (Q1〜Q36)

Step 2. 推定された正答確率とテストデータより、残差得点 \(d\) を算出

• \(d_j = u_j - \hat{P_j(\theta)}\) の値を計算して \(d\) として保存 (1 ≦ j ≦36)
  
• P1$fと指定することで，正答確率の推定値が抽出される
  ⇒ テストデータ (df2, 2:37) から正答確率を引いた残差得点を \(d\) に保存

d <- df2[, 2:37] - P1$f # P1$f と指定して正答確率の推定値を得る  

• 220人の受験生の残差得点を \(d\)に関して、最初の2人分の計算結果を表示してみる

head(d, 2)

         Q1         Q2         Q3         Q4         Q5         Q6        Q7
1 0.3764938 -0.8020783  0.1631724 -0.1601448 -0.3495679 -0.5014304 0.2435011
2 0.3413914 -0.8394026 -0.8595188 -0.1913116 -0.3623292  0.4799635 0.2159193
          Q8         Q9        Q10        Q11       Q12        Q13        Q14
1 -0.5156387 -0.3963863 -0.2938277 -0.2133097 0.5882151  0.4192013 -0.1394617
2 -0.5255719  0.5577902  0.6729046 -0.2192435 0.5543485 -0.5846924  0.8477109
         Q15        Q16        Q17       Q18       Q19       Q20        Q21
1  0.1400500  0.2853420 -0.2009591 0.3995311 0.5077143 0.6791613  0.3724295
2 -0.9091997 -0.7505344  0.7731181 0.3395556 0.4853587 0.5949887 -0.6836084
         Q22       Q23        Q24        Q25        Q26        Q27         Q28
1  0.4285537 0.8251662 -0.2088083 -0.3641840 -0.3244871 0.10250561 0.013926729
2 -0.6141168 0.8273212 -0.2098727 -0.3511968 -0.3256762 0.07911754 0.006752633
         Q29        Q30        Q31        Q32        Q33         Q34        Q35
1 0.08805911  0.7458928 -0.5750497 -0.3306807 -0.7512059 -0.93752709 -0.7731963
2 0.07738319 -0.2766025  0.3661200  0.6631874  0.2133131  0.05396073  0.2023193
         Q36
1 -0.7868388
2  0.1822417

受験者1のQ1 に関する残差得点 \(d_{11}\)をチェック ・例えば、受験者1のQ1 に関する残差得点 \(d_1\) は 0.3764938
・受験者1の Q1 への回答 \(u_{ij} = u_{11}\) を確認してみる

df2[1,1]

[1] 1

• 受験者1の Q1 への回答は 1（正答）
  
• 受験者1の Q1 に関する正答確率の推定値 \(\hat{P_{11}(\theta)}\) は 0.3764938 （上の表を参照）
  →　受験者1の残差得点 \(d_{11}\) は

df2[1,1] - 0.3764938 

[1] 0.6235062

Step 3. cor関数を利用して \(Q_3\) 統計量の値を計算

Q3 <- cor(x = d, 
  y = d, 
  use = "pairwise.complete.obs")

Q3

              Q1           Q2           Q3          Q4           Q5
Q1   1.000000000  0.009886907  0.034121554  0.03204650 -0.022418404
Q2   0.009886907  1.000000000  0.317732473 -0.08574863  0.020644169
Q3   0.034121554  0.317732473  1.000000000  0.01146813  0.045591880
Q4   0.032046504 -0.085748632  0.011468127  1.00000000  0.075721918
Q5  -0.022418404  0.020644169  0.045591880  0.07572192  1.000000000
Q6   0.007891653 -0.020275002 -0.017514937  0.08361975  0.171859994
Q7   0.259282965 -0.133855055  0.193419751  0.05437894  0.055906854
Q8   0.023542930 -0.006175132  0.133085698  0.11080293  0.107717436
Q9   0.208310513  0.019541557  0.066953602 -0.05410054  0.034347951
Q10 -0.053825803  0.025219185 -0.050113633 -0.01950422 -0.025465340
Q11  0.104439715  0.059988111 -0.059023864 -0.01287251 -0.115062269
Q12 -0.069880684 -0.121055322  0.015923937  0.02265305 -0.065214653
Q13  0.006596707 -0.067929669  0.065307629  0.04603092 -0.024677731
Q14 -0.123471683  0.038626994 -0.088242065  0.05058035  0.044342004
Q15 -0.187588925 -0.192223406 -0.112931026 -0.03449322  0.047847993
Q16 -0.065964071  0.148826052  0.059084916 -0.10640883  0.048125580
Q17 -0.125988256 -0.138459632 -0.030708951  0.02089458  0.105634259
Q18 -0.117934891 -0.215137661 -0.050055431  0.02008157 -0.034400696
Q19 -0.042349813 -0.074286510 -0.013183026  0.08005636 -0.044317227
Q20  0.089713531 -0.100337431 -0.137220556 -0.05228693 -0.026598320
Q21  0.206744556  0.045103252  0.067822287 -0.14381815 -0.109519482
Q22 -0.224841515 -0.088713793  0.007106703  0.00526165  0.015919306
Q23  0.043509223 -0.070608449  0.035077077 -0.06147393 -0.106388690
Q24  0.012253545  0.011727148  0.001136900  0.05220205 -0.033701332
Q25 -0.184731860  0.079612724 -0.071704436  0.03408949 -0.087477341
Q26 -0.154967373 -0.086031330 -0.062776133  0.06173078  0.093433018
Q27 -0.163197267 -0.137776365  0.002154665 -0.08487991 -0.122926616
Q28 -0.187230268  0.061033239 -0.111605793 -0.03488348 -0.166869359
Q29 -0.068083714 -0.109336493 -0.007862368 -0.07246145 -0.026042571
Q30  0.030106542 -0.201782640 -0.105756970  0.01136671  0.023234319
Q31  0.030706412  0.079325789 -0.217054018 -0.17931599 -0.023370849
Q32  0.023349267 -0.039722366 -0.103998766  0.10300000  0.001096352
Q33 -0.075554059 -0.075920713 -0.087435195 -0.01217830 -0.052734491
Q34 -0.068405742 -0.013609937 -0.048448079 -0.05144967 -0.174253595
Q35 -0.230976571  0.038050329  0.012540627 -0.11322961 -0.032161768
Q36 -0.088632485  0.180440721 -0.233642934 -0.23914696 -0.053835199
              Q6           Q7            Q8           Q9          Q10
Q1   0.007891653  0.259282965  0.0235429299  0.208310513 -0.053825803
Q2  -0.020275002 -0.133855055 -0.0061751319  0.019541557  0.025219185
Q3  -0.017514937  0.193419751  0.1330856982  0.066953602 -0.050113633
Q4   0.083619752  0.054378939  0.1108029337 -0.054100540 -0.019504222
Q5   0.171859994  0.055906854  0.1077174358  0.034347951 -0.025465340
Q6   1.000000000 -0.140294600  0.0699840399  0.087207759 -0.020695027
Q7  -0.140294600  1.000000000  0.0948314628  0.064404871 -0.023663912
Q8   0.069984040  0.094831463  1.0000000000 -0.085170351 -0.117841219
Q9   0.087207759  0.064404871 -0.0851703507  1.000000000 -0.078668991
Q10 -0.020695027 -0.023663912 -0.1178412186 -0.078668991  1.000000000
Q11 -0.164143717 -0.110141493 -0.1548236524 -0.078552088 -0.017609863
Q12 -0.145540548 -0.126604378 -0.0159544118  0.010083255 -0.003034357
Q13  0.018133965  0.117912749 -0.0684955662  0.008986200 -0.060879465
Q14 -0.106378735 -0.066357408 -0.1049457200  0.051559876  0.208205192
Q15 -0.132096698 -0.023591737 -0.0471786107 -0.182609213 -0.078644845
Q16 -0.112140877  0.048761829  0.0315694218 -0.028769065  0.008935241
Q17 -0.108526818  0.006101086 -0.0428148063 -0.061125795  0.004969302
Q18  0.065298568 -0.066404898  0.0180952265 -0.057961006 -0.190419720
Q19  0.051631430  0.090103398 -0.1364626145 -0.042495498 -0.061669890
Q20  0.060383968 -0.071822947  0.0915433737 -0.169671204 -0.145621562
Q21 -0.040293453  0.154152753 -0.0092376612  0.072168602 -0.158168578
Q22  0.099799568 -0.060743930 -0.0811958426 -0.143751454  0.088217877
Q23 -0.081254334  0.068457904 -0.0568037216 -0.062129096  0.046325247
Q24  0.011589912  0.031428589  0.0029310995  0.056204002 -0.130926120
Q25  0.023226209 -0.102309489 -0.0301042299 -0.105153360  0.034539943
Q26  0.075275843 -0.023846631  0.0155331299  0.103550844 -0.103240811
Q27 -0.180994174 -0.051723433 -0.0863186246 -0.143704130  0.003949235
Q28  0.079706784 -0.058429068  0.0471657005  0.030711936  0.056478690
Q29 -0.038908008 -0.116288848  0.0173776712 -0.051279941 -0.151326181
Q30 -0.092357134 -0.032637477  0.0438337661  0.075657369 -0.137503610
Q31 -0.092200886 -0.160488266 -0.1031664510 -0.032744618  0.093132185
Q32  0.009361335 -0.027806697 -0.0225877502  0.028363739  0.071103164
Q33  0.029218601 -0.176797110  0.0002102769 -0.114417709  0.059726742
Q34  0.049427170  0.046074375 -0.0552229930 -0.001893341  0.007064854
Q35 -0.005126877 -0.211961241  0.0205869432 -0.081962235  0.020231377
Q36  0.046336997 -0.120572563 -0.1167839009 -0.186826982  0.030329647
             Q11          Q12          Q13          Q14         Q15
Q1   0.104439715 -0.069880684  0.006596707 -0.123471683 -0.18758892
Q2   0.059988111 -0.121055322 -0.067929669  0.038626994 -0.19222341
Q3  -0.059023864  0.015923937  0.065307629 -0.088242065 -0.11293103
Q4  -0.012872506  0.022653050  0.046030922  0.050580349 -0.03449322
Q5  -0.115062269 -0.065214653 -0.024677731  0.044342004  0.04784799
Q6  -0.164143717 -0.145540548  0.018133965 -0.106378735 -0.13209670
Q7  -0.110141493 -0.126604378  0.117912749 -0.066357408 -0.02359174
Q8  -0.154823652 -0.015954412 -0.068495566 -0.104945720 -0.04717861
Q9  -0.078552088  0.010083255  0.008986200  0.051559876 -0.18260921
Q10 -0.017609863 -0.003034357 -0.060879465  0.208205192 -0.07864484
Q11  1.000000000  0.027320010  0.065748164 -0.094486634  0.06554182
Q12  0.027320010  1.000000000  0.052617014  0.108252984 -0.01863240
Q13  0.065748164  0.052617014  1.000000000 -0.157073025  0.09797507
Q14 -0.094486634  0.108252984 -0.157073025  1.000000000 -0.11071151
Q15  0.065541823 -0.018632396  0.097975071 -0.110711508  1.00000000
Q16  0.093484341 -0.066244049  0.006336439 -0.107227284  0.03668920
Q17 -0.028490224  0.038794226 -0.104642371  0.130970274  0.08648302
Q18 -0.065776876 -0.036146306  0.016786472 -0.146250990  0.19868159
Q19  0.127463148 -0.101394841  0.309542389 -0.245450028  0.02848397
Q20  0.110835939 -0.025467038  0.025532361 -0.155068734 -0.10609646
Q21  0.001575685 -0.077785074 -0.019877492 -0.132610610 -0.03704507
Q22 -0.057145264 -0.084625555 -0.110013501  0.043419251  0.05240642
Q23  0.007281967 -0.104003285 -0.008896153 -0.022232327 -0.10305158
Q24 -0.139887556 -0.154906672  0.040934718 -0.041065175  0.03302357
Q25 -0.033000970  0.097170242  0.124061691 -0.013516025  0.12064576
Q26  0.016757236  0.073185284  0.202544894 -0.070267105  0.04444988
Q27  0.068355839  0.034136586 -0.055848146  0.055358525  0.10029509
Q28 -0.135444034  0.021877107  0.070767023  0.050244680 -0.21210904
Q29  0.012468339  0.108038104 -0.005736872 -0.089265772 -0.05536671
Q30 -0.181588938  0.082122681 -0.009434220 -0.066045810  0.09182519
Q31  0.072655107 -0.043095849 -0.194768117  0.147985843 -0.17054781
Q32  0.085372989 -0.042099013  0.195068549  0.139047718 -0.11628357
Q33 -0.006694599 -0.009169009  0.004001498  0.051388511 -0.04769245
Q34 -0.062391398 -0.103522685  0.052817175 -0.027267785 -0.05604047
Q35 -0.059905855 -0.063281306 -0.001131171  0.088081655 -0.05083122
Q36  0.021177183 -0.048643109 -0.093935058 -0.009293623 -0.02955211
             Q16          Q17          Q18           Q19         Q20
Q1  -0.065964071 -0.125988256 -0.117934891 -4.234981e-02  0.08971353
Q2   0.148826052 -0.138459632 -0.215137661 -7.428651e-02 -0.10033743
Q3   0.059084916 -0.030708951 -0.050055431 -1.318303e-02 -0.13722056
Q4  -0.106408833  0.020894576  0.020081569  8.005636e-02 -0.05228693
Q5   0.048125580  0.105634259 -0.034400696 -4.431723e-02 -0.02659832
Q6  -0.112140877 -0.108526818  0.065298568  5.163143e-02  0.06038397
Q7   0.048761829  0.006101086 -0.066404898  9.010340e-02 -0.07182295
Q8   0.031569422 -0.042814806  0.018095227 -1.364626e-01  0.09154337
Q9  -0.028769065 -0.061125795 -0.057961006 -4.249550e-02 -0.16967120
Q10  0.008935241  0.004969302 -0.190419720 -6.166989e-02 -0.14562156
Q11  0.093484341 -0.028490224 -0.065776876  1.274631e-01  0.11083594
Q12 -0.066244049  0.038794226 -0.036146306 -1.013948e-01 -0.02546704
Q13  0.006336439 -0.104642371  0.016786472  3.095424e-01  0.02553236
Q14 -0.107227284  0.130970274 -0.146250990 -2.454500e-01 -0.15506873
Q15  0.036689196  0.086483024  0.198681588  2.848397e-02 -0.10609646
Q16  1.000000000 -0.029698577  0.145523248 -1.081232e-01 -0.11500640
Q17 -0.029698577  1.000000000  0.035456809 -6.642719e-02 -0.13744892
Q18  0.145523248  0.035456809  1.000000000 -6.741981e-03 -0.07064684
Q19 -0.108123224 -0.066427194 -0.006741981  1.000000e+00  0.03555511
Q20 -0.115006396 -0.137448920 -0.070646841  3.555511e-02  1.00000000
Q21  0.118040723 -0.017578954  0.025189527 -1.138519e-01 -0.07202984
Q22 -0.084077233  0.026103238 -0.052789910  1.129208e-01 -0.09103218
Q23 -0.009491892  0.048660728  0.048793354 -1.044595e-06  0.12882612
Q24 -0.009684228 -0.041791310  0.178026755  1.464077e-02 -0.02636510
Q25 -0.020821781 -0.028339081  0.058300278  2.162773e-01 -0.03639152
Q26 -0.049645990 -0.169120813 -0.047079606  2.150177e-01 -0.03184346
Q27 -0.226380967  0.023552522 -0.074090869 -3.556374e-02 -0.02166301
Q28  0.221335455 -0.048922973  0.023114263  8.122669e-02 -0.06094474
Q29 -0.195272050 -0.083864562 -0.155523458 -7.272298e-02  0.08757880
Q30 -0.061708953 -0.022775827  0.037891617  3.245349e-02  0.04028998
Q31 -0.108526170 -0.032491871 -0.194752854 -3.878954e-02 -0.06453145
Q32 -0.121514450 -0.181096283 -0.110383633  2.090223e-01  0.03828574
Q33 -0.306257349 -0.074984765 -0.117806238  9.977867e-02 -0.02147670
Q34 -0.111903611 -0.043359241 -0.111431600 -3.941760e-02 -0.14672501
Q35 -0.061998877  0.068280768 -0.040908606 -2.429961e-01 -0.06392623
Q36 -0.100979020 -0.148475914 -0.128791389  4.311872e-03 -0.06787787
             Q21          Q22           Q23          Q24         Q25
Q1   0.206744556 -0.224841515  4.350922e-02  0.012253545 -0.18473186
Q2   0.045103252 -0.088713793 -7.060845e-02  0.011727148  0.07961272
Q3   0.067822287  0.007106703  3.507708e-02  0.001136900 -0.07170444
Q4  -0.143818151  0.005261650 -6.147393e-02  0.052202053  0.03408949
Q5  -0.109519482  0.015919306 -1.063887e-01 -0.033701332 -0.08747734
Q6  -0.040293453  0.099799568 -8.125433e-02  0.011589912  0.02322621
Q7   0.154152753 -0.060743930  6.845790e-02  0.031428589 -0.10230949
Q8  -0.009237661 -0.081195843 -5.680372e-02  0.002931100 -0.03010423
Q9   0.072168602 -0.143751454 -6.212910e-02  0.056204002 -0.10515336
Q10 -0.158168578  0.088217877  4.632525e-02 -0.130926120  0.03453994
Q11  0.001575685 -0.057145264  7.281967e-03 -0.139887556 -0.03300097
Q12 -0.077785074 -0.084625555 -1.040033e-01 -0.154906672  0.09717024
Q13 -0.019877492 -0.110013501 -8.896153e-03  0.040934718  0.12406169
Q14 -0.132610610  0.043419251 -2.223233e-02 -0.041065175 -0.01351602
Q15 -0.037045069  0.052406420 -1.030516e-01  0.033023569  0.12064576
Q16  0.118040723 -0.084077233 -9.491892e-03 -0.009684228 -0.02082178
Q17 -0.017578954  0.026103238  4.866073e-02 -0.041791310 -0.02833908
Q18  0.025189527 -0.052789910  4.879335e-02  0.178026755  0.05830028
Q19 -0.113851909  0.112920847 -1.044595e-06  0.014640773  0.21627725
Q20 -0.072029835 -0.091032179  1.288261e-01 -0.026365101 -0.03639152
Q21  1.000000000 -0.007084475  1.291636e-01  0.220622731 -0.20658704
Q22 -0.007084475  1.000000000  1.042582e-01  0.038872276  0.04213179
Q23  0.129163598  0.104258237  1.000000e+00  0.027340762 -0.08885207
Q24  0.220622731  0.038872276  2.734076e-02  1.000000000 -0.09906947
Q25 -0.206587038  0.042131788 -8.885207e-02 -0.099069469  1.00000000
Q26 -0.049832890  0.156483142 -4.340173e-02 -0.040852936  0.27704549
Q27 -0.137378364  0.083337809  6.909193e-02 -0.044156533  0.05776937
Q28 -0.182944964 -0.042875722  8.423069e-02 -0.038877642  0.03732878
Q29 -0.080250133 -0.108991511 -9.687969e-02 -0.032683769  0.01261836
Q30 -0.034790291  0.002390046 -2.732291e-02 -0.006065601 -0.07934241
Q31 -0.189729575 -0.070280614 -6.048095e-02 -0.189601051  0.07145086
Q32 -0.081892040  0.040794625 -5.175579e-02  0.057092141  0.19845010
Q33 -0.258843944 -0.027085811 -4.765537e-02 -0.088544154  0.03801166
Q34 -0.027582299 -0.047692198  3.227390e-02  0.042848359  0.11581128
Q35 -0.116247412 -0.130932779 -5.525026e-02 -0.009563297  0.08013536
Q36  0.069817858  0.011721526 -1.282015e-01 -0.013320761  0.10553556
            Q26           Q27          Q28          Q29          Q30
Q1  -0.15496737 -0.1631972675 -0.187230268 -0.068083714  0.030106542
Q2  -0.08603133 -0.1377763648  0.061033239 -0.109336493 -0.201782640
Q3  -0.06277613  0.0021546646 -0.111605793 -0.007862368 -0.105756970
Q4   0.06173078 -0.0848799054 -0.034883479 -0.072461455  0.011366709
Q5   0.09343302 -0.1229266158 -0.166869359 -0.026042571  0.023234319
Q6   0.07527584 -0.1809941742  0.079706784 -0.038908008 -0.092357134
Q7  -0.02384663 -0.0517234333 -0.058429068 -0.116288848 -0.032637477
Q8   0.01553313 -0.0863186246  0.047165700  0.017377671  0.043833766
Q9   0.10355084 -0.1437041297  0.030711936 -0.051279941  0.075657369
Q10 -0.10324081  0.0039492348  0.056478690 -0.151326181 -0.137503610
Q11  0.01675724  0.0683558387 -0.135444034  0.012468339 -0.181588938
Q12  0.07318528  0.0341365857  0.021877107  0.108038104  0.082122681
Q13  0.20254489 -0.0558481458  0.070767023 -0.005736872 -0.009434220
Q14 -0.07026710  0.0553585252  0.050244680 -0.089265772 -0.066045810
Q15  0.04444988  0.1002950893 -0.212109042 -0.055366713  0.091825185
Q16 -0.04964599 -0.2263809674  0.221335455 -0.195272050 -0.061708953
Q17 -0.16912081  0.0235525218 -0.048922973 -0.083864562 -0.022775827
Q18 -0.04707961 -0.0740908686  0.023114263 -0.155523458  0.037891617
Q19  0.21501774 -0.0355637431  0.081226691 -0.072722979  0.032453493
Q20 -0.03184346 -0.0216630087 -0.060944743  0.087578804  0.040289983
Q21 -0.04983289 -0.1373783637 -0.182944964 -0.080250133 -0.034790291
Q22  0.15648314  0.0833378094 -0.042875722 -0.108991511  0.002390046
Q23 -0.04340173  0.0690919303  0.084230694 -0.096879687 -0.027322912
Q24 -0.04085294 -0.0441565328 -0.038877642 -0.032683769 -0.006065601
Q25  0.27704549  0.0577693669  0.037328776  0.012618358 -0.079342413
Q26  1.00000000 -0.0145508673  0.016309250  0.032833672 -0.190435566
Q27 -0.01455087  1.0000000000 -0.081683223  0.375445898  0.001190998
Q28  0.01630925 -0.0816832231  1.000000000 -0.059032516  0.001690726
Q29  0.03283367  0.3754458979 -0.059032516  1.000000000 -0.033983811
Q30 -0.19043557  0.0011909979  0.001690726 -0.033983811  1.000000000
Q31  0.01677513  0.0009586121 -0.078585833  0.051233901 -0.060512878
Q32  0.28889853 -0.0044039029 -0.086876659 -0.040610335 -0.140414451
Q33  0.05734486  0.0177362381 -0.003594532  0.173366165 -0.062617974
Q34  0.09163023  0.0618966462  0.107842826  0.211232135 -0.019771677
Q35  0.09166825  0.1513400780 -0.004121764  0.181270500 -0.198914612
Q36  0.07734392  0.0192065930 -0.109933507  0.138863822 -0.090929431
              Q31          Q32           Q33          Q34          Q35
Q1   0.0307064124  0.023349267 -0.0755540587 -0.068405742 -0.230976571
Q2   0.0793257887 -0.039722366 -0.0759207128 -0.013609937  0.038050329
Q3  -0.2170540179 -0.103998766 -0.0874351952 -0.048448079  0.012540627
Q4  -0.1793159942  0.102999997 -0.0121783004 -0.051449673 -0.113229607
Q5  -0.0233708486  0.001096352 -0.0527344909 -0.174253595 -0.032161768
Q6  -0.0922008864  0.009361335  0.0292186008  0.049427170 -0.005126877
Q7  -0.1604882657 -0.027806697 -0.1767971096  0.046074375 -0.211961241
Q8  -0.1031664510 -0.022587750  0.0002102769 -0.055222993  0.020586943
Q9  -0.0327446183  0.028363739 -0.1144177091 -0.001893341 -0.081962235
Q10  0.0931321850  0.071103164  0.0597267420  0.007064854  0.020231377
Q11  0.0726551066  0.085372989 -0.0066945992 -0.062391398 -0.059905855
Q12 -0.0430958491 -0.042099013 -0.0091690092 -0.103522685 -0.063281306
Q13 -0.1947681175  0.195068549  0.0040014982  0.052817175 -0.001131171
Q14  0.1479858426  0.139047718  0.0513885109 -0.027267785  0.088081655
Q15 -0.1705478144 -0.116283566 -0.0476924532 -0.056040472 -0.050831216
Q16 -0.1085261701 -0.121514450 -0.3062573494 -0.111903611 -0.061998877
Q17 -0.0324918706 -0.181096283 -0.0749847654 -0.043359241  0.068280768
Q18 -0.1947528543 -0.110383633 -0.1178062376 -0.111431600 -0.040908606
Q19 -0.0387895432  0.209022256  0.0997786738 -0.039417597 -0.242996112
Q20 -0.0645314462  0.038285737 -0.0214767049 -0.146725007 -0.063926230
Q21 -0.1897295750 -0.081892040 -0.2588439443 -0.027582299 -0.116247412
Q22 -0.0702806137  0.040794625 -0.0270858113 -0.047692198 -0.130932779
Q23 -0.0604809526 -0.051755787 -0.0476553743  0.032273904 -0.055250265
Q24 -0.1896010509  0.057092141 -0.0885441542  0.042848359 -0.009563297
Q25  0.0714508552  0.198450101  0.0380116588  0.115811279  0.080135363
Q26  0.0167751285  0.288898533  0.0573448559  0.091630226  0.091668248
Q27  0.0009586121 -0.004403903  0.0177362381  0.061896646  0.151340078
Q28 -0.0785858330 -0.086876659 -0.0035945316  0.107842826 -0.004121764
Q29  0.0512339006 -0.040610335  0.1733661646  0.211232135  0.181270500
Q30 -0.0605128776 -0.140414451 -0.0626179745 -0.019771677 -0.198914612
Q31  1.0000000000  0.204063923  0.1625927249  0.077579164  0.125417309
Q32  0.2040639234  1.000000000  0.0392426162 -0.119599333  0.196978889
Q33  0.1625927249  0.039242616  1.0000000000  0.197027886  0.179135980
Q34  0.0775791637 -0.119599333  0.1970278863  1.000000000  0.129715649
Q35  0.1254173086  0.196978889  0.1791359799  0.129715649  1.000000000
Q36  0.2078304454  0.045166826  0.1064891260  0.138106635  0.190943254
             Q36
Q1  -0.088632485
Q2   0.180440721
Q3  -0.233642934
Q4  -0.239146959
Q5  -0.053835199
Q6   0.046336997
Q7  -0.120572563
Q8  -0.116783901
Q9  -0.186826982
Q10  0.030329647
Q11  0.021177183
Q12 -0.048643109
Q13 -0.093935058
Q14 -0.009293623
Q15 -0.029552114
Q16 -0.100979020
Q17 -0.148475914
Q18 -0.128791389
Q19  0.004311872
Q20 -0.067877872
Q21  0.069817858
Q22  0.011721526
Q23 -0.128201527
Q24 -0.013320761
Q25  0.105535561
Q26  0.077343920
Q27  0.019206593
Q28 -0.109933507
Q29  0.138863822
Q30 -0.090929431
Q31  0.207830445
Q32  0.045166826
Q33  0.106489126
Q34  0.138106635
Q35  0.190943254
Q36  1.000000000

Step 5. \(Q_3\) 統計量の値に関するヒートマップを作成する

• Q3 の絶対値をとり
  
• 対角成分（自己相関1.0）をNAにして除外し
  
• 絶対値が0.2以上の個数を数え、全体の割合を計算してみる

# Q3の絶対値をとる
Q3_abs <- abs(Q3)

# 対角成分（自己相関1.0）をNAにして除外する
diag(Q3_abs) <- NA

# 絶対値が0.2以上の個数を数える
count_0.2_or_more <- sum(Q3_abs >= 0.2, na.rm = TRUE)

# 全体の要素数（対角除外後の数）
total_elements <- sum(!is.na(Q3_abs))

# 割合を計算
proportion_0.2_or_more <- count_0.2_or_more / total_elements

# 結果を表示
cat("絶対値が0.2以上の個数:", count_0.2_or_more, "\n")

絶対値が0.2以上の個数: 66 

cat("割合:", proportion_0.2_or_more, "\n")

割合: 0.05238095 

• 𝑄3の値の絶対値が 0.2以上のケースは66（全体の5%程度）
  
• 得られた結果をヒートマップで可視化してみる

library(ggplot2)
library(reshape2)

# Q3行列を絶対値に
Q3_abs <- abs(Q3)
diag(Q3_abs) <- NA  # 対角成分除外

# meltしてロング形式に
Q3_long <- melt(Q3_abs, varnames = c("Item1", "Item2"), value.name = "Q3_value")

# 0.2超えのフラグをつける
Q3_long$Violation <- Q3_long$Q3_value > 0.2

# ヒートマップ作成
ggplot(Q3_long, aes(x = Item1, y = Item2, fill = Violation)) +
  geom_tile(color = "white") +
  scale_fill_manual(values = c("white", "deeppink")) +
  labs(
    title = "局所独立性違反ヒートマップ（ピンク＝違反）",
    x = "項目",
    y = "項目"
  ) +
  theme_bw(base_family = "HiraKakuProN-W3") +   # 先に白背景＋日本語フォント
  theme(
    axis.text.x = element_text(angle = 90, hjust = 1, vjust = 1),  # ←ここで傾ける
    axis.text.y = element_text(hjust = 1)
  )

• どの項目同士のペアのQ3値が高いのかを一覧にしてみる

library(dplyr)
library(reshape2)

# --- Q3絶対値＋小数点第3位四捨五入
Q3_abs <- abs(Q3)
Q3_abs <- round(Q3_abs, 3)
diag(Q3_abs) <- NA  # 対角成分はNA

# --- ロング形式に変換
Q3_long <- melt(Q3_abs, varnames = c("Item1", "Item2"), value.name = "Q3_abs")

# --- ダブりを除去するため、Item1とItem2を並び替え
Q3_long <- Q3_long %>%
  filter(!is.na(Q3_abs)) %>%
  mutate(Item_min = pmin(as.character(Item1), as.character(Item2)),
         Item_max = pmax(as.character(Item1), as.character(Item2))) %>%
  select(Item1 = Item_min, Item2 = Item_max, Q3_abs)

# --- さらに重複を除去
Q3_long <- distinct(Q3_long)

# --- 0.2以上違反だけ取り出して、大きい順に並べる
violations <- Q3_long %>%
  filter(Q3_abs >= 0.2) %>%
  arrange(desc(Q3_abs))

DT::datatable(violations)

局所独立性の検討のポイント

各項目の正誤が、潜在特性値 \(\theta\) の値の大小によってのみばらつくかどうか
・Q_3 の絶対値が 0 に近いほど、項目反応間に局所独立性を仮定できる

・\(Q_3\) の値の絶対値が 0.2以上の場合　→　問題あり
→　局所依存性の疑い（Chen & Thissen, 1997）

・\(Q_3\) の値の絶対値が 0.2以下の場合　→　問題なし

・\(Q_3\) の値の絶対値の値を計算してみた
→　\(Q_3\) の値の絶対値が0.2以上の個数が 66 　→　問題あり
→　局所独立性が成立していないことを示唆
→　対策が必要

• Q3行列は「IRTモデルで推定された正答確率」と「実際の回答データ」 のズレ（= 残差）から作った残差相関行列
  
• 従って、Q3が0.2以上というのは「局所独立性の仮定が破れている」ということ →　次のいずれかの可能性が高いことを示唆している

🔵 学力θで説明しきれない強い依存関係が残っている
🔵 何か別の理由で項目同士が結びついている

・\(Q_3\)はあくまで一つの基準に過ぎないので注意が必要

References

• Lord, F. M. (1952). A theory of test scores. Psychometric Monograph.
  
• 教育心理学会チュートリアル2018資料
  
• 項目反応理論(IRT)の考え方と実践(2019)
  
• 芝祐順編 (1991)『項目反応理論: 基礎と応用』、東京大学出版会