# 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.0.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_Dich is the estimation function when all items are dichotomously scored.

• IRTest_Poly is the estimation function when all items are polytomously scored.

• IRTest_Cont is the estimation function when all items are continuously scored.

• IRTest_Mix is the estimation function for a mixed-format test, a test comprising both dichotomous item(s) and polytomous item(s).

• factor_score estimates factor scores of examinees.

• coef_se returns standard errors of item parameter estimates.

• best_model selects the best model using an evaluation criterion.

• item_fit tests the statistical fit of all items individually.

• inform_f_item calculates the information value(s) of an item.

• inform_f_test calculates the information value(s) of a test.

• plot_item draws item response function(s) of an item.

• reliability calculates marginal reliability coefficient of IRT.

• latent_distribution returns evaluated PDF value(s) of an estimated latent distribution.

• DataGeneration generates several objects that can be useful for computer simulation studies. Among these are simulated item parameters, ability parameters and the corresponding item-response data.

• dist2 is a probability density function of two-component Gaussian mixture distribution.

• original_par_2GM converts re-parameterized parameters of two-component Gaussian mixture distribution into original parameters.

• cat_clps recommends category collapsing based on item parameters (or, equivalently, item response functions).

• recategorize implements the category collapsing.

• For S3 methods, anova, coef, logLik, plot, print, and summary are 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)
Mod1 <- IRTest_Dich(data = data,
model = 2,
latent_dist = "LLS",
h=4)

• 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) ### 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() • 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)
Mod1 <- IRTest_Poly(data = data,
model = "GRM",
latent_dist = "KDE")

• 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.6927174  0.9969194  1.06205110         NA          NA         NA
#> item2  1.7800947 -1.2013643  1.82111903         NA          NA         NA
#> item3  0.7669264 -3.5001678 -1.61657371         NA          NA         NA
#> item4  1.0820871 -1.6034036 -0.32237336         NA          NA         NA
#> item5  2.4860088  1.1361240  2.59875555         NA          NA         NA
#> item6  2.2625003 -0.5346921 -0.32758612         NA          NA         NA
#> item7  2.2570703  1.0684211  1.21711094         NA          NA         NA
#> item8  2.2387209 -1.2063932 -0.76041872 -0.4095813  0.51061883 0.61617892
#> item9  2.0585433  0.9851166  1.18638827  1.2777290  2.46807985 2.94760698
#> item10 1.9954821 -1.6062525 -1.13918150 -0.6967430 -0.60024205 0.60670105
#> item11 1.2312826 -2.7664081 -1.57784994 -0.9757298 -0.02755462 0.02436854
#> item12 1.0241866 -1.8145487 -1.21604456 -0.1067477  0.22886929 0.39358421
#> item13 0.9936170 -1.0079401  0.02655751  0.4302467  0.82603078 1.43435312
#> item14 2.3405434 -0.7925176  0.09595055  0.1125307  0.36411266 0.69098174
#>               b_6
#> item1          NA
#> item2          NA
#> item3          NA
#> item4          NA
#> item5          NA
#> item6          NA
#> item7          NA
#> item8  0.90461124
#> item9          NA
#> item10 0.83694606
#> item11 0.02869944
#> item12 2.11605587
#> item13 2.29126888
#> item14 1.05494694

### Standard errors of the item parameter estimates
coef_se(Mod1)
#>                 a        b_1        b_2        b_3        b_4        b_5
#> item1  0.12114675 0.06010133 0.06256077         NA         NA         NA
#> item2  0.10631563 0.06036759 0.08507624         NA         NA         NA
#> item3  0.08252368 0.35675759 0.16752036         NA         NA         NA
#> item4  0.07641223 0.11281918 0.06629506         NA         NA         NA
#> item5  0.17634306 0.04730486 0.13279360         NA         NA         NA
#> item6  0.12921027 0.04052895 0.03888963         NA         NA         NA
#> item7  0.16026198 0.04918148 0.05451057         NA         NA         NA
#> item8  0.09300996 0.04915604 0.04128802 0.03845833 0.03946031 0.04049672
#> item9  0.13200530 0.04875219 0.05452913 0.05771128 0.12313423 0.17342420
#> item10 0.08786155 0.06542056 0.05126873 0.04391077 0.04293820 0.04416089
#> item11 0.07477824 0.16876607 0.09624610 0.07246504 0.05842578 0.05868945
#> item12 0.06406466 0.12223803 0.09381491 0.06729118 0.06842015 0.07071505
#> item13 0.06399375 0.09006960 0.06887473 0.07265587 0.08262169 0.10680169
#> item14 0.09756599 0.04153892 0.03643059 0.03642742 0.03702934 0.03977123
#>               b_6
#> item1          NA
#> item2          NA
#> item3          NA
#> item4          NA
#> item5          NA
#> item6          NA
#> item7          NA
#> item8  0.04487770
#> item9          NA
#> item10 0.04784901
#> item11 0.05871696
#> item12 0.13902993
#> item13 0.15381263
#> item14 0.04599722

