第14章 Rによる等化
1 パッケージ.
- plinkパッケージを読み込む.
library(lattice)
library(sm)## Package 'sm', version 2.2-5.7: type help(sm) for summary informationlibrary(msm)
library(MASS)##
## Attaching package: 'MASS'## The following object is masked from 'package:sm':
##
## musclelibrary(polycor)
library(ltm)
library(irtoys)
library(plink)##
## Attaching package: 'plink'## The following objects are masked from 'package:ltm':
##
## gpcm, grm2 ２パラメータモデル.
パッケージの関数とその働き.
- as.irt.pars：データの準備.
- plink：等化の実行.
- link.pars：等化後のパラメータを出力.
2.1 データ.
共通項目15個、独自項目5個の２種類の試験フォームTとFを考える.
Tを基準としてFを等化する。
2.2 関数 as.irt.pars.
- x：項目パラメータが収められたリスト.
- common：共通項目が含まれるデータフレーム.
- cat：各項目のカテゴリー数を表現するベクトル.
2.3 項目パラメタの準備
- a：識別力.
- b：困難度.
- c：当て推量での正答確率.
（ 2パタメータモデルでは、c=0 ）
#フォームTの項目パラメタの読み込み
paramT <- read.csv("data/paramT.csv")
head(paramT,5)| a | b | c |
|---|---|---|
| 0.8235682 | -0.7774502 | 0 |
| 1.4799587 | 0.1996423 | 0 |
| 0.9487968 | 0.4171269 | 0 |
| 1.1796335 | -1.0853898 | 0 |
| 0.5040311 | -0.0541100 | 0 |
#フォームFの項目パラメタの読み込み
paramF <- read.csv("data/paramF.csv")
head(paramF,5)| a | b | c |
|---|---|---|
| 0.6936499 | -2.7364145 | 0 |
| 0.8282025 | -0.2048634 | 0 |
| 0.7362474 | -0.9427309 | 0 |
| 1.1313340 | -0.4131648 | 0 |
| 0.8762318 | -0.9030303 | 0 |
2.3.1 xの指定.
項目パラメタ行列を要素としてもつリスト.
pm <- list(paramT,paramF)
str(pm)## List of 2
## $ :'data.frame': 20 obs. of 3 variables:
## ..$ a: num [1:20] 0.824 1.48 0.949 1.18 0.504 ...
## ..$ b: num [1:20] -0.7775 0.1996 0.4171 -1.0854 -0.0541 ...
## ..$ c: int [1:20] 0 0 0 0 0 0 0 0 0 0 ...
## $ :'data.frame': 20 obs. of 3 variables:
## ..$ a: num [1:20] 0.694 0.828 0.736 1.131 0.876 ...
## ..$ b: num [1:20] -2.736 -0.205 -0.943 -0.413 -0.903 ...
## ..$ c: int [1:20] 0 0 0 0 0 0 0 0 0 0 ...この pmを x=pm として as.irt.pars で使用します。
2.3.2 等化係数の推定（手計算）.
第11章「等化」では、共通項目を利用した等化の理論が説明されている。 要点は次のとおり。
＜基本＞（p251参照）.
F社の作成したテスト、T社の作成したテストがある。
集団Aが２つのテストを受けたとき、それぞれのテストでの集団Aの能力パラメータを
母平均： \(\mu_{\theta F},~\mu_{\theta T}\) 母標準偏差：\(\sigma_{\theta F},~\sigma_{\theta T}\) とするとき
\(\dfrac{\theta_{iF}-\mu_{\theta F}}{\sigma_{\theta F}}=\dfrac{\theta_{iT}-\mu_{\theta T}}{\sigma_{\theta T}}\).
\(A=\dfrac{\sigma_{\theta
T}}{\sigma_{\theta F}}~,~~K=\mu_{\theta T}-A\mu_{\theta F}\)
とおくと、 \(\theta_{iT}=A
\theta_{iF}+K\).
（A,Kを等化係数という）.
この式により、Fでの能力得点が試験Tでの得点に換算される。つまり、Tを基準とした変換である。
実際の場面では、等化係数Aを求めるのに.
- 困難度パラメータの標準偏差sigmaを用いる（Mean-Sigma法）.
- 識別力パラメータの平均Meanを用いる（Mean-Mean法）.
などが用いられる。
この方法を用いて、等化係数を手計算してみよう。
共通項目は、Tの６〜20番、Fの1〜15番である。
Mean-Sigma法.
- Aは困難度bの標準偏差の比.
- Kは能力の平均値の一次式.
\(A=\dfrac{\sigma_{b T}}{\sigma_{b F}}~,~~K=\mu_{\theta T}-A\mu_{\theta F}\) \(\theta_{iT}=A \theta_{iF}+K\).
これらを利用する。
A <- sd(paramT$b[6:20])/sd(paramF$b[1:15])
A## [1] 1.013564K <- mean(paramT$b[6:20])-A*mean(paramF$b[1:15])
K## [1] 0.2730375Mean-Mean法.
- Aは識別力aの平均の比.
\(A=\dfrac{\mu_{aF}}{\mu_{aT}}~,~~K=\mu_{\theta T}-A\mu_{\theta F}\) \(\theta_{iT}=A \theta_{iF}+K\).
これらを利用する。
A <- mean(paramF$a[1:15])/mean(paramT$a[6:20])
A## [1] 1.06812K <- mean(paramT$b[6:20])-A*mean(paramF$b[1:15])
K## [1] 0.3093083これは、plinkによる計算（後出）に一致している。
2.3.3 能力パラメタ（事前）pretheta.
推定される能力パラメータはTとFで、似てはいるが異なった値を取っている。
#両フォームの能力パラメタの読み込み
pretheta <- read.csv("data/paramtheta.csv")
str(pretheta)## 'data.frame': 100 obs. of 2 variables:
## $ thetaT: num 2.301 0.351 -0.908 1.316 -0.443 ...
## $ thetaF: num 2.542 0.441 -1.306 1.695 -0.49 ...thetaT <- pretheta[,1]
thetaF <- pretheta[,2]
round(apply(pretheta,2,mean),2)## thetaT thetaF
## 0.05 0.02plot(pretheta[,c(2,1)])
abline(v=0)
abline(h=0)
abline(a = 0, b = 1, col = 2)
両フォームの能力パラメタは、ほぼ近い値を取っていることがわかる。
ヒストグラムで比べてみる。
par(mfrow=c(1,2))
hist(thetaT,prob=T)
hist(thetaF,prob=T)
#能力パラメタのリスト化（能力パラメタの等化が必要な場合）
theta <- list(thetaT,thetaF)
str(theta)## List of 2
## $ : num [1:100] 2.301 0.351 -0.908 1.316 -0.443 ...
## $ : num [1:100] 2.542 0.441 -1.306 1.695 -0.49 ...横軸の能力パラメータは標準化された能力因子であるので、どちらも(-3,3)の範囲になっている。
2.4 共通項目common
JT <- 20 # フォームTの項目数
JF <- 20 # フォームFの項目数
comx<- data.frame("T"=6:20,"F"=1:15)
str(comx)## 'data.frame': 15 obs. of 2 variables:
## $ T: int 6 7 8 9 10 11 12 13 14 15 ...
## $ F: int 1 2 3 4 5 6 7 8 9 10 ...2.5 catの定義
rescat <- list(rep(2,JT),rep(2,JF))
rescat## [[1]]
## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## [[2]]
## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 22.6 plink用オブジェクト.
pmT <- as.poly.mod(n=JT,model="drm",items=1:JT)
pmT## An object of class "poly.mod"
## Slot "model":
## [1] "drm"
##
## Slot "items":
## $drm
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20pmF <- as.poly.mod(n=JF,model="drm",items=1:JF)
pmF## An object of class "poly.mod"
## Slot "model":
## [1] "drm"
##
## Slot "items":
## $drm
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20p.mod <- list(pmT,pmF)
str(p.mod)## List of 2
## $ :Formal class 'poly.mod' [package "plink"] with 2 slots
## .. ..@ model: chr "drm"
## .. ..@ items:List of 1
## .. .. ..$ drm: int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
## $ :Formal class 'poly.mod' [package "plink"] with 2 slots
## .. ..@ model: chr "drm"
## .. ..@ items:List of 1
## .. .. ..$ drm: int [1:20] 1 2 3 4 5 6 7 8 9 10 ...2.7 plinkによる等化
＜基本＞（p251参照）.
F社の作成したテスト、T社の作成したテストがある。
集団Aが２つのテストを受けたとき、それぞれのテストでの集団Aの能力パラメータを
母平均： \(\mu_{\theta F},~\mu_{\theta T}\) 母標準偏差：\(\sigma_{\theta F},~\sigma_{\theta T}\) とするとき
\(\dfrac{\theta_{iF}-\mu_{\theta F}}{\sigma_{\theta F}}=\dfrac{\theta_{iT}-\mu_{\theta T}}{\sigma_{\theta T}}\).
\(A=\dfrac{\sigma_{\theta
T}}{\sigma_{\theta F}}~,~~K=\mu_{\theta T}-A\mu_{\theta F}\)
とおくと、 \(\theta_{iT}=A
\theta_{iF}+K\).
（A,Kを等化係数という）.
p138の３行目によれば、表14.2のデータは.
\(\mu_{\theta F}=0,~\mu_{\theta T}=0.3\) 、 \(\sigma_{\theta F}=1,~\sigma_{\theta T}=1\) であるので.
\(A=\dfrac{\sigma_{\theta T}}{\sigma_{\theta F}}~=1,~~~~~K=\mu_{\theta T}-A\mu_{\theta F}=0.3\).
これは次にplinkにより出力される等化係数とほぼ一致している。
2.8 等化後のパラメータ
res <- as.irt.pars(x=pm,common=comx,cat=rescat,
poly.mod=p.mod)
out <- plink(x=res,rescale="MS",ability=theta,
base.grp=1)
summary(out)## ------- group2/group1* -------
## Linking Constants
##
## A B
## Mean/Mean 1.068120 0.309309
## Mean/Sigma 1.013564 0.273038
## Haebara 1.027976 0.273266
## Stocking-Lord 1.031035 0.285172
##
## Ability Descriptive Statistics
##
## group1 group2
## Mean 0.0504 0.2923
## SD 1.1313 1.1597
## Min -2.5125 -2.4236
## Max 2.3006 2.8497Fによる能力パラメータは、等化により 0.3 近くになっている。
能力パラメータは等化
\(\theta_{iT}=1\cdot\theta_{iF}+0.3\)
により、
Fでの能力得点が試験Tでの得点に換算されていることがわかる。
2.9 項目パラメータ.
link.pars(out)## $group1
## [,1] [,2] [,3]
## [1,] 0.8235682 -0.77745018 0
## [2,] 1.4799587 0.19964226 0
## [3,] 0.9487968 0.41712686 0
## [4,] 1.1796335 -1.08538984 0
## [5,] 0.5040311 -0.05410996 0
## [6,] 0.9490341 -2.13148039 0
## [7,] 0.6907286 0.07438215 0
## [8,] 0.6572211 -0.79298132 0
## [9,] 1.1863945 -0.02735968 0
## [10,] 1.2571493 -0.53645166 0
## [11,] 0.5816125 1.86097028 0
## [12,] 1.3900551 -0.68467007 0
## [13,] 0.6074409 -2.47362082 0
## [14,] 1.1359687 -1.40995586 0
## [15,] 0.6955948 0.68738854 0
## [16,] 0.6181592 0.07021495 0
## [17,] 0.8582456 0.97570347 0
## [18,] 1.0994677 -1.80039250 0
## [19,] 1.0843885 -0.56542585 0
## [20,] 1.3439980 0.74135662 0
##
## $group2
## [,1] [,2] [,3]
## [1,] 0.6843671 -2.50049327 0
## [2,] 0.8171191 0.06539582 0
## [3,] 0.7263946 -0.68248014 0
## [4,] 1.1161940 -0.14573097 0
## [5,] 0.8645056 -0.64224098 0
## [6,] 0.5066901 1.82385878 0
## [7,] 1.4205946 -0.71017019 0
## [8,] 0.6988747 -2.31636182 0
## [9,] 1.1790087 -1.37401911 0
## [10,] 0.7390692 0.68701757 0
## [11,] 0.5305547 0.05919524 0
## [12,] 0.8862056 0.97198435 0
## [13,] 1.5578066 -1.60656375 0
## [14,] 1.5837780 -0.42872692 0
## [15,] 1.6062216 0.78702105 0
## [16,] 0.9429131 2.21582553 0
## [17,] 1.3121671 -0.65073523 0
## [18,] 0.6806354 -0.37173119 0
## [19,] 0.8737066 0.78215204 0
## [20,] 1.8231106 0.83078201 0#フォームFの等化前項目パラメタと等化後項目パラメタの比較
equateparam <- cbind(link.pars(out)[[2]],paramF)
colnames(equateparam) <- c(paste0("等化後",c("a","b","c")),
paste0("等化前",c("a","b","c")))
head(equateparam)| 等化後a | 等化後b | 等化後c | 等化前a | 等化前b | 等化前c |
|---|---|---|---|---|---|
| 0.6843671 | -2.5004933 | 0 | 0.6936499 | -2.7364145 | 0 |
| 0.8171191 | 0.0653958 | 0 | 0.8282025 | -0.2048634 | 0 |
| 0.7263946 | -0.6824801 | 0 | 0.7362474 | -0.9427309 | 0 |
| 1.1161940 | -0.1457310 | 0 | 1.1313340 | -0.4131648 | 0 |
| 0.8645056 | -0.6422410 | 0 | 0.8762318 | -0.9030303 | 0 |
| 0.5066901 | 1.8238588 | 0 | 0.5135628 | 1.5300670 | 0 |
plot(paramF[,2],link.pars(out)[[2]][,2])
abline(v=0)
abline(h=0)
abline(h=0.3,col="lightblue")
困難度パラメータは等化
\(\mu_{b T}=1\cdot\mu_{b F}+0.3\) により変換されていることがわかる。
2.9.1 能力パラメタ.
group1：試験Tの能力 group2：試験Fの能力.
link.ability(out)## $group1
## [1] 2.300581855 0.351470489 -0.908344926 1.316451557 -0.442681422
## [6] -0.670353040 0.173006072 0.475440458 2.296876237 -0.863990424
## [11] 0.452188816 -1.536988031 -0.434656795 0.006139786 0.477823516
## [16] 0.646069873 -1.229310168 1.488837511 0.671549757 0.998304430
## [21] -1.282438326 1.336051534 0.692976884 -2.512515824 0.645898498
## [26] 0.559568706 2.027548371 1.340145992 1.286602838 -1.057109358
## [31] -0.172594134 -0.290922555 -0.913281614 0.664547880 0.969736098
## [36] -0.179154827 -0.692934199 0.490560243 -0.370426214 0.712006618
## [41] -0.564688804 -0.384532764 1.111401239 -0.250923448 -1.131197730
## [46] -1.446580461 1.589429640 0.978381743 -2.320455080 -0.534119038
## [51] 0.093715863 0.849666051 0.652537065 0.090070769 1.088445261
## [56] 0.642246082 0.981971641 0.344282008 1.320813580 -2.448326030
## [61] -0.507813807 0.033640523 2.120044662 1.430538482 -2.426986879
## [66] 1.131930495 -0.209885086 1.230970968 0.809685607 -0.814427905
## [71] -0.524818283 0.609579653 -0.926794239 0.626950948 0.329053549
## [76] 0.867391764 0.552755333 -1.677640439 -1.033272438 -1.461700215
## [81] -0.480708353 0.882404018 -1.269185348 -1.996944032 2.282984519
## [86] -0.376894827 -1.629116927 1.706774755 0.565526268 1.476068256
## [91] -0.527736874 -0.583519488 -1.429734141 0.822872367 -1.041281595
## [96] -0.942136179 0.516953605 -1.396086381 -0.746583520 -0.439422990
##
## $group2
## [1] 2.84972006 0.71958243 -1.05071919 1.99081619 -0.22370024 -0.40716859
## [7] -0.17167033 1.22796262 2.63382974 -0.42506138 0.73573264 -1.23760648
## [13] -0.02956701 0.55049625 1.05715461 0.85323106 -1.05817476 1.21893335
## [19] 0.37586055 1.63974016 -1.00934114 1.47359099 1.05553187 -1.90263867
## [25] 0.55310921 0.49327197 2.34369012 1.38920976 1.38728960 -0.83657243
## [31] 0.43956110 0.10417957 -1.32083048 0.92980799 1.43146661 -0.04899259
## [37] -0.53829139 0.91659804 0.32408504 1.01405783 0.04679022 0.24155201
## [43] 1.51441351 0.26054072 -0.94330551 -0.91058807 1.90921735 1.13680219
## [49] -2.08804515 -0.82139903 0.15773357 1.08792563 0.69742024 0.26121735
## [55] 1.45959631 0.41585248 0.69917338 0.57126408 2.10775854 -1.84270612
## [61] -0.12397921 0.43566216 2.24203705 1.84805309 -1.75435973 1.27304572
## [67] -0.39665505 1.26133321 0.98213635 -0.58616291 -0.33881111 0.74011837
## [73] -0.82760957 0.96778546 0.95261622 1.13601517 0.98766136 -1.49590485
## [79] -0.14362406 -1.62726854 -0.37423960 1.33233530 -1.43304237 -2.42358767
## [85] 2.31492648 -0.36144135 -1.66058919 1.49507520 0.88498842 1.34004486
## [91] -0.23004026 -0.52123006 -0.87744284 1.18890562 -0.89434286 -1.03067350
## [97] 1.21005797 -1.00636923 -0.67415775 0.010464132.9.2 Fの等化前能力パラメタと等化後能力パラメタの比較.
eqtheta <- cbind(link.ability(out)[[2]],thetaF)
colnames(eqtheta) <- c("等化後theta","等化前theta")
head(eqtheta)## 等化後theta 等化前theta
## [1,] 2.8497201 2.5421997
## [2,] 0.7195824 0.4405686
## [3,] -1.0507192 -1.3060420
## [4,] 1.9908162 1.6947901
## [5,] -0.2237002 -0.4900907
## [6,] -0.4071686 -0.6711037ヒストグラムを比べてみる。
par(mfrow=c(1,2))
hist(eqtheta[,2],prob=T,xlim=c(-3,3))
hist(eqtheta[,1],prob=T,xlim=c(-3,3))
par(family = "HiraKakuProN-W3") #日本語フォントの指定
plot(eqtheta[,c(2,1)])
abline(v=0)
abline(h=0)
abline(h=0.3,col="lightblue")
ここでも \(\theta_{iT}=1\cdot\theta_{iF}+0.3\) であり、
等化係数は能力パラメータにおいても共通であることがわかる。
以下、参考までに段階反応モデルのコードを出力する。
詳細はテキスト p339 を参照のこと。
3 段階反応モデル.
# フォームTの項目パラメタの読み込み
paramT <- read.csv("data/pgrmT.csv")
paramT2 <- read.csv("data/pgrmT2.csv")
# フォームFの項目パラメタの読み込み
paramF <- read.csv("data/pgrmF.csv")
paramF2 <- read.csv("data/pgrmF2.csv")
# 項目パラメタのリスト化
pm <- list(paramT,paramF)
pm2 <- list(paramT2,paramF2)
#共通項目の指定
comx <- data.frame("T"=6:20,"F"=1:15)
#項目数の指定
JT <- 20
JF <- 20
#各項目のカテゴリ数の指定
rescat <- list(rep(4,JT),rep(4,JF))
#項目反応モデルの指定
pmT <- as.poly.mod(n=JT,model="grm",items=1:JT)
pmF <- as.poly.mod(n=JF,model="grm",items=1:JF)
#irt.parsオブジェクトの作成
res <- as.irt.pars(x=pm,common=comx,cat=rescat,
poly.mod=list(pmT,pmF),location=FALSE)
out <- plink(res,rescale="MS",base.grp=1)
summary(out)## ------- group2/group1* -------
## Linking Constants
##
## A B
## Mean/Mean 0.999203 0.483078
## Mean/Sigma 1.000900 0.484094
## Haebara 0.998447 0.482014
## Stocking-Lord 1.000406 0.483557link.pars(out)## $group1
## b1 b2 b3
## [1,] 1.1873582 -2.495332 0.59924015 2.356994
## [2,] 0.7664458 -1.906957 -0.54211197 1.863036
## [3,] 1.3390791 -1.533934 0.12388197 1.218443
## [4,] 1.0394187 -1.077727 0.10785853 1.802290
## [5,] 0.9806650 -1.970519 -0.27742341 2.771048
## [6,] 1.3311846 -2.390960 -0.51818141 2.187834
## [7,] 0.8947741 -2.504355 1.14339000 1.591959
## [8,] 0.9069090 -1.947751 1.18527240 2.309181
## [9,] 1.1764619 -2.142846 0.01122222 1.478733
## [10,] 0.9194169 -2.055250 -1.01594680 2.502540
## [11,] 0.6264719 -2.369268 -0.81232591 1.944044
## [12,] 0.8368884 -2.093057 0.28629346 2.037840
## [13,] 0.8115216 -3.416342 0.26796251 2.041544
## [14,] 1.1008019 -2.051720 0.07205164 0.966141
## [15,] 1.1209480 -1.441426 -0.46015352 2.612207
## [16,] 0.9915649 -2.155555 -0.12739010 2.779434
## [17,] 1.1664833 -2.752221 -0.33086454 1.760433
## [18,] 0.8872730 -3.062815 -0.27524680 1.740629
## [19,] 0.7122494 -2.982257 -0.05811460 2.810124
## [20,] 1.2714353 -2.690109 0.06326856 2.666690
##
## $group2
## b1 b2 b3
## [1,] 1.2680879 -2.3875149 -0.46972645 2.1319689
## [2,] 0.8616094 -2.5672297 1.15049198 1.5471137
## [3,] 0.8219173 -2.0204187 1.21455494 2.2993512
## [4,] 1.1898218 -2.1720377 0.02307474 1.4464065
## [5,] 0.9037854 -2.0021411 -1.02459438 2.4790931
## [6,] 0.6343291 -2.3040428 -0.86230841 1.8807042
## [7,] 0.8607964 -2.0908819 0.32865582 2.0698986
## [8,] 0.7622907 -3.4796129 0.20169171 2.0787362
## [9,] 1.0910909 -2.0317450 0.02427113 0.9574358
## [10,] 1.1133905 -1.3111032 -0.41435228 2.6929688
## [11,] 1.0477009 -2.2390460 -0.04392256 2.6943811
## [12,] 1.0982732 -2.6906397 -0.33764168 1.7284376
## [13,] 0.9546348 -3.0264528 -0.30976993 1.8051420
## [14,] 0.7954141 -2.9808572 -0.01542877 2.7603745
## [15,] 1.3262338 -2.7479965 0.07685362 2.7424996
## [16,] 1.3796237 -1.5176059 0.42474507 2.6189078
## [17,] 0.7957041 -2.3716011 0.32292716 1.7418252
## [18,] 1.3325309 -1.9936431 1.26071159 2.9049788
## [19,] 1.3127814 -1.6488565 0.26076194 1.8537470
## [20,] 0.5953500 -0.9832675 -0.17441459 2.8890304