1. Dichotomous items
The function DataGeneration can be used in a preparation
step. This function returns a set of artificial data and the true
parameters underlying the data.
Alldata <- DataGeneration(model_D = 2,
N=1000,
nitem_D = 15,
latent_dist = "2NM",
d = 1.664,
sd_ratio = 2,
prob = 0.3)
data <- Alldata$data_D
theta <- Alldata$theta
colnames(data) <- paste0("item", 1:15)
- Analysis (parameter estimation)
Mod1 <- IRTest_Dich(data = data,
model = 2,
latent_dist = "LLS",
h=4)
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 1e-04 on 45th iterations.
#>
#> Model Fit:
#> log-likeli -7662.596
#> deviance 15325.19
#> AIC 15393.19
#> BIC 15560.06
#> HQ 15456.61
#>
#> The Number of Parameters:
#> item 30
#> dist 4
#> total 34
#>
#> The Number of Items: 15
#>
#> The Estimated Latent Distribution:
#> method - LLS
#> ----------------------------------------
#>
#>
#>
#>
#> . . . @ @ .
#> . @ @ @ @ @ . . @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> +---------+---------+---------+---------+
#> -2 -1 0 1 2
### Log-likelihood
logLik(Mod1)
#> [1] -7662.596
### The estimated item parameters
coef(Mod1)
#> a b c
#> item1 0.9836460 1.329453969 0
#> item2 2.2856088 -0.687404304 0
#> item3 1.1690163 -0.215261793 0
#> item4 0.8122027 0.003225493 0
#> item5 1.6372745 -1.189646841 0
#> item6 1.2152175 0.121197281 0
#> item7 1.5656469 0.360962838 0
#> item8 2.5239584 1.182579688 0
#> item9 2.3468154 0.148729175 0
#> item10 1.0642603 -0.894474941 0
#> item11 2.2604198 1.540381124 0
#> item12 1.6180704 -0.263752903 0
#> item13 1.5673423 0.147437141 0
#> item14 1.8951604 -1.107805291 0
#> item15 1.5037221 -0.179279324 0
### Standard errors of the item parameter estimates
coef_se(Mod1)
#> a b c
#> item1 0.09101569 0.11521202 NA
#> item2 0.14467868 0.04196280 NA
#> item3 0.08298967 0.06304486 NA
#> item4 0.07293209 0.08380740 NA
#> item5 0.12363246 0.06673970 NA
#> item6 0.08463917 0.06036465 NA
#> item7 0.10018749 0.05100714 NA
#> item8 0.20149922 0.04677130 NA
#> item9 0.13758907 0.03887234 NA
#> item10 0.08506922 0.08394286 NA
#> item11 0.21591435 0.07116718 NA
#> item12 0.10017918 0.04990388 NA
#> item13 0.09811743 0.05014991 NA
#> item14 0.13792130 0.05618842 NA
#> item15 0.09501085 0.05207528 NA
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "MLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
### Standard errors of ability parameter estimates
plot(fscore$theta, fscore$theta_se)
- The result of latent distribution estimation
plot(Mod1, mapping = aes(colour="Estimated"), linewidth = 1) +
lims(y = c(0, .75))+
geom_line(
mapping=aes(
x=seq(-6,6,length=121),
y=dist2(seq(-6,6,length=121), prob = .3, d = 1.664, sd_ratio = 2),
colour="True"),
linewidth = 1)+
labs(title="The estimated latent density using '2NM'", colour= "Type")+
theme_bw()
p1 <- plot_item(Mod1,1)
p2 <- plot_item(Mod1,4)
p3 <- plot_item(Mod1,8)
p4 <- plot_item(Mod1,10)
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
item_fit(Mod1)
#> stat df p.value
#> item1 11.96160 7 0.1018
#> item2 30.48842 7 0.0001
#> item3 12.97163 7 0.0728
#> item4 11.04379 7 0.1367
#> item5 16.82329 7 0.0186
#> item6 10.73184 7 0.1508
#> item7 15.91906 7 0.0259
#> item8 35.92517 7 0.0000
#> item9 19.94981 7 0.0057
#> item10 16.02558 7 0.0249
#> item11 24.35734 7 0.0010
#> item12 18.96927 7 0.0083
#> item13 18.77397 7 0.0089
#> item14 17.79823 7 0.0129
#> item15 13.57151 7 0.0593
reliability(Mod1)
#> $summed.score.scale
#> $summed.score.scale$test
#> test reliability
#> 0.8550052
#>
#> $summed.score.scale$item
#> item1 item2 item3 item4 item5 item6 item7 item8
#> 0.1412426 0.4667651 0.2394863 0.1366549 0.2791276 0.2522459 0.3378448 0.3719255
#> item9 item10 item11 item12 item13 item14 item15
#> 0.5174065 0.1884686 0.2572114 0.3619343 0.3482940 0.3337077 0.3339890
#>
#>
#> $theta.scale
#> test reliability
#> 0.8404062
- Posterior distributions for the examinees
Each examinee’s posterior distribution is identified in the E-step of
the estimation algorithm (i.e., EM algorithm). Posterior distributions
can be found in Mod1$Pk.
set.seed(1)
selected_examinees <- sample(1:1000,6)
post_sample <-
data.frame(
X = rep(seq(-6,6, length.out=121),6),
prior = rep(Mod1$Ak/(Mod1$quad[2]-Mod1$quad[1]), 6),
posterior = 10*c(t(Mod1$Pk[selected_examinees,])),
ID = rep(paste("examinee", selected_examinees), each=121)
)
ggplot(data=post_sample, mapping=aes(x=X))+
geom_line(mapping=aes(y=posterior, group=ID, color='Posterior'))+
geom_line(mapping=aes(y=prior, group=ID, color='Prior'))+
labs(title="Posterior densities for selected examinees", x=expression(theta), y='density')+
facet_wrap(~ID, ncol=2)+
theme_bw()
- Test information function
ggplot()+
stat_function(
fun = inform_f_test,
args = list(Mod1)
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 1),
mapping = aes(color="Item 1")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 2),
mapping = aes(color="Item 2")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 3),
mapping = aes(color="Item 3")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 4),
mapping = aes(color="Item 4")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 5),
mapping = aes(color="Item 5")
)+
lims(x=c(-6,6))+
labs(title="Test information function", x=expression(theta), y='information')+
theme_bw()
2. Polytomous items
Alldata <- DataGeneration(model_P = "GRM",
categ = rep(c(3,7), each = 7),
N=1000,
nitem_P = 14,
latent_dist = "2NM",
d = 1.664,
sd_ratio = 2,
prob = 0.3)
data <- Alldata$data_P
theta <- Alldata$theta
colnames(data) <- paste0("item", 1:14)
- Analysis (parameter estimation)
Mod1 <- IRTest_Poly(data = data,
model = "GRM",
latent_dist = "KDE")
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 1e-04 on 30th iterations.
#>
#> Model Fit:
#> log-likeli -11867.05
#> deviance 23734.1
#> AIC 23870.1
#> BIC 24203.83
#> HQ 23996.94
#>
#> The Number of Parameters:
#> item 67
#> dist 1
#> total 68
#>
#> The Number of Items: 14
#>
#> The Estimated Latent Distribution:
#> method - KDE
#> ----------------------------------------
#>
#>
#>
#> .
#> . @ @ @ @ @ @ @ @ .
#> . @ @ @ @ @ @ @ @ @ @ @
#> . @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> +---------+---------+---------+---------+
#> -2 -1 0 1 2
### Log-likelihood
logLik(Mod1)
#> [1] -11867.05
### The estimated item parameters
coef(Mod1)
#> a b_1 b_2 b_3 b_4 b_5
#> item1 1.7988485 0.05401123 0.09398055 NA NA NA
#> item2 2.0043067 -0.49151815 0.88695638 NA NA NA
#> item3 0.8364980 -3.18445026 -2.38926385 NA NA NA
#> item4 0.9158959 -1.89519730 -1.11485493 NA NA NA
#> item5 2.0544736 1.21477800 1.86998613 NA NA NA
#> item6 2.2364252 -0.86509330 -0.76042262 NA NA NA
#> item7 2.5872397 1.62484703 1.71397551 NA NA NA
#> item8 2.4386132 -1.31540968 -1.08825163 -0.93434347 -0.5066227 -0.4651258
#> item9 2.0036899 2.34331711 2.53589908 3.16351996 3.4377154 NA
#> item10 2.0581555 -0.76315133 -0.52287005 -0.31563388 -0.2768377 0.3010941
#> item11 1.3617716 -2.04232380 -1.52909952 -1.12964290 -0.6051969 -0.5874733
#> item12 1.0856304 -0.85140919 -0.44231746 0.05233322 0.1786312 0.2497041
#> item13 0.9913453 -0.28218061 0.24473148 0.44191560 0.6480541 0.9343655
#> item14 2.1127859 -0.17731670 0.28797638 0.29842595 0.3696553 0.5204740
#> b_6
#> item1 NA
#> item2 NA
#> item3 NA
#> item4 NA
#> item5 NA
#> item6 NA
#> item7 NA
#> item8 -0.3070251
#> item9 NA
#> item10 0.4299311
#> item11 NA
#> item12 1.0180891
#> item13 1.4177310
#> item14 0.6766842
### Standard errors of the item parameter estimates
coef_se(Mod1)
#> a b_1 b_2 b_3 b_4 b_5
#> item1 0.11073114 0.04486570 0.04497786 NA NA NA
#> item2 0.09858216 0.04347355 0.04920590 NA NA NA
#> item3 0.10469227 0.35380581 0.25735550 NA NA NA
#> item4 0.08207781 0.15972272 0.10612244 NA NA NA
#> item5 0.14580557 0.05639000 0.08780972 NA NA NA
#> item6 0.14441077 0.04469104 0.04257786 NA NA NA
#> item7 0.23648116 0.06790597 0.07362468 NA NA NA
#> item8 0.12154940 0.05036114 0.04419347 0.04113448 0.03667716 0.03650767
#> item9 0.24603699 0.15352565 0.17566410 0.26414336 0.31559599 NA
#> item10 0.09610058 0.04531634 0.04187806 0.04018262 0.03998476 0.04105944
#> item11 0.08788259 0.11824810 0.08970124 0.07247437 0.05788855 0.05758618
#> item12 0.06830391 0.07887899 0.06808232 0.06406234 0.06478591 0.06550493
#> item13 0.06828836 0.07211232 0.07048538 0.07355080 0.07864225 0.08840144
#> item14 0.10387197 0.04007432 0.03968341 0.03974503 0.04025363 0.04184932
#> b_6
#> item1 NA
#> item2 NA
#> item3 NA
#> item4 NA
#> item5 NA
#> item6 NA
#> item7 NA
#> item8 0.03634449
#> item9 NA
#> item10 0.04246064
#> item11 NA
#> item12 0.08534119
#> item13 0.11032712
#> item14 0.04428654
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "MLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
### Standard errors of ability parameter estimates
plot(fscore$theta, fscore$theta_se)
- The result of latent distribution estimation
plot(Mod1, mapping = aes(colour="Estimated"), linewidth = 1) +
stat_function(
fun = dist2,
args = list(prob = .3, d = 1.664, sd_ratio = 2),
mapping = aes(colour = "True"),
linewidth = 1) +
lims(y = c(0, .75)) +
labs(title="The estimated latent density using '2NM'", colour= "Type")+
theme_bw()
p1 <- plot_item(Mod1,1)
p2 <- plot_item(Mod1,4)
p3 <- plot_item(Mod1,8)
p4 <- plot_item(Mod1,10)
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
item_fit(Mod1)
#> stat df p.value
#> item1 14.35691 15 0.4987
#> item2 18.95744 15 0.2157
#> item3 13.47278 15 0.5658
#> item4 18.21689 15 0.2514
#> item5 22.98208 15 0.0845
#> item6 23.89178 15 0.0670
#> item7 11.01423 15 0.7516
#> item8 51.99993 47 0.2855
#> item9 20.19194 31 0.9316
#> item10 47.39453 47 0.4565
#> item11 38.71699 39 0.4827
#> item12 44.63516 47 0.5710
#> item13 51.69728 47 0.2955
#> item14 45.54566 47 0.5329
reliability(Mod1)
#> $summed.score.scale
#> $summed.score.scale$test
#> test reliability
#> 0.8423786
#>
#> $summed.score.scale$item
#> item1 item2 item3 item4 item5 item6 item7
#> 0.38547011 0.49427544 0.07314142 0.13961946 0.36918780 0.42828060 0.36663417
#> item8 item9 item10 item11 item12 item13 item14
#> 0.54822704 0.17552217 0.51886901 0.28678300 0.24874696 0.21063051 0.49749765
#>
#>
#> $theta.scale
#> test reliability
#> 0.8699627
- Posterior distributions for the examinees
Each examinee’s posterior distribution is identified in the E-step of
the estimation algorithm (i.e., EM algorithm). Posterior distributions
can be found in Mod1$Pk.
set.seed(1)
selected_examinees <- sample(1:1000,6)
post_sample <-
data.frame(
X = rep(seq(-6,6, length.out=121),6),
prior = rep(Mod1$Ak/(Mod1$quad[2]-Mod1$quad[1]), 6),
posterior = 10*c(t(Mod1$Pk[selected_examinees,])),
ID = rep(paste("examinee", selected_examinees), each=121)
)
ggplot(data=post_sample, mapping=aes(x=X))+
geom_line(mapping=aes(y=posterior, group=ID, color='Posterior'))+
geom_line(mapping=aes(y=prior, group=ID, color='Prior'))+
labs(title="Posterior densities for selected examinees", x=expression(theta), y='density')+
facet_wrap(~ID, ncol=2)+
theme_bw()
- Test information function
ggplot()+
stat_function(
fun = inform_f_test,
args = list(Mod1)
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 1),
mapping = aes(color="Item 1 (3 cats)")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 2),
mapping = aes(color="Item 2 (3 cats)")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 3),
mapping = aes(color="Item 3 (3 cats)")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 8),
mapping = aes(color="Item 8 (7 cats)")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 9),
mapping = aes(color="Item 9 (7 cats)")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 10, "p"),
mapping = aes(color="Item10 (7 cats)")
)+
lims(x=c(-6,6))+
labs(title="Test information function", x=expression(theta), y='information')+
theme_bw()
3. Continuous items
- Statistical details of the continuous IRT
Beta distribution (click)
\[
\begin{align}
f(x) &= \frac{1}{Beta(\alpha, \beta)}x^{\alpha-1}(1-x)^{(\beta-1)}
\\
&=
\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{(\beta-1)}
\end{align}
\]
\(E(x)=\frac{\alpha}{\alpha+\beta}\)
and \(Var(x)=\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta=1)}\)
If we reparameterize \(\mu=\frac{\alpha}{\alpha+\beta}\) and \(\nu=\alpha+\beta\),
\[
f(x) =
\frac{\Gamma(\nu)}{\Gamma(\mu\nu)\Gamma(\nu(1-\mu))}x^{\mu\nu-1}(1-x)^{(\nu(1-\mu)-1)}
\] No Jacobian transformation required since \(\mu\) and \(\nu\) are parameters of the \(f(x)\), not variables.
Useful equations (click)
\(\psi(\bullet)\) and \(\psi_1(\bullet)\) denote for digamma and
trigamma functions, respectively.
\[
\begin{align}
E[\log{x}] &= \int_{0}^{1}{\log{x}f(x) \,dx} \\
&= \int_{0}^{1}{\log{x} \frac{1}{Beta(\alpha,
\beta)}x^{\alpha-1}(1-x)^{(\beta-1)} \,dx} \\
&= \frac{1}{Beta(\alpha, \beta)} \int_{0}^{1}{\log{(x)}
x^{\alpha-1}(1-x)^{(\beta-1)} \,dx} \\
&= \frac{1}{Beta(\alpha, \beta)} \int_{0}^{1}{\frac{\partial
x^{\alpha-1}(1-x)^{(\beta-1)}}{\partial \alpha} \,dx} \\
&= \frac{1}{Beta(\alpha, \beta)} \frac{\partial}{\partial
\alpha}\int_{0}^{1}{ x^{\alpha-1}(1-x)^{(\beta-1)} \,dx} \\
&= \frac{1}{Beta(\alpha, \beta)} \frac{\partial Beta(\alpha,
\beta)}{\partial \alpha} \\
&= \frac{\partial \log{[Beta(\alpha, \beta)]}}{\partial \alpha} \\
&= \frac{\partial \log{[\Gamma(\alpha)]}}{\partial \alpha} -
\frac{\partial \log{[\Gamma(\alpha + \beta)]}}{\partial \alpha} \\
&= \psi(\alpha) - \psi(\alpha+\beta)
\end{align}
\]
Similarly, \(E[\log{(1-x)}]=\psi(\beta) -
\psi(\alpha+\beta)\).
Furthermore, using \(\frac{\partial
Beta(\alpha,\beta)}{\partial \alpha} =
Beta(\alpha,\beta)\left(\psi(\alpha)-\psi(\alpha+\beta)\right)\)
and \(\frac{\partial^2
Beta(\alpha,\beta)}{\partial \alpha^2} =
Beta(\alpha,\beta)\left(\psi(\alpha)-\psi(\alpha+\beta)\right)^2 +
Beta(\alpha,\beta)\left(\psi_1(\alpha)-\psi_1(\alpha+\beta)\right)\),
\[
\begin{align}
E\left[(\log{x})^2\right] &= \frac{1}{Beta(\alpha, \beta)}
\frac{\partial^2 Beta(\alpha, \beta)}{\partial \alpha^2} \\
&= \left(\psi(\alpha)-\psi(\alpha+\beta)\right)^2 +
\left(\psi_1(\alpha)-\psi_1(\alpha+\beta)\right)
\end{align}
\]
This leads to,
\[
\begin{align}
Var\left[\log{x}\right] &= E\left[(\log{x})^2\right] -
E\left[\log{x}\right]^2 \\
&=\psi_1(\alpha)-\psi_1(\alpha+\beta)
\end{align}
\]
Continuous IRT (click)
\[
\mu = \frac{e^{a(\theta -b)}}{1+e^{a(\theta -b)}} \\
\frac{\partial \mu}{\partial \theta} = a\mu(1-\mu) \\
\frac{\partial \mu}{\partial a} = (\theta - b)\mu(1-\mu) \\
\frac{\partial \mu}{\partial b} = -a\mu(1-\mu) \\
\frac{\partial \mu}{\partial \nu} = 0
\]
- Probability of a response
\[
f(x)=P(x|\, \theta, a, b, \nu) =
\frac{\Gamma(\nu)}{\Gamma(\mu\nu)\Gamma(\nu(1-\mu))} x^{\mu\nu-1}
(1-x)^{\nu(1-\mu)-1} \\
\]
\[
\log{f} =
\log{\Gamma(\nu)}-\log{\Gamma(\mu\nu)}-\log{\Gamma(\nu(1-\mu))} +
(\mu\nu-1)\log{x} + (\nu(1-\mu)-1) \log{(1-x)}
\]
\[
\frac{\partial \log{f}}{\partial \theta} =
a\nu\mu(1-\mu)\left[-\psi{(\mu\nu)}+\psi{(\nu(1-\mu))}+
\log{\left(\frac{x}{1-x}\right)}\right]
\]
\[
E\left[ \left(\frac{\partial \log{f}}{\partial \theta}\right)^2 \right]
=
(a\nu\mu(1-\mu))^2\left[E\left[
\log{\left(\frac{x}{1-x}\right)^2}\right]
-2 \left(\psi{(\mu\nu)}-\psi{(\nu(1-\mu))}\right )E\left[
\log{\left(\frac{x}{1-x}\right)}\right]
+\left(\psi{(\mu\nu)}-\psi{(\nu(1-\mu))}\right )^2
\right]
\]
\[
\begin{align}
E\left[ \log{\left(\frac{x}{1-x}\right)^2}\right]
&= E\left[ \log{\left(x\right)^2}\right]
-2 E\left[ \log{\left(x\right)}\log{\left(1-x\right)}\right]
+ E\left[ \log{\left(1-x\right)^2}\right] \\
&= Var\left[ \log{\left(x\right)}\right]+E\left[
\log{\left(x\right)}\right]^2 \\
&\qquad -2 Cov\left[
\log{\left(x\right)}\log{\left(1-x\right)}\right]-2E\left[
\log{\left(x\right)}\right]E\left[ \log{\left(1-x\right)}\right] \\
&\qquad + Var\left[ \log{\left(1-x\right)}\right]+E\left[
\log{\left(1-x\right)}\right]^2 \\
&= \psi_{1}(\alpha)-\psi_{1}(\alpha+\beta) +E\left[
\log{\left(x\right)}\right]^2 \\
&\qquad +2 \psi_{1}(\alpha+\beta)-2E\left[
\log{\left(x\right)}\right]E\left[ \log{\left(1-x\right)}\right] \\
&\qquad + \psi_{1}(\beta)-\psi_{1}(\alpha+\beta)+E\left[
\log{\left(1-x\right)}\right]^2 \\
&= \psi_{1}(\alpha)-\psi_{1}(\alpha+\beta) +\left[
\psi(\alpha)-\psi(\alpha+\beta)\right]^2 \\
&\qquad +2 \psi_{1}(\alpha+\beta)-2
\left(\psi(\alpha)-\psi(\alpha+\beta)\right)\left(\psi(\beta)-\psi(\alpha+\beta)\right) \\
&\qquad +
\psi_{1}(\beta)-\psi_{1}(\alpha+\beta)+\left[\psi(\beta)-\psi(\alpha+\beta)\right]^2
\\
&= \psi_{1}(\alpha) +\psi_{1}(\beta)
+\left[\psi(\alpha)-\psi(\beta)\right]^2
\end{align}
\]
\[
\begin{align}
E\left[\left(\frac{\partial \log{f}}{\partial \theta}\right)^2 \right]
& =
(a\nu\mu(1-\mu))^2\left[E\left[
\log{\left(\frac{x}{1-x}\right)^2}\right]
-2 \left(\psi{(\mu\nu)}-\psi{(\nu(1-\mu))}\right )E\left[
\log{\left(\frac{x}{1-x}\right)}\right]
+\left(\psi{(\mu\nu)}-\psi{(\nu(1-\mu))}\right )^2
\right] \\
&= (a\nu\mu(1-\mu))^2\left[\psi_{1}(\alpha) +\psi_{1}(\beta)
+\left[\psi(\alpha)-\psi(\beta)\right]^2
-2 \left(\psi{(\alpha)}-\psi{(\beta)}\right
)\left(\psi{(\alpha)}-\psi{(\beta)}\right )
+\left(\psi{(\alpha)}-\psi{(\beta)}\right )^2
\right] \\
&= (a\nu\mu(1-\mu))^2\left[\psi_{1}(\alpha) +\psi_{1}(\beta)
\right] \\
\end{align}
\]
\[
I(\theta) = E\left[\left(\frac{\partial \log{f}}{\partial
\theta}\right)^2 \right] = (a\nu\mu(1-\mu))^2\left[\psi_{1}(\alpha)
+\psi_{1}(\beta)\right]
\]
- Log-likelihood in the M-step of the MML-EM procedure
Marginal log-likelihood of an item can be expressed as follows:
\[
\ell = \sum_{j} \sum_{q}\gamma_{jq}\log{L_{jq}},
\]
where \(\gamma_{jq}=E\left[\Pr\left(\theta_j \in
\theta_{q}^{*}\right)\right]\) is the expected probability of
respondent \(j\)’s ability (\(\theta_j\)) belonging to the \(\theta_{q}^{*}\) of the quadrature scheme
and is calculated at the E-step of the MML-EM procedure, and \(L_{jq}\) is the likelihood of respondent
\(j\)’s response at \(\theta_{q}^{*}\) for the item of current
interest.
- First derivative of the likelihood
\[
\frac{\partial \ell}{\partial a} = \sum_{q}
\left(\theta_{q}-b\right)\nu\mu_{q}\left(1-\mu_{q}\right)\left[
S_{1q}-S_{2q}-f_q\left[ \psi(\mu_{q}\nu)-\psi(\nu(1-\mu_{q})) \right]
\right] \\
\frac{\partial \ell}{\partial b} =
-a\sum_{q}\nu\mu_{q}\left(1-\mu_{q}\right)\left[
S_{1q}-S_{2q}-f_q\left[ \psi(\mu_{q}\nu)-\psi(\nu(1-\mu_{q})) \right]
\right] \\
\frac{\partial \ell}{\partial \nu} = N\psi(\nu) +\sum_{q}\left[
\mu_{q}(S_{1q}-f_q\psi(\mu_{q}\nu)) +
(1-\mu_{q})(S_{2q}-f_q\psi(\nu(1-\mu_{q})))
\right]
\]
where \(S_{1q} =
\sum_{j}{\gamma_{jq}\log{x_j}}\) and \(S_{2q} =
\sum_{j}{\gamma_{jq}\log{(1-x_j)}}\). Since \(E_q[S_{1q}]=f_q\left[\psi(\mu_{q}\nu))-\psi(\nu)\right]\)
and \(E_q[S_{2q}]=f_q\left[\psi(\nu(1-\mu_{q})))-\psi(\nu)\right]\),
the expected values of the first derivatives are 0.
To keep \(\nu\) positive, let \(\nu = \exp{\xi}\); \(\frac{\partial\nu}{\partial\xi}=\exp{\xi}=\nu\).
\[
\frac{\partial \ell}{\partial \xi} = N\nu\psi(\nu) +\nu\sum_{q}\left[
\mu_{q}(S_{1q}-f_q\psi(\mu_{q}\nu)) +
(1-\mu_{q})(S_{2q}-f_q\psi(\nu(1-\mu_{q})))
\right]
\]
- Second derivative of the likelihood
\[
E\left( \frac{\partial^2\ell}{\partial a^2}\right) = -\sum_{q}
\left\{\left(\theta_{q}-b\right)\nu\mu_{q}\left(1-\mu_{q}\right)\right\}^{2}f_{q}\left[
\psi_{1}(\mu_{q}\nu)+\psi_{1}(\nu(1-\mu_{q})) \right] \\
E\left(\frac{\partial^2\ell}{\partial a \partial b}\right) = a\sum_{q}
\left(\theta_{q}-b\right)\left\{\nu\mu_{q}\left(1-\mu_{q}\right)\right\}^{2}f_{q}\left[
\psi_{1}(\mu_{q}\nu)+\psi_{1}(\nu(1-\mu_{q})) \right] \\
E\left(\frac{\partial^2\ell}{\partial a \partial \nu}\right) = -\sum_{q}
\left(\theta_{q}-b\right)\nu\mu_{q}\left(1-\mu_{q}\right)f_q\left[
\mu_{q}\psi_{1}(\mu_{q}\nu)-(1-\mu_{q})\psi_{1}(\nu(1-\mu_{q})) \right]
\\
E\left(\frac{\partial^2\ell}{\partial b^2}\right) = -a^{2}\sum_{q}
\left\{\nu\mu_{q}\left(1-\mu_{q}\right)\right\}^{2}f_{q}\left[
\psi_{1}(\mu_{q}\nu)+\psi_{1}(\nu(1-\mu_{q})) \right] \\
E\left(\frac{\partial^2\ell}{\partial b \partial \nu}\right) = a\sum_{q}
\nu\mu_{q}\left(1-\mu_{q}\right)f_q\left[
\mu_{q}\psi_{1}(\mu_{q}\nu)-(1-\mu_{q})\psi_{1}(\nu(1-\mu_{q})) \right]
\\
E\left(\frac{\partial^2\ell}{\partial \nu^2}\right) = N\psi_{1}(\nu) -
\sum_{q}f_q\left[
\mu_{q}^{2}\psi_{1}(\mu_{q}\nu)+(1-\mu_{q})^{2}\psi_{1}(\nu(1-\mu_{q}))
\right]
\]
If we use \(\xi\) instead of \(\nu\),
\[
E\left(\frac{\partial^2\ell}{\partial a \partial \xi}\right) = -\sum_{q}
\left(\theta_{q}-b\right)\nu^{2}\mu_{q}\left(1-\mu_{q}\right)f_q\left[
\mu_{q}\psi_{1}(\mu_{q}\nu)-(1-\mu_{q})\psi_{1}(\nu(1-\mu_{q})) \right]
\\
E\left(\frac{\partial^2\ell}{\partial b \partial \xi}\right) = a\sum_{q}
\nu^{2}\mu_{q}\left(1-\mu_{q}\right)f_q\left[
\mu_{q}\psi_{1}(\mu_{q}\nu)-(1-\mu_{q})\psi_{1}(\nu(1-\mu_{q})) \right]
\\
E\left(\frac{\partial^2\ell}{\partial \xi^2}\right) =
N\nu^{2}\psi_{1}(\nu) - \nu^{2}\sum_{q}f_q\left[
\mu_{q}^{2}\psi_{1}(\mu_{q}\nu)+(1-\mu_{q})^{2}\psi_{1}(\nu(1-\mu_{q}))
\right]
\]
The function DataGeneration can be used in a preparation
step. This function returns a set of artificial data and the true
parameters underlying the data.
Alldata <- DataGeneration(N=1000,
nitem_C = 8,
latent_dist = "2NM",
a_l = .3,
a_u = .7,
d = 1.664,
sd_ratio = 2,
prob = 0.3)
data <- Alldata$data_C
theta <- Alldata$theta
colnames(data) <- paste0("item", 1:8)
- Analysis (parameter estimation)
Mod1 <- IRTest_Cont(data = data,
latent_dist = "KDE")
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 1e-04 on 31st iterations.
#>
#> Model Fit:
#> log-likeli 2825.34
#> deviance -5650.679
#> AIC -5600.679
#> BIC -5477.985
#> HQ -5554.047
#>
#> The Number of Parameters:
#> item 24
#> dist 1
#> total 25
#>
#> The Number of Items: 8
#>
#> The Estimated Latent Distribution:
#> method - KDE
#> ----------------------------------------
#>
#>
#>
#> . . . . .
#> @ @ @ @ @ @ @ .
#> @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ .
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> +---------+---------+---------+---------+
#> -2 -1 0 1 2
### Log-likelihood
logLik(Mod1)
#> [1] 2825.34
### The estimated item parameters
coef(Mod1)
#> a b nu
#> item1 0.5585796 1.2924295 8.470660
#> item2 0.3406646 -0.6946243 5.423534
#> item3 0.4047930 -0.1765202 11.279811
#> item4 0.4839899 -0.1150749 9.674475
#> item5 0.3109872 -1.1347071 9.879757
#> item6 0.4666123 0.1573268 9.476279
#> item7 0.6554120 0.3054941 7.947572
#> item8 0.3897490 1.3419728 4.469948
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "WLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
### Standard errors of ability parameter estimates
plot(fscore$theta, fscore$theta_se)
- The result of latent distribution estimation
plot(Mod1, mapping = aes(colour="Estimated"), linewidth = 1) +
lims(y = c(0, .75))+
geom_line(
mapping=aes(
x=seq(-6,6,length=121),
y=dist2(seq(-6,6,length=121), prob = .3, d = 1.664, sd_ratio = 2),
colour="True"),
linewidth = 1)+
labs(title="The estimated latent density using '2NM'", colour= "Type")+
theme_bw()
p1 <- plot_item(Mod1,1)
p2 <- plot_item(Mod1,2)
p3 <- plot_item(Mod1,3)
p4 <- plot_item(Mod1,4)
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
reliability(Mod1)
#> $summed.score.scale
#> $summed.score.scale$test
#> NULL
#>
#> $summed.score.scale$item
#> NULL
#>
#>
#> $theta.scale
#> test reliability
#> 0.7923939
- Posterior distributions for the examinees
Each examinee’s posterior distribution is identified in the E-step of
the estimation algorithm (i.e., EM algorithm). Posterior distributions
can be found in Mod1$Pk.
set.seed(1)
selected_examinees <- sample(1:1000,6)
post_sample <-
data.frame(
X = rep(seq(-6,6, length.out=121),6),
prior = rep(Mod1$Ak/(Mod1$quad[2]-Mod1$quad[1]), 6),
posterior = 10*c(t(Mod1$Pk[selected_examinees,])),
ID = rep(paste("examinee", selected_examinees), each=121)
)
ggplot(data=post_sample, mapping=aes(x=X))+
geom_line(mapping=aes(y=posterior, group=ID, color='Posterior'))+
geom_line(mapping=aes(y=prior, group=ID, color='Prior'))+
labs(title="Posterior densities for selected examinees", x=expression(theta), y='density')+
facet_wrap(~ID, ncol=2)+
theme_bw()
- Test information function
ggplot()+
stat_function(
fun = inform_f_test,
args = list(Mod1)
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 1),
mapping = aes(color="Item 1")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 2),
mapping = aes(color="Item 2")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 3),
mapping = aes(color="Item 3")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 4),
mapping = aes(color="Item 4")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 5),
mapping = aes(color="Item 5")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 6),
mapping = aes(color="Item 6")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 7),
mapping = aes(color="Item 7")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 8),
mapping = aes(color="Item 8")
)+
lims(x=c(-6,6))+
labs(title="Test information function", x=expression(theta), y='information')+
theme_bw()
4. Mixed-format test
As in the cases of dichotomous and polytomous items, the function
DataGeneration can be used in the preparation step. This
function returns artificial data and some useful objects for analysis
(i.e., theta, data_D, item_D,
data_P, & item_P).
Alldata <- DataGeneration(model_D = 2,
model_P = "GRM",
N=1000,
nitem_D = 10,
nitem_P = 5,
latent_dist = "2NM",
d = 1.664,
sd_ratio = 1,
prob = 0.5)
DataD <- Alldata$data_D
DataP <- Alldata$data_P
theta <- Alldata$theta
colnames(DataD) <- paste0("item", 1:10)
colnames(DataP) <- paste0("item", 1:5)
- Analysis (parameter estimation)
Mod1 <- IRTest_Mix(data_D = DataD,
data_P = DataP,
model_D = "2PL",
model_P = "GRM",
latent_dist = "KDE")
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 1e-04 on 43rd iterations.
#>
#> Model Fit:
#> log-likeli -2832829
#> deviance 5665658
#> AIC 5665750
#> BIC 5665976
#> HQ 5665836
#>
#> The Number of Parameters:
#> item 45
#> dist 1
#> total 46
#>
#> The Number of Items:
#> dichotomous 10
#> polyotomous 5
#>
#> The Estimated Latent Distribution:
#> method - KDE
#> ----------------------------------------
#>
#>
#>
#>
#> . @ @ @ . . . . . .
#> . @ @ @ @ @ @ @ @ @ @ @ .
#> . @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> +---------+---------+---------+---------+
#> -2 -1 0 1 2
### Log-likelihood
logLik(Mod1)
#> [1] -2832829
### The estimated item parameters
coef(Mod1)
#> $Dichotomous
#> a b c
#> item1 0.7891202 1.2625485 0
#> item2 1.2499965 -0.6495540 0
#> item3 1.8413035 -0.2451646 0
#> item4 1.3725921 -0.2179645 0
#> item5 2.3290534 -1.2186082 0
#> item6 1.2979304 0.1663487 0
#> item7 1.1538293 0.3064596 0
#> item8 1.0723793 1.2980064 0
#> item9 2.3294905 0.1395925 0
#> item10 2.6566442 -0.9381758 0
#>
#> $Polytomous
#> a b_1 b_2 b_3 b_4
#> item1 1.7960156 -1.3772222 -0.4490053 -0.07115149 -0.02269190
#> item2 2.6132883 -1.2435968 -0.5217730 -0.07802599 0.64758495
#> item3 0.9616471 -0.9965864 -0.9735484 -0.39246415 0.01840465
#> item4 1.2488737 -1.1075645 -0.1732780 0.11281676 0.76505236
#> item5 1.8737454 -1.8756187 -1.3588250 -0.31849429 -0.25998056
### Standard errors of the item parameter estimates
coef_se(Mod1)
#> $Dichotomous
#> a b c
#> item1 0.07934206 0.13525483 NA
#> item2 0.09168355 0.06571107 NA
#> item3 0.11346457 0.04489424 NA
#> item4 0.09237972 0.05525207 NA
#> item5 0.18356284 0.05217865 NA
#> item6 0.08899798 0.05764454 NA
#> item7 0.08430875 0.06452792 NA
#> item8 0.09291723 0.10476652 NA
#> item9 0.13860432 0.03870397 NA
#> item10 0.18921268 0.04016531 NA
#>
#> $Polytomous
#> a b_1 b_2 b_3 b_4
#> item1 0.09208419 0.06469868 0.04524516 0.04436168 0.04457621
#> item2 0.10609490 0.04453496 0.03528593 0.03442853 0.03772448
#> item3 0.07092293 0.09430095 0.09328353 0.07468456 0.07197619
#> item4 0.07098634 0.07674239 0.05756408 0.05754074 0.06819791
#> item5 0.09945962 0.08399190 0.06140695 0.04355022 0.04347977
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "MLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
### Standard errors of ability parameter estimates
plot(fscore$theta, fscore$theta_se)
- The result of latent distribution estimation
plot(Mod1, mapping = aes(colour="Estimated"), linewidth = 1) +
stat_function(
fun = dist2,
args = list(prob = .5, d = 1.664, sd_ratio = 1),
mapping = aes(colour = "True"),
linewidth = 1) +
lims(y = c(0, .75)) +
labs(title="The estimated latent density using '2NM'", colour= "Type")+
theme_bw()
p1 <- plot_item(Mod1,1, type="d")
p2 <- plot_item(Mod1,4, type="d")
p3 <- plot_item(Mod1,8, type="d")
p4 <- plot_item(Mod1,10, type="d")
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
p1 <- plot_item(Mod1,1, type="p")
p2 <- plot_item(Mod1,2, type="p")
p3 <- plot_item(Mod1,3, type="p")
p4 <- plot_item(Mod1,4, type="p")
grid.arrange(p1, p2, p3, p4, ncol=2, nrow=2)
item_fit(Mod1)
#> $Dichotomous
#> stat df p.value
#> item1 7.789464 7 0.3515
#> item2 4.200670 7 0.7564
#> item3 9.025276 7 0.2508
#> item4 10.302423 7 0.1721
#> item5 8.308049 7 0.3062
#> item6 16.019905 7 0.0249
#> item7 29.980665 7 0.0001
#> item8 19.251959 7 0.0074
#> item9 13.524989 7 0.0603
#> item10 9.733639 7 0.2042
#>
#> $Polytomous
#> stat df p.value
#> item1 31.18450 31 0.4569
#> item2 49.30446 31 0.0196
#> item3 43.98313 31 0.0612
#> item4 43.01137 31 0.0741
#> item5 46.51346 31 0.0363
reliability(Mod1)
#> $summed.score.scale
#> $summed.score.scale$test
#> test reliability
#> 0.8521683
#>
#> $summed.score.scale$item
#> item1_D item2_D item3_D item4_D item5_D item6_D item7_D item8_D
#> 0.1105553 0.2355015 0.3930709 0.2830257 0.3620972 0.2659396 0.2252809 0.1663649
#> item9_D item10_D item1_P item2_P item3_P item4_P item5_P
#> 0.4929316 0.4546042 0.4444386 0.6534434 0.1920587 0.3118755 0.4410517
#>
#>
#> $theta.scale
#> test reliability
#> 0.8529496
- Posterior distributions for the examinees
Each examinee’s posterior distribution is identified in the E-step of
the estimation algorithm (i.e., EM algorithm). Posterior distributions
can be found in Mod1$Pk.
set.seed(1)
selected_examinees <- sample(1:1000,6)
post_sample <-
data.frame(
X = rep(seq(-6,6, length.out=121),6),
prior = rep(Mod1$Ak/(Mod1$quad[2]-Mod1$quad[1]), 6),
posterior = 10*c(t(Mod1$Pk[selected_examinees,])),
ID = rep(paste("examinee", selected_examinees), each=121)
)
ggplot(data=post_sample, mapping=aes(x=X))+
geom_line(mapping=aes(y=posterior, group=ID, color='Posterior'))+
geom_line(mapping=aes(y=prior, group=ID, color='Prior'))+
labs(title="Posterior densities for selected examinees", x=expression(theta), y='density')+
facet_wrap(~ID, ncol=2)+
theme_bw()
- Test information function
ggplot()+
stat_function(
fun = inform_f_test,
args = list(Mod1)
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 1, "d"),
mapping = aes(color="Dichotomous Item 1")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 2, "d"),
mapping = aes(color="Dichotomous Item 2")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 3, "d"),
mapping = aes(color="Dichotomous Item 3")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 1, "p"),
mapping = aes(color="Polytomous Item 1")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 2, "p"),
mapping = aes(color="Polytomous Item 2")
)+
stat_function(
fun=inform_f_item,
args = list(Mod1, 3, "p"),
mapping = aes(color="Polytomous Item 3")
)+
lims(x=c(-6,6))+
labs(title="Test information function", x=expression(theta), y='information')+
theme_bw()