### 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() • 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.29337 15 0.1277 #> item2 21.13029 15 0.1327 #> item3 16.15910 15 0.3716 #> item4 15.77327 15 0.3973 #> item5 16.28395 15 0.3634 #> item6 19.00944 15 0.2133 #> item7 11.45525 15 0.7197 #> item8 52.70511 47 0.2629 #> item9 45.09418 39 0.2322 #> item10 57.70921 47 0.1361 #> item11 36.09294 47 0.8761 #> item12 63.70416 47 0.0526 #> item13 44.67282 47 0.5695 #> item14 54.55534 47 0.2093 • Reliability reliability(Mod1) #>$summed.score.scale
#> $summed.score.scale$test
#> test reliability
#>        0.8677075
#>
#> $summed.score.scale$item
#>      item1      item2      item3      item4      item5      item6      item7
#> 0.31133676 0.34645125 0.09162922 0.21941856 0.42950210 0.48845041 0.40308028
#>      item8      item9     item10     item11     item12     item13     item14
#> 0.58554508 0.38834451 0.51985244 0.28783877 0.24223806 0.22846845 0.59089272
#>
#>
#> $theta.scale #> test reliability #> 0.8895091 • 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}

\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)$

• 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]$

• Information

$E\left[\frac{\partial \log{f}}{\partial \theta} \right]^2 = (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]$

• 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=500,
nitem_C = 4,
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:4)

• Analysis (parameter estimation)
Mod1 <- IRTest_Cont(data = data,
latent_dist = "2NM",
range = c(-4,4),
q = 81,
threshold = 0.01)

• Results
### Summary
summary(Mod1)
#> Convergence:
#> Successfully converged below the threshold of 0.01 on 8th iterations.
#>
#> Model Fit:
#>  log-likeli   795.9182
#>    deviance   -1591.836
#>         AIC   -1561.836
#>         BIC   -1498.617
#>          HQ   -1537.029
#>
#> The Number of Parameters:
#>        item   12
#>        dist   3
#>       total   15
#>
#> The Number of Items:  4
#>
#> The Estimated Latent Distribution:
#> method - 2NM
#> ----------------------------------------
#>
#>
#>           .               . @
#>         @ @             . @ @ @
#>         @ @ @           @ @ @ @
#>       @ @ @ @         @ @ @ @ @ @
#>       @ @ @ @ @     . @ @ @ @ @ @
#>     . @ @ @ @ @ . . @ @ @ @ @ @ @ @
#>     @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#>   . @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ @ .
#> +---------+---------+---------+---------+
#> -2        -1        0         1         2

### Log-likelihood
logLik(Mod1)
#> [1] 795.9182

### The estimated item parameters
coef(Mod1)
#>               a          b        nu
#> item1 0.7477393  1.4019014  8.928225
#> item2 0.8252087 -0.7099602  8.695513
#> item3 1.1819349 -0.1936382 10.999908
#> item4 1.0115812 -0.1104953  9.667993

