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)
#> Warning: package 'ggplot2' was built under R version 4.3.3
library(gridExtra)

0. Introduction

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.


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)


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.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)

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.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(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

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()

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)


Mod1 <- IRTest_Poly(data = data,
                    model = "GRM",
                    latent_dist = "KDE")


### 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)

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  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(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

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()

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

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 \]

\[ 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] \]

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.

\[ \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] \]

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


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.296 
#>    deviance   -5650.592 
#>         AIC   -5600.592 
#>         BIC   -5477.898 
#>          HQ   -5553.96 
#> 
#> 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.296

### The estimated item parameters
coef(Mod1)
#>               a          b        nu
#> item1 0.5585459  1.2924969  8.470383
#> item2 0.3406422 -0.6946817  5.423495
#> item3 0.4047894 -0.1765341 11.280134
#> item4 0.4839787 -0.1150889  9.674621
#> item5 0.3109717 -1.1347742  9.879742
#> item6 0.4666147  0.1573137  9.476702
#> item7 0.6553871  0.3054921  7.947426
#> item8 0.3897131  1.3420839  4.469861

### 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)

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.7923788

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()

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)




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

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.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(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

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()

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 <- 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))
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.369661e-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.221764e+00  0.0223
#> DC6  -5382.533 10765.07 10817.07 10944.67 10865.56     26  1.159148e+01  0.0007
#> DC7  -5383.124 10766.25 10820.25 10952.76 10870.61     27 -1.181958e+00  1.0000
#> DC8  -5395.022 10790.04 10846.04 10983.46 10898.27     28 -2.379597e+01  1.0000
#> DC9  -5380.459 10760.92 10818.92 10961.24 10873.01     29  2.912513e+01  0.0000
#> DC10 -5387.178 10774.36 10834.36 10981.59 10890.31     30 -1.343746e+01  1.0000
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.56
#> DC7  10870.61
#> DC8  10898.27
#> DC9  10873.01
#> DC10 10890.31
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