IRT without the normality assumption
library(IRTest)
#> Thank you for using IRTest!
#> Please cite the package as:
#> Li, S. (2024). IRTest: Parameter estimation of item response theory with estimation of latent distribution (Version 2.1.0). R package.
#> URL: https://CRAN.R-project.org/package=IRTest
library(ggplot2)
library(gridExtra)0. Introduction
IRTest is a useful tool for $\mathcal{\color{red}{IRT}}$ (item response theory) parameter $\mathcal{\color{red}{est}}$imation, especially when the violation of normality assumption on latent distribution is suspected.
IRTest deals with uni-dimensional latent variable.
For missing values, IRTest adopts full information maximum likelihood (FIML) approach.
In IRTest, including the conventional usage of Gaussian distribution, several methods are available for estimation of latent distribution:
- empirical histogram method,
- two-component Gaussian mixture distribution,
- Davidian curve,
- kernel density estimation,
- log-linear smoothing.
Installation
The CRAN version of IRTest can be installed on R-console with:
install.packages("IRTest")For the development version, it can be installed on R-console with:
devtools::install_github("SeewooLi/IRTest")Functions
Followings are the functions of IRTest.
IRTest_Dichis the estimation function when all items are dichotomously scored.IRTest_Polyis the estimation function when all items are polytomously scored.IRTest_Contis the estimation function when all items are continuously scored.IRTest_Mixis the estimation function for a mixed-format test, a test comprising both dichotomous item(s) and polytomous item(s).factor_scoreestimates factor scores of examinees.coef_sereturns standard errors of item parameter estimates.best_modelselects the best model using an evaluation criterion.item_fittests the statistical fit of all items individually.inform_f_itemcalculates the information value(s) of an item.inform_f_testcalculates the information value(s) of a test.plot_itemdraws item response function(s) of an item.reliabilitycalculates marginal reliability coefficient of IRT.latent_distributionreturns evaluated PDF value(s) of an estimated latent distribution.DataGenerationgenerates several objects that can be useful for computer simulation studies. Among these are simulated item parameters, ability parameters and the corresponding item-response data.dist2is a probability density function of two-component Gaussian mixture distribution.original_par_2GMconverts re-parameterized parameters of two-component Gaussian mixture distribution into original parameters.cat_clpsrecommends category collapsing based on item parameters (or, equivalently, item response functions).recategorizeimplements the category collapsing.adaptive_testestimates ability parameter(s) which can be utilized for computer adaptive testing (CAT).For S3 methods,
anova,coef,logLik,plot,print, andsummaryare available.
1. Dichotomous items
- Preparing data
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)
- Results
### 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.329453946 0
#> item2 2.2856085 -0.687404395 0
#> item3 1.1690163 -0.215261808 0
#> item4 0.8122027 0.003225478 0
#> item5 1.6372745 -1.189646902 0
#> item6 1.2152174 0.121197279 0
#> item7 1.5656468 0.360962860 0
#> item8 2.5239591 1.182579616 0
#> item9 2.3468151 0.148729212 0
#> item10 1.0642602 -0.894474997 0
#> item11 2.2604206 1.540380888 0
#> item12 1.6180702 -0.263752931 0
#> item13 1.5673422 0.147437154 0
#> item14 1.8951603 -1.107805359 0
#> item15 1.5037220 -0.179279341 0
### Standard errors of the item parameter estimates
coef_se(Mod1)
#> a b c
#> item1 0.09101569 0.11521202 NA
#> item2 0.14467865 0.04196281 NA
#> item3 0.08298966 0.06304486 NA
#> item4 0.07293208 0.08380740 NA
#> item5 0.12363246 0.06673970 NA
#> item6 0.08463917 0.06036465 NA
#> item7 0.10018748 0.05100714 NA
#> item8 0.20149930 0.04677128 NA
#> item9 0.13758904 0.03887235 NA
#> item10 0.08506922 0.08394287 NA
#> item11 0.21591447 0.07116716 NA
#> item12 0.10017916 0.04990388 NA
#> item13 0.09811742 0.05014992 NA
#> item14 0.13792130 0.05618842 NA
#> item15 0.09501084 0.05207528 NA
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "MLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
- 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()
- Item response function
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
item_fit(Mod1)
#> stat df p.value
#> item1 11.96161 7 0.1018
#> item2 30.48845 7 0.0001
#> item3 12.97163 7 0.0728
#> item4 11.04379 7 0.1367
#> item5 16.82330 7 0.0186
#> item6 10.73185 7 0.1508
#> item7 15.91906 7 0.0259
#> item8 35.92519 7 0.0000
#> item9 19.94982 7 0.0057
#> item10 16.02559 7 0.0249
#> item11 24.35735 7 0.0010
#> item12 18.96929 7 0.0083
#> item13 18.77398 7 0.0089
#> item14 17.79824 7 0.0129
#> item15 13.57152 7 0.0593- Reliability
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.2791275 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
- Preparing data
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)
- Results
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 1e-04 on 48th iterations.
#>
#> Model Fit:
#> log-likeli -14017.91
#> deviance 28035.82
#> AIC 28175.82
#> BIC 28519.36
#> HQ 28306.39
#>
#> The Number of Parameters:
#> item 69
#> dist 1
#> total 70
#>
#> The Number of Items: 14
#>
#> The Estimated Latent Distribution:
#> method - KDE
#> ----------------------------------------
#>
#>
#>
#> .
#> . . . . . @ @ @ @ .
#> . @ @ @ @ @ @ @ @ @ @ @
#> . @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @
#> @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> +---------+---------+---------+---------+
#> -2 -1 0 1 2
### Log-likelihood
logLik(Mod1)
#> [1] -14017.91
### The estimated item parameters
coef(Mod1)
#> a b_1 b_2 b_3 b_4 b_5
#> item1 1.6927195 0.9969061 1.06203769 NA NA NA
#> item2 1.7801070 -1.2013789 1.82109854 NA NA NA
#> item3 0.7669308 -3.5001643 -1.61658138 NA NA NA
#> item4 1.0820874 -1.6034180 -0.32238872 NA NA NA
#> item5 2.4860192 1.1361107 2.59872923 NA NA NA
#> item6 2.2624780 -0.5347136 -0.32760542 NA NA NA
#> item7 2.2570721 1.0684089 1.21709835 NA NA NA
#> item8 2.2387170 -1.2064121 -0.76043872 -0.4095999 0.51060630 0.61616670
#> item9 2.0585506 0.9851031 1.18637400 1.2777143 2.46805763 2.94758072
#> item10 1.9954823 -1.6062659 -1.13919859 -0.6967613 -0.60026029 0.60668805
#> item11 1.2312860 -2.7664141 -1.57786173 -0.9757435 -0.02756969 0.02435341
#> item12 1.0241874 -1.8145620 -1.21605900 -0.1067629 0.22885408 0.39356896
#> item13 0.9936194 -1.0079542 0.02654204 0.4302308 0.82601434 1.43433573
#> item14 2.3405366 -0.7925391 0.09593596 0.1125162 0.36409964 0.69096989
#> b_6
#> item1 NA
#> item2 NA
#> item3 NA
#> item4 NA
#> item5 NA
#> item6 NA
#> item7 NA
#> item8 0.90459947
#> item9 NA
#> item10 0.83693332
#> item11 0.02868431
#> item12 2.11603958
#> item13 2.29124955
#> item14 1.05493556
### Standard errors of the item parameter estimates
coef_se(Mod1)
#> a b_1 b_2 b_3 b_4 b_5
#> item1 0.12114689 0.06010110 0.06256052 NA NA NA
#> item2 0.10631664 0.06036684 0.08507554 NA NA NA
#> item3 0.08252412 0.35675536 0.16751913 NA NA NA
#> item4 0.07641219 0.11281891 0.06629505 NA NA NA
#> item5 0.17634471 0.04730457 0.13279301 NA NA NA
#> item6 0.12920803 0.04052923 0.03889003 NA NA NA
#> item7 0.16026266 0.04918132 0.05451039 NA NA NA
#> item8 0.09300971 0.04915573 0.04128798 0.03845849 0.03946038 0.04049676
#> item9 0.13200630 0.04875196 0.05452885 0.05771099 0.12313391 0.17342377
#> item10 0.08786163 0.06542019 0.05126847 0.04391073 0.04293822 0.04416089
#> item11 0.07477836 0.16876562 0.09624564 0.07246472 0.05842571 0.05868938
#> item12 0.06406471 0.12223779 0.09381473 0.06729116 0.06842014 0.07071503
#> item13 0.06399374 0.09006920 0.06887461 0.07265572 0.08262144 0.10680126
#> item14 0.09756547 0.04153885 0.03643082 0.03642764 0.03702948 0.03977127
#> b_6
#> item1 NA
#> item2 NA
#> item3 NA
#> item4 NA
#> item5 NA
#> item6 NA
#> item7 NA
#> item8 0.04487767
#> item9 NA
#> item10 0.04784895
#> item11 0.05871689
#> item12 0.13902974
#> item13 0.15381193
#> item14 0.04599720
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "MLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
- 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()
- Item response function
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
item_fit(Mod1)
#> stat df p.value
#> item1 21.29236 15 0.1277
#> item2 21.13152 15 0.1327
#> item3 16.15904 15 0.3716
#> item4 15.77401 15 0.3972
#> item5 16.28425 15 0.3634
#> item6 19.01024 15 0.2133
#> item7 11.45574 15 0.7197
#> item8 52.70762 47 0.2628
#> item9 45.09359 39 0.2322
#> item10 57.71112 47 0.1361
#> item11 36.09235 47 0.8761
#> item12 63.70451 47 0.0526
#> item13 44.67274 47 0.5695
#> item14 54.55589 47 0.2092- Reliability
reliability(Mod1)
#> $summed.score.scale
#> $summed.score.scale$test
#> test reliability
#> 0.8677097
#>
#> $summed.score.scale$item
#> item1 item2 item3 item4 item5 item6 item7
#> 0.31133980 0.34645148 0.09162912 0.21941976 0.42950113 0.48845602 0.40308075
#> item8 item9 item10 item11 item12 item13 item14
#> 0.58555025 0.38834543 0.51985586 0.28784254 0.24224088 0.22847219 0.59089951
#>
#>
#> $theta.scale
#> test reliability
#> 0.88951- 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 ν = α + β,
$$ 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 μ and ν are parameters of the f(x), not variables.
Useful equations (click)
ψ(•) and ψ1(•) 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)] = ψ(β) − ψ(α + β).
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)
- Expected item response
$$ \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} \\ $$
- Gradient
log f = log Γ(ν) − log Γ(μν) − log Γ(ν(1 − μ)) + (μν − 1)log x + (ν(1 − μ) − 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] $$
- Information
$$ 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:
ℓ = ∑j∑qγjqlog Ljq,
where γjq = E[Pr (θj ∈ θq*)] is the expected probability of respondent j’s ability (θj) belonging to the θq* of the quadrature scheme and is calculated at the E-step of the MML-EM procedure, and Ljq is the likelihood of respondent j’s response at θ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 S1q = ∑jγjqlog xj and S2q = ∑jγjqlog (1 − xj). Since Eq[S1q] = fq[ψ(μqν)) − ψ(ν)] and Eq[S2q] = fq[ψ(ν(1 − μq))) − ψ(ν)], the expected values of the first derivatives are 0.
To keep ν positive, let ν = exp ξ; $\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 ξ instead of ν,
$$ 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] $$
- Preparing data
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)
- Results
### 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)
- 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()
- Item response function
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
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
- Preparing data
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")
- Results
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 1e-04 on 36th iterations.
#>
#> Model Fit:
#> log-likeli -2780566
#> deviance 5561131
#> AIC 5561223
#> BIC 5561449
#> HQ 5561309
#>
#> 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] -2780566
### The estimated item parameters
coef(Mod1)
#> $Dichotomous
#> a b c
#> item1 0.7683683 1.2939304 0
#> item2 1.2294183 -0.6494932 0
#> item3 1.7672442 -0.2390758 0
#> item4 1.3960448 -0.2077453 0
#> item5 2.2445345 -1.2316052 0
#> item6 1.2985066 0.1763732 0
#> item7 1.1526450 0.3158040 0
#> item8 1.0885917 1.2914036 0
#> item9 2.2773847 0.1572647 0
#> item10 2.5495176 -0.9461665 0
#>
#> $Polytomous
#> a b_1 b_2 b_3 b_4
#> item1 1.8272272 -1.783438 0.1984669 0.9922946 1.0446000
#> item2 2.5845616 -2.530349 -1.0438747 -0.2306363 1.1981656
#> item3 0.8657566 -1.627537 -1.5699170 -0.4078166 0.4323842
#> item4 1.3558507 -1.933205 -0.2443156 0.2910723 1.8017200
#> item5 1.8299748 -2.515338 -1.6680715 0.4303525 0.5690980
### Standard errors of the item parameter estimates
coef_se(Mod1)
#> $Dichotomous
#> a b c
#> item1 0.07937424 0.14077175 NA
#> item2 0.08999887 0.06686880 NA
#> item3 0.10830526 0.04645660 NA
#> item4 0.09220313 0.05474921 NA
#> item5 0.17969290 0.05450810 NA
#> item6 0.08828121 0.05763423 NA
#> item7 0.08384635 0.06445304 NA
#> item8 0.09462973 0.10261824 NA
#> item9 0.13326661 0.03955610 NA
#> item10 0.18435943 0.04145412 NA
#>
#> $Polytomous
#> a b_1 b_2 b_3 b_4
#> item1 0.08767258 0.07747476 0.04418122 0.05281440 0.05406337
#> item2 0.10897155 0.11465706 0.04087710 0.03586851 0.04359641
#> item3 0.06511079 0.13486219 0.13143549 0.08208761 0.08430779
#> item4 0.07128506 0.10271357 0.05481408 0.05504136 0.09618179
#> item5 0.09350612 0.12345899 0.07273097 0.04574168 0.04719542
### The estimated ability parameters
fscore <- factor_score(Mod1, ability_method = "MLE")
plot(theta, fscore$theta)
abline(b=1, a=0)
- 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()
- Item response function
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
item_fit(Mod1)
#> $Dichotomous
#> stat df p.value
#> item1 7.541907 7 0.3747
#> item2 6.020983 7 0.5373
#> item3 14.437397 7 0.0439
#> item4 8.472302 7 0.2928
#> item5 14.428052 7 0.0441
#> item6 8.709380 7 0.2742
#> item7 20.171976 7 0.0052
#> item8 14.488597 7 0.0431
#> item9 11.621903 7 0.1137
#> item10 8.916880 7 0.2587
#>
#> $Polytomous
#> stat df p.value
#> item1 35.50281 31 0.2643
#> item2 48.39715 31 0.0241
#> item3 24.25691 31 0.7999
#> item4 39.90793 31 0.1311
#> item5 35.89161 31 0.2498- Reliability
reliability(Mod1)
#> $summed.score.scale
#> $summed.score.scale$test
#> test reliability
#> 0.8532466
#>
#> $summed.score.scale$item
#> item1_D item2_D item3_D item4_D item5_D item6_D item7_D item8_D
#> 0.1041052 0.2343230 0.3845557 0.2940426 0.3500123 0.2687530 0.2261287 0.1659941
#> item9_D item10_D item1_P item2_P item3_P item4_P item5_P
#> 0.4932213 0.4398023 0.4471148 0.6238627 0.1737248 0.3539266 0.4372403
#>
#>
#> $theta.scale
#> test reliability
#> 0.8754227- 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()
5. Model comparison
- Data generation and model fitting
data <- DataGeneration(N=1000,
nitem_D = 10,
latent_dist = "2NM",
d = 1.664,
sd_ratio = 2,
prob = 0.3)$data_Dmodel_fits <- list()
model_fits[[1]] <- IRTest_Dich(data)
model_fits[[2]] <- IRTest_Dich(data, latent_dist = "EHM")
model_fits[[3]] <- IRTest_Dich(data, latent_dist = "2NM")
model_fits[[4]] <- IRTest_Dich(data, latent_dist = "KDE")
for(i in 1:10){
model_fits[[i+4]] <- IRTest_Dich(data, latent_dist = "DC", h = i)
}
names(model_fits) <- c("Normal", "EHM", "2NM", "KDM", paste0("DC", 1:10))- The best Davidian-curve model
do.call(what = "anova", args = model_fits[5:14])
#> Result of model comparison
#>
#> logLik deviance AIC BIC HQ n_pars chi p_value
#> DC1 -5390.940 10781.88 10823.88 10926.94 10863.05 21 NA NA
#> DC2 -5390.940 10781.88 10825.88 10933.85 10866.92 22 -9.369596e-05 1.0000
#> DC3 -5390.843 10781.69 10827.69 10940.56 10870.59 23 1.931828e-01 0.6603
#> DC4 -5390.940 10781.88 10829.88 10947.67 10874.65 24 -1.930907e-01 1.0000
#> DC5 -5388.329 10776.66 10826.66 10949.35 10873.29 25 5.221786e+00 0.0223
#> DC6 -5386.540 10773.08 10825.08 10952.68 10873.58 26 3.577166e+00 0.0586
#> DC7 -5384.421 10768.84 10822.84 10955.35 10873.21 27 4.237581e+00 0.0395
#> DC8 -5451.342 10902.68 10958.68 11096.10 11010.91 28 -1.338423e+02 1.0000
#> DC9 -5380.744 10761.49 10819.49 10961.81 10873.58 29 1.411962e+02 0.0000
#> DC10 -5377.908 10755.82 10815.82 10963.05 10871.78 30 5.672238e+00 0.0172
do.call(what = "best_model", args = model_fits[5:14])
#> The best model: DC1
#>
#> HQ
#> DC1 10863.05
#> DC2 10866.92
#> DC3 10870.59
#> DC4 10874.65
#> DC5 10873.29
#> DC6 10873.58
#> DC7 10873.21
#> DC8 11010.91
#> DC9 10873.58
#> DC10 10871.78- The best model overall