第6章 項目反応理論
6 項目反応理論
6.1 項目反応理論入門
No example.
6.2 Rにおける実装
No example.
6.3 分析準備
## パッケージのロード
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.4 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()## 冊子情報の読み込み
form1 <- read.csv("data/forminfo_F0001.csv")
## テスト得点・採点済みデータ読み込み
## 項目部分のみをselect + all_ofで取り出す（第4章参照）
scored_21g1 <-
read.csv("data/scored_data_2021_grade1.csv") %>%
as_tibble() %>%
select(all_of(form1$item_id)) %>%
print()## # A tibble: 1,000 × 50
## M100025 M100021 M100015 M100074 M100004 M100055 M100059 M100100 M100054
## <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 0 0 0 0 1 1 1 1 1
## 2 0 0 0 0 1 1 0 1 0
## 3 1 1 1 0 1 1 1 1 0
## 4 0 1 1 0 1 1 1 1 1
## 5 0 1 1 0 1 1 1 1 1
## 6 0 1 0 0 1 1 1 1 1
## 7 1 1 0 0 1 1 0 1 1
## 8 0 0 0 0 1 1 1 1 0
## 9 0 1 0 0 1 0 0 1 0
## 10 0 0 1 0 0 1 1 1 1
## # … with 990 more rows, and 41 more variables: M100053 <int>, M100062 <int>,
## # M100073 <int>, M100016 <int>, M100077 <int>, M100098 <int>, M100031 <int>,
## # M100008 <int>, M100079 <int>, M100048 <int>, M100024 <int>, M100044 <int>,
## # M100009 <int>, M100057 <int>, M100002 <int>, M100096 <int>, M100010 <int>,
## # M100067 <int>, M100087 <int>, M100083 <int>, M100056 <int>, M100040 <int>,
## # M100085 <int>, M100099 <int>, M100081 <int>, M100017 <int>, M100006 <int>,
## # M100019 <int>, M100042 <int>, M100027 <int>, M100080 <int>, … write.csv(x = scored_21g1, file ="data/scored_21g1.csv")6.4 項目パラメータの推定
6.4.1 モデルの指定
library(mirt)## Loading required package: stats4## Loading required package: latticefit_rasch <- mirt(data = scored_21g1, model = 1, itemtype = "Rasch",
SE = TRUE, verbose = FALSE)
fit_2pl <- mirt(scored_21g1, 1, "2PL", SE = TRUE, verbose = FALSE)
fit_3pl <- mirt(scored_21g1, 1, "3PL", SE = TRUE, verbose = FALSE)## EM cycles terminated after 500 iterations.6.4.2 モデル選択
## ラッシュモデル vs 2PLモデル
anova(fit_rasch, fit_2pl) ## AIC SABIC HQ BIC logLik X2 df p
## fit_rasch 47573.66 47661.97 47668.79 47823.95 -23735.83
## fit_2pl 47437.95 47611.12 47624.48 47928.72 -23618.97 233.709 49 0## 2PLモデル vs 3PLモデル
anova(fit_2pl, fit_3pl) ## AIC SABIC HQ BIC logLik X2 df p
## fit_2pl 47437.95 47611.12 47624.48 47928.72 -23618.97
## fit_3pl 47507.43 47767.19 47787.23 48243.60 -23603.72 30.514 50 0.9876.4.3 項目パラメータの推定値の取得
## 項目パラメータ + 母集団パラメータ
## 結果は項目ごとのリスト形式（1項目のみ表示）
class(coef(fit_2pl, printSE = FALSE))## [1] "mirt_list" "list"coef(fit_2pl, printSE = FALSE)[[1]] # 推定値 + 信頼区間（デフォルト）## a1 d g u
## par 1.587539 -0.174018138 0 1
## CI_2.5 1.350747 -0.353971678 NA NA
## CI_97.5 1.824332 0.005935402 NA NAcoef(fit_2pl, printSE = TRUE)[[1]] # 推定値 + 標準誤差## a1 d logit(g) logit(u)
## par 1.5875394 -0.17401814 -999 999
## SE 0.1208149 0.09181472 NA NAcoef(fit_2pl, printSE = TRUE, IRTpars = TRUE)[[1]] # 識別力・困難度形式## a b g u
## par 1.5875394 0.10961500 0 1
## SE 0.1208149 0.05790081 NA NAlist.
[[1]]はリストの１番目のデータフレーム等.
## データフレーム形式（行が項目 x パラメータ）
## 識別力・困難度パラメータのみ抽出
est_21g1_df <-
coef(fit_2pl, IRTpars = TRUE, simplify = TRUE) %>%
pluck("items") # pluck: リストの要素を取り出すpurrrパッケージの関数
## 抽出直後はgやuも含まれているのに加え，
## 項目IDが変数でなく行名となっていて使いにくい
head(est_21g1_df, 3)## a b g u
## M100025 1.587539 0.1096150 0 1
## M100021 1.720443 -1.2204766 0 1
## M100015 1.541955 -0.4766847 0 1## 項目IDを明示的に変数として持ち，かつaとbのみを抽出する
ip_est_21g1 <-
coef(fit_2pl, IRTpars = TRUE, simplify = TRUE) %>%
pluck("items") %>%
as_tibble(rownames = "item_id") %>%
select(item_id, a, b) %>%
print()## # A tibble: 50 × 3
## item_id a b
## <chr> <dbl> <dbl>
## 1 M100025 1.59 0.110
## 2 M100021 1.72 -1.22
## 3 M100015 1.54 -0.477
## 4 M100074 1.49 1.13
## 5 M100004 1.98 -2.16
## 6 M100055 1.89 -1.60
## 7 M100059 1.43 -0.296
## 8 M100100 1.95 -2.12
## 9 M100054 1.70 -0.376
## 10 M100053 1.60 -1.03
## # … with 40 more rows6.5 受検者パラメータの推定
## 事後平均（EAP; デフォルト）
theta_21g1 <-
fscores(fit_2pl, method = "EAP", full.scores.SE = TRUE) %>%
as_tibble() %>%
set_names(c("theta_est", "theta_se")) %>% # リストやデータフレームの名前を設定
print() # theta_est: 推定値, theta_se: 標準誤差## # A tibble: 1,000 × 2
## theta_est theta_se
## <dbl> <dbl>
## 1 -0.701 0.225
## 2 -0.862 0.226
## 3 0.313 0.238
## 4 -0.251 0.226
## 5 0.834 0.265
## 6 0.301 0.238
## 7 0.267 0.237
## 8 -1.61 0.251
## 9 -1.34 0.238
## 10 -1.02 0.229
## # … with 990 more rowsstr(theta_21g1)## tibble [1,000 × 2] (S3: tbl_df/tbl/data.frame)
## $ theta_est: num [1:1000] -0.701 -0.862 0.313 -0.251 0.834 ...
## $ theta_se : num [1:1000] 0.225 0.226 0.238 0.226 0.265 ...summary(theta_21g1$theta_est)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.9390508 -0.6703030 0.0223903 -0.0003268 0.6482894 2.5753640library(beeswarm)
boxplot(theta_21g1$theta_est)
beeswarm(theta_21g1$theta_est,pch=16,cex=0.5,add=T)
sd(theta_21g1$theta_est)## [1] 0.96677076.6 IRTのグラフ
6.6.1 項目特性曲線（ICC）
## 前半25項目を項目ごとにプロット
plot(fit_2pl, type = "trace", which.item = 1:25, facet_items = TRUE) 
## 後半25項目を1つの図にプロット
plot(fit_2pl, type = "trace", which.item = 26:50, facet_items = FALSE)
6.6.2 テスト特性曲線（TCC）
## TCC
plot(fit_2pl) # type = "score"がデフォルト
6.6.3 テスト情報曲線（TIC）
plot(fit_2pl, type = "info") # テスト情報量
plot(fit_2pl, type = "SE") # 標準誤差
plot(fit_2pl, type = "infoSE") # テスト情報量 + 標準誤差
6.7 検証的分析
6.7.1 探索的分析と検証的分析
No example.
6.7.2 検証的分析を使った等化の原理
No example.
6.7.3 検証的分析による等化の実行
## 冊子情報の読み込み (Form 4, 2022)
form4 <-
read.csv("data/forminfo_F0004.csv") %>%
as_tibble() %>%
print()## # A tibble: 50 × 6
## item_id form_id order key content_domain cognitive_domain
## <chr> <chr> <int> <int> <chr> <chr>
## 1 M100029 F0004 1 1 Number Applying
## 2 M100007 F0004 2 1 Algebra Knowing
## 3 M100089 F0004 3 5 Number Applying
## 4 M100013 F0004 4 5 Number Knowing
## 5 M100003 F0004 5 5 Algebra Reasoning
## 6 M100055 F0004 6 5 Algebra Knowing
## 7 M100059 F0004 7 5 Algebra Applying
## 8 M100100 F0004 8 1 Statistics Knowing
## 9 M100076 F0004 9 2 Statistics Knowing
## 10 M100053 F0004 10 3 Number Applying
## # … with 40 more rows## テスト得点・採点済みデータ読み込み (2022, 1年生)
## 項目部分のみを取り出す
scored_22g1 <-
read.csv("data/scored_data_2022_grade1.csv") %>%
as_tibble() %>%
select(all_of(form4$item_id)) %>%
print(width = 50)## # A tibble: 1,000 × 50
## M100029 M100007 M100089 M100013 M100003 M100055
## <int> <int> <int> <int> <int> <int>
## 1 1 1 1 1 1 1
## 2 0 1 0 0 0 1
## 3 1 0 1 0 0 1
## 4 1 1 1 0 1 1
## 5 1 1 1 1 0 1
## 6 0 1 1 1 1 1
## 7 1 1 1 1 0 1
## 8 0 0 1 0 0 0
## 9 1 1 1 1 1 1
## 10 0 1 1 1 0 1
## # … with 990 more rows, and 44 more variables:
## # M100059 <int>, M100100 <int>, M100076 <int>,
## # M100053 <int>, M100062 <int>, M100064 <int>,
## # M100016 <int>, M100077 <int>, M100032 <int>,
## # M100031 <int>, M100008 <int>, M100023 <int>,
## # M100048 <int>, M100024 <int>, M100090 <int>,
## # M100009 <int>, M100043 <int>, …## パラメータテーブルの取得（pars = "values"がポイント）
partab <- mirt(scored_22g1, 1, "2PL", pars = "values")
head(partab)# これがパラメータテーブル## group item class name parnum value lbound ubound est prior.type
## 1 all M100029 dich a1 1 0.8510000 -Inf Inf TRUE none
## 2 all M100029 dich d 2 0.8808717 -Inf Inf TRUE none
## 3 all M100029 dich g 3 0.0000000 0 1 FALSE none
## 4 all M100029 dich u 4 1.0000000 0 1 FALSE none
## 5 all M100007 dich a1 5 0.8510000 -Inf Inf TRUE none
## 6 all M100007 dich d 6 0.6106930 -Inf Inf TRUE none
## prior_1 prior_2
## 1 NaN NaN
## 2 NaN NaN
## 3 NaN NaN
## 4 NaN NaN
## 5 NaN NaN
## 6 NaN NaN## 冊子情報の統合
anchor <-
inner_join(form1, form4, suffix = c("_F1", "_F4"), by = "item_id") %>%
as_tibble() %>%
print()## # A tibble: 25 × 11
## item_id form_…¹ order…² key_F1 conte…³ cogni…⁴ form_…⁵ order…⁶ key_F4 conte…⁷
## <chr> <chr> <int> <int> <chr> <chr> <chr> <int> <int> <chr>
## 1 M100055 F0001 6 5 Algebra Knowing F0004 6 5 Algebra
## 2 M100059 F0001 7 5 Algebra Applyi… F0004 7 5 Algebra
## 3 M100100 F0001 8 1 Statis… Knowing F0004 8 1 Statis…
## 4 M100053 F0001 10 3 Number Applyi… F0004 10 3 Number
## 5 M100062 F0001 11 2 Geomet… Applyi… F0004 11 2 Geomet…
## 6 M100016 F0001 13 4 Statis… Knowing F0004 13 4 Statis…
## 7 M100077 F0001 14 2 Number Applyi… F0004 14 2 Number
## 8 M100031 F0001 16 4 Algebra Knowing F0004 16 4 Algebra
## 9 M100008 F0001 17 1 Statis… Applyi… F0004 17 1 Statis…
## 10 M100048 F0001 19 5 Statis… Reason… F0004 19 5 Statis…
## # … with 15 more rows, 1 more variable: cognitive_domain_F4 <chr>, and
## # abbreviated variable names ¹form_id_F1, ²order_F1, ³content_domain_F1,
## # ⁴cognitive_domain_F1, ⁵form_id_F4, ⁶order_F4, ⁷content_domain_F4## 項目パラメータの推定値を結合（前のスクリプトに追記）
anchor <-
inner_join(form1, form4, suffix = c("_F1", "_F4"),
by = c("item_id", "content_domain", "cognitive_domain")) %>%
as_tibble() %>%
inner_join(ip_est_21g1, by = "item_id") %>%
print()## # A tibble: 25 × 11
## item_id form_id…¹ order…² key_F1 conte…³ cogni…⁴ form_…⁵ order…⁶ key_F4 a
## <chr> <chr> <int> <int> <chr> <chr> <chr> <int> <int> <dbl>
## 1 M100055 F0001 6 5 Algebra Knowing F0004 6 5 1.89
## 2 M100059 F0001 7 5 Algebra Applyi… F0004 7 5 1.43
## 3 M100100 F0001 8 1 Statis… Knowing F0004 8 1 1.95
## 4 M100053 F0001 10 3 Number Applyi… F0004 10 3 1.60
## 5 M100062 F0001 11 2 Geomet… Applyi… F0004 11 2 1.49
## 6 M100016 F0001 13 4 Statis… Knowing F0004 13 4 1.61
## 7 M100077 F0001 14 2 Number Applyi… F0004 14 2 1.33
## 8 M100031 F0001 16 4 Algebra Knowing F0004 16 4 1.69
## 9 M100008 F0001 17 1 Statis… Applyi… F0004 17 1 1.76
## 10 M100048 F0001 19 5 Statis… Reason… F0004 19 5 1.13
## # … with 15 more rows, 1 more variable: b <dbl>, and abbreviated variable names
## # ¹form_id_F1, ²order_F1, ³content_domain, ⁴cognitive_domain, ⁵form_id_F4,
## # ⁶order_F4## 識別力・困難度形式から傾き・切片形式へ変換（前のスクリプトに追記）
anchor <-
inner_join(form1, form4, suffix = c("_F1", "_F4"),
by = c("item_id", "content_domain", "cognitive_domain")) %>%
as_tibble() %>%
inner_join(ip_est_21g1, by = "item_id") %>%
mutate(d = - a * b)
select(anchor, a, b, d) # 変換の結果の確認## # A tibble: 25 × 3
## a b d
## <dbl> <dbl> <dbl>
## 1 1.89 -1.60 3.02
## 2 1.43 -0.296 0.423
## 3 1.95 -2.12 4.14
## 4 1.60 -1.03 1.65
## 5 1.49 -0.497 0.742
## 6 1.61 -0.304 0.490
## 7 1.33 0.715 -0.947
## 8 1.69 0.426 -0.718
## 9 1.76 0.399 -0.701
## 10 1.13 -0.138 0.157
## # … with 15 more rows## パラメータテーブル内での項目IDとパラメータ名の組合せ
item_name_partab <- paste0(partab$item, "_", partab$name)
head(item_name_partab, 3)## [1] "M100029_a1" "M100029_d" "M100029_g"## 共通項目の項目IDとパラメータ名の組合せ
item_name_anchor_a <- paste0(anchor$item_id, "_a1") # 識別力
head(item_name_anchor_a)## [1] "M100055_a1" "M100059_a1" "M100100_a1" "M100053_a1" "M100062_a1"
## [6] "M100016_a1"item_name_anchor_d <- paste0(anchor$item_id, "_d") # 切片
head(item_name_anchor_d, 3)## [1] "M100055_d" "M100059_d" "M100100_d"## 固定した識別力パラメータのパラメータテーブル内での位置
idx_a <- match(item_name_anchor_a, item_name_partab)
head(idx_a, 3)## [1] 21 25 29idx_d <- match(item_name_anchor_d, item_name_partab)
head(idx_d, 3)## [1] 22 26 30## 対応する位置のみanchor内の値を代入（上書き）
partab1 <- partab
partab1$value[idx_a] <- anchor$a
partab1$value[idx_d] <- anchor$d
## 対応する位置のみestをFALSEに変更し，固定する
partab1$est[idx_a] <- FALSE
partab1$est[idx_d] <- FALSE## 項目側を固定するので受検者側を開放する
partab1$est[partab1$class == "GroupPars"] <- TRUE## tidyverseを用いたパラメータテーブル修正の例
## anchorをロング型に変換
anchor_long <-
anchor %>%
select(item = item_id, a1 = a, d) %>%
pivot_longer(c(a1, d), values_to = "fixed_value")
## オリジナルのpartabと行順を揃えるために連番seqを一時的に追加
partab2 <-
partab %>%
mutate(seq = 1:n()) %>%
left_join(anchor_long, by = c("item", "name")) %>%
mutate(value = if_else(is.na(fixed_value), value, fixed_value),
est = is.na(fixed_value),
est = if_else(class == "GroupPars", TRUE, est),
est = if_else(name %in% c("g", "u"), FALSE, est)) %>%
arrange(seq) %>% # 行順を合わせる
select(all_of(names(partab))) # 列順を合わせる
str(partab2) # tibbleになっていないことを確認## 'data.frame': 202 obs. of 12 variables:
## $ group : chr "all" "all" "all" "all" ...
## $ item : chr "M100029" "M100029" "M100029" "M100029" ...
## $ class : chr "dich" "dich" "dich" "dich" ...
## $ name : chr "a1" "d" "g" "u" ...
## $ parnum : int 1 2 3 4 5 6 7 8 9 10 ...
## $ value : num 0.851 0.881 0 1 0.851 ...
## $ lbound : num -Inf -Inf 0 0 -Inf ...
## $ ubound : num Inf Inf 1 1 Inf ...
## $ est : logi TRUE TRUE FALSE FALSE TRUE TRUE ...
## $ prior.type: chr "none" "none" "none" "none" ...
## $ prior_1 : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...
## $ prior_2 : num NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN ...## 固定推定の実行
fit_22g1_fixed <- mirt(scored_22g1, 1, pars = partab1, verbose = FALSE)
# fit_22g1_fixed <- mirt(scored_22g1, 1, pars = partab2, verbose = FALSE)
## 項目パラメータの推定値を抽出し，加工する
ip_est_22g1 <-
coef(fit_22g1_fixed, simplify = TRUE, IRTpars = TRUE) %>%
pluck("items") %>%
as_tibble(rownames = "item_id") %>%
select(item_id, a, b) %>%
print()## # A tibble: 50 × 3
## item_id a b
## <chr> <dbl> <dbl>
## 1 M100029 1.30 -0.693
## 2 M100007 1.20 -0.496
## 3 M100089 1.53 -1.83
## 4 M100013 1.02 0.353
## 5 M100003 1.32 0.668
## 6 M100055 1.89 -1.60
## 7 M100059 1.43 -0.296
## 8 M100100 1.95 -2.12
## 9 M100076 1.45 -1.53
## 10 M100053 1.60 -1.03
## # … with 40 more rows## 共通項目の推定値をForm1とForm4でマージする
## inner_joinなので共通項目以外は落とされる
ip_est_fixed_2122 <-
inner_join(ip_est_21g1, ip_est_22g1,
by = "item_id", suffix = c("_F1", "_F4")) %>%
print()## # A tibble: 25 × 5
## item_id a_F1 b_F1 a_F4 b_F4
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 M100055 1.89 -1.60 1.89 -1.60
## 2 M100059 1.43 -0.296 1.43 -0.296
## 3 M100100 1.95 -2.12 1.95 -2.12
## 4 M100053 1.60 -1.03 1.60 -1.03
## 5 M100062 1.49 -0.497 1.49 -0.497
## 6 M100016 1.61 -0.304 1.61 -0.304
## 7 M100077 1.33 0.715 1.33 0.715
## 8 M100031 1.69 0.426 1.69 0.426
## 9 M100008 1.76 0.399 1.76 0.399
## 10 M100048 1.13 -0.138 1.13 -0.138
## # … with 15 more rows## 同じ値かどうかをチェックする
with(ip_est_fixed_2122, all.equal(a_F1, a_F4)) # 識別力## [1] TRUEwith(ip_est_fixed_2122, all.equal(b_F1, b_F4)) # 困難度## [1] TRUE## 固定推定時に母集団パラメータが推定されているかチェック
coef(fit_22g1_fixed)$GroupPars # Mean_1 != 0, COV_11 != 1であることを確認## MEAN_1 COV_11
## par 0.02842774 0.9453686write.csv(ip_est_21g1, "data/item_parameter_estimate_2021_grade1.csv", row.names = FALSE)
write.csv(ip_est_22g1, "data/item_parameter_estimate_2022_grade1.csv", row.names = FALSE)
write.csv(theta_21g1, "data/theta_estimate_2021_grade1.csv", row.names = FALSE)6.8 その他の機能
No example.
6.9 章のまとめ
No example.