### 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() • 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.8900839 • 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:500,6)
post_sample <-
data.frame(
X = rep(seq(-4,4, length.out=81),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=81)
)

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")
)+
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 -2780534 #> deviance 5561068 #> AIC 5561160 #> BIC 5561386 #> HQ 5561246 #> #> 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] -2780534 ### The estimated item parameters coef(Mod1) #>$Dichotomous
#>                a          b c
#> item1  0.7683505  1.2938369 0
#> item2  1.2294115 -0.6496226 0
#> item3  1.7672321 -0.2392062 0
#> item4  1.3960361 -0.2078735 0
#> item5  2.2445332 -1.2317321 0
#> item6  1.2984903  0.1762492 0
#> item7  1.1526331  0.3156826 0
#> item8  1.0885619  1.2913119 0
#> item9  2.2773673  0.1571321 0
#> item10 2.5495257 -0.9462915 0
#>
#> $Polytomous #> a b_1 b_2 b_3 b_4 #> item1 1.8271920 -1.783584 0.1983415 0.9921879 1.0444947 #> item2 2.5845230 -2.530500 -1.0440054 -0.2307682 1.1980661 #> item3 0.8657413 -1.627685 -1.5700642 -0.4079450 0.4322696 #> item4 1.3558303 -1.933352 -0.2444436 0.2909519 1.8016323 #> item5 1.8299428 -2.515494 -1.6682141 0.4302294 0.5689776 ### Standard errors of the item parameter estimates coef_se(Mod1) #>$Dichotomous
#>                 a          b  c
#> item1  0.07937222 0.14077527 NA
#> item2  0.08999901 0.06686910 NA
#> item3  0.10830579 0.04645674 NA
#> item4  0.09220318 0.05474947 NA
#> item5  0.17969265 0.05450845 NA
#> item6  0.08828036 0.05763504 NA
#> item7  0.08384538 0.06445392 NA
#> item8  0.09462655 0.10262131 NA
#> item9  0.13326629 0.03955644 NA
#> item10 0.18436089 0.04145419 NA
#>
#> $Polytomous #> a b_1 b_2 b_3 b_4 #> item1 0.08767064 0.07747637 0.04418195 0.05281588 0.05406488 #> item2 0.10897017 0.11465869 0.04087772 0.03586859 0.04359767 #> item3 0.06510997 0.13486488 0.13143809 0.08208896 0.08430948 #> item4 0.07128387 0.10271517 0.05481474 0.05504227 0.09618334 #> item5 0.09350432 0.12346125 0.07273241 0.04574267 0.04719659 ### 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()

• 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.541279 7 0.3748 #> item2 6.020176 7 0.5374 #> item3 14.439173 7 0.0439 #> item4 8.471500 7 0.2929 #> item5 14.427626 7 0.0441 #> item6 8.708037 7 0.2743 #> item7 20.169505 7 0.0052 #> item8 14.487466 7 0.0432 #> item9 11.620862 7 0.1137 #> item10 8.916796 7 0.2587 #> #>$Polytomous
#>           stat df p.value
#> item1 35.50107 31  0.2644
#> item2 48.39362 31  0.0241
#> item3 24.25769 31  0.7999
#> item4 39.90463 31  0.1312
#> item5 35.88964 31  0.2499
• Reliability
reliability(Mod1)
#> $summed.score.scale #>$summed.score.scale$test #> test reliability #> 0.853241 #> #>$summed.score.scale$item #> item1_D item2_D item3_D item4_D item5_D item6_D item7_D item8_D #> 0.1041039 0.2343127 0.3845401 0.2940323 0.3499981 0.2687461 0.2261250 0.1659947 #> item9_D item10_D item1_P item2_P item3_P item4_P item5_P #> 0.4932106 0.4397854 0.4471110 0.6238494 0.1737165 0.3539179 0.4372312 #> #> #>$theta.scale
#> test reliability
#>        0.8754165
• 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_D
model_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.369668e-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.221771e+00  0.0223
#> DC6  -5382.711 10765.42 10817.42 10945.02 10865.92     26  1.123518e+01  0.0008
#> DC7  -5382.815 10765.63 10819.63 10952.14 10869.99     27 -2.071151e-01  1.0000
#> DC8  -5386.439 10772.88 10828.88 10966.30 10881.11     28 -7.248650e+00  1.0000
#> DC9  -5380.575 10761.15 10819.15 10961.47 10873.24     29  1.172887e+01  0.0006
#> DC10 -5378.354 10756.71 10816.71 10963.94 10872.67     30  4.440263e+00  0.0351
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  10865.92
#> DC7  10869.99
#> DC8  10881.11
#> DC9  10873.24
#> DC10 10872.67
• The best model overall
do.call(what = "best_model", args = c(model_fits[c(1:4,5)], criterion ="AIC"))
#> The best model: 2NM
#>
#>             AIC
#> Normal 10821.88
#> EHM    11030.50
#> 2NM    10801.16
#> KDM    10807.41
#> DC1    10823.88