Note. The title comes from a click-bait headline generator, and I apologize for that. It was just an experiment.
Recently, I had to compute the population level mean and covariances for a set of items given the item parameters when the underlying response process is based on a multidimensional item response model. One can always approximate these values by simulating a dataset with a considerable sample size (e.g., 100,000) using the response model and then calculate the descriptive statistics. Given the large sample size, these numbers will be very close to population values. Instead, I was trying to calculate them directly from the model by integrating out the latent variables. The nice thing about this approach is that you can compute population-level parameters (e.g., SEM fit indices) if you are doing a simulation, or you can examine the plausibility of your soon-to-be generated datasets. Below, I will present the approach for compensatory and partially compensatory multidimensional IRT models.
Compensatory Models
Suppose you are working with a three-dimensional 4-PL compensatory item response model. In this model, the probability of a correct response can be written as,
\[ \begin{equation} P(Y=1 | \theta_1,\theta_2,\theta_3,a_1,a_2,a_3,d,g,u) = g + (u-g)\frac{1}{1 + e^{-(a_1\theta_1 + a_2 \theta_2+a_3\theta_3+d)}}\tag{1}, \end{equation} \]
in which \(\theta_1\), \(\theta_2\), and \(\theta_3\) are the coordinates for an individual on three latent dimensions; \(a_1\), \(a_2\), and \(a_3\) are the item discrimination parameters, \(g\) is the item guessing parameter (lower bound), \(u\) is the slipping parameter (upper bound), and \(d\) is the intercept parameter.
Below is a slightly modified table of item parameters from Reckase(2009) to generate data (Table 6.1).
ipar <- read.csv('data/ipar.csv')
ipar
a1 a2 a3 d g u
1 1.01 0.01 0.05 0.01 0.09 0.99
2 1.35 0.01 0.54 -0.58 0.07 0.98
3 1.38 0.09 0.47 -1.04 0.11 0.95
4 0.41 1.11 0.20 -0.37 0.08 0.97
5 0.29 1.24 0.02 0.49 0.18 0.96
6 0.05 0.48 0.00 0.29 0.12 0.99
7 0.00 0.19 1.19 -0.19 0.09 0.99
8 0.02 0.33 2.14 -1.84 0.01 1.00
9 0.00 0.22 0.90 0.51 0.05 1.00Download the item parameter file
Population-Level Means
We can compute the expected value (mean) for any single item by integrating the probability function over the joint distribution of \(\theta_1\), \(\theta_2\), and \(\theta_3\).
\[\begin{equation} E[Y]=\int_{\theta_1}\int_{\theta_2}\int_{\theta_3}P(Y=1)f(\theta_1,\theta_2,\theta_3)d_{\theta_1}d_{\theta_2}d_{\theta_3}\tag{2} \end{equation}\]
where \(f(\theta_1,\theta_2,\theta_3)\) is the joint density function for \(\theta_1\), \(\theta_2\), and \(\theta_3\) and \(P(Y=1)\) is given in Eq(1). In this post, it is assumed that \(\theta_1\), \(\theta_2\), and \(\theta_3\) follow a multivariate normal distribution with
I was recently experimenting with R packages for numerical integration, and I found the mvQuad package very easy to use and super fast. To do this integration using the mvQuad function, we first write a function that takes item parameters, person parameters, and μ and Σ for the multivariate normal distribution as arguments and returns the value of the following product. \[P(Y=1)f(\theta_1,\theta_2,\theta_3).\]
In this function:
thetais a 3-column matrix with columns representing \(\theta_1\), \(\theta_2\), \(\theta_3\), respectively.- \(a_1\), \(a_2\), and \(a_3\) are the discrimination parameters
- \(d\) is the intercept parameter
- \(g\) is the lower bound parameter
- \(u\) is the upper bound parameter
- M is the mean vector with a length of 3 for \(\theta_1\), \(\theta_2\), \(\theta_3\), respectively
- S is the 3 x 3 covariance matrix for \(\theta_1\), \(\theta_2\), \(\theta_3\)
M and S are internally passed to the dmvrnorm() function from the mvtnorm package to compute the joint density for \(\theta_1\), \(\theta_2\), \(\theta_3\).
Next, we generate a multivariate grid of Gauss-Hermite quadratures using the createNIGrid() function from the mvQuad package. This multivariate grid will be used for argument theta in the p.m4pl function. As we have three dimensions in our model and calculating a triple integration, the dim argument should be set to three. The type argument is set to GHe for the Gauss-Hermite rule. The level argument is the number of grid points you want to approximate the latent continuum. For instance, if we choose 3, three points on the continuum will be selected to represent the whole \(\theta\) scale. When the number of grid points increases, the precision of approximation to integration also increases at the expense of increased computational time. However, I found this package super fast, although I use 40 grid points, as you will see later. For now, I set it to three to demonstrate the structure. Also, even though it sounds few, three grid points do a fine job of approximation.
require(mvQuad)
nw <- createNIGrid(dim=3,
type="GHe",
level=3)
round(getNodes(nw),3)
[,1] [,2] [,3]
[1,] -1.732 -1.732 -1.732
[2,] 0.000 -1.732 -1.732
[3,] 1.732 -1.732 -1.732
[4,] -1.732 0.000 -1.732
[5,] 0.000 0.000 -1.732
[6,] 1.732 0.000 -1.732
[7,] -1.732 1.732 -1.732
[8,] 0.000 1.732 -1.732
[9,] 1.732 1.732 -1.732
[10,] -1.732 -1.732 0.000
[11,] 0.000 -1.732 0.000
[12,] 1.732 -1.732 0.000
[13,] -1.732 0.000 0.000
[14,] 0.000 0.000 0.000
[15,] 1.732 0.000 0.000
[16,] -1.732 1.732 0.000
[17,] 0.000 1.732 0.000
[18,] 1.732 1.732 0.000
[19,] -1.732 -1.732 1.732
[20,] 0.000 -1.732 1.732
[21,] 1.732 -1.732 1.732
[22,] -1.732 0.000 1.732
[23,] 0.000 0.000 1.732
[24,] 1.732 0.000 1.732
[25,] -1.732 1.732 1.732
[26,] 0.000 1.732 1.732
[27,] 1.732 1.732 1.732When we set the level argument to three, it picked three points representing each latent dimension: -1.732, 0, and 1.732. nw is a 27 x 3 matrix. Each column represents a dimension, and 27 (\(3^3\)) is the number of all possible combinations for these three grid points. If we set the level argument to 40, nw would be a \(40^3\) x 3 matrix.
You can also get the assigned weights for each combination of grid points.
round(getWeights(nw),3)
[,1]
[1,] 6.564
[2,] 5.858
[3,] 6.564
[4,] 5.858
[5,] 5.229
[6,] 5.858
[7,] 6.564
[8,] 5.858
[9,] 6.564
[10,] 5.858
[11,] 5.229
[12,] 5.858
[13,] 5.229
[14,] 4.667
[15,] 5.229
[16,] 5.858
[17,] 5.229
[18,] 5.858
[19,] 6.564
[20,] 5.858
[21,] 6.564
[22,] 5.858
[23,] 5.229
[24,] 5.858
[25,] 6.564
[26,] 5.858
[27,] 6.564After we create the function object (p.m4pl) to integrate over and an object for multivariate grid points (nw), we can now use the quadrature() function to compute the expected value. The first argument for the quadrature() function is the name of the function object. The second argument for the quadrature() function is the object for the multivariate grid which is passed to the function to be integrated as the first argument (theta). Then, all arguments expected by the function to be integrated are given to the quadrature() function. Below, I assume the latent dimensions do not correlate to each other, but we can simply modify the Sigma if we want to assume the latent dimensions are correlated.
mu <- c(0,0,0)
Sigma <- diag(3)
quadrature(p.m4pl,
grid = nw,
a1=1.01,
a2=0.01,
a3=0.05,
d=0.01,
g=0.09,
u=0.99,
M=mu,
S=Sigma)
[1] 0.5419This indicates that the expected value for an item with parameters of \((a_1,a_2,a_3,d,g,u) = (1.01, 0.01, 0.05, 0.01, 0.09, 0.99)\) would be 0.5419 assuming that three latent dimensions follow a multivariate normal distribution with zero means, unit variances, and zero correlations.
To make this easier, we wrap everything in a single function, and call this function e.m4pl.
e.m4pl <- function(a1,a2,a3,d,g,u,M,S,nquad){
p.m4pl <- function(theta,a1,a2,a3,d,g,u,M,S){
z = a1*theta[,1] + a2*theta[,2] + a3*theta[,3] + d
p = g + (u-g)*(1/(1+exp(-z)))
f = dmvnorm(theta,mean=M,sigma=S)
p*f
}
nw <- createNIGrid(dim=3, type="GHe", level=nquad)
quadrature(p.m4pl,
grid = nw,
a1=a1,
a2=a2,
a3=a3,
d=d,
g=g,
u=u,
M=M,
S=S)
}
Now, we can compute the expected value for any given set of item parameters when the responses are generated based on a multidimensional 4PL model as shown in Eq(1)
mu <- c(0,0,0)
Sigma <- diag(3)
i = 2 # item number
e.m4pl(a1 = ipar[i,]$a1,
a2 = ipar[i,]$a2,
a3 = ipar[i,]$a3,
d = ipar[i,]$d,
g = ipar[i,]$g,
u = ipar[i,]$u,
M = mu,
S = Sigma,
nquad = 40)
[1] 0.4315# Compute the population level means for all 9 items
exp.means <- c()
for(i in 1:9){
exp.means[i]= e.m4pl(a1 = ipar[i,]$a1,
a2 = ipar[i,]$a2,
a3 = ipar[i,]$a3,
d = ipar[i,]$d,
g = ipar[i,]$g,
u = ipar[i,]$u,
M = mu,
S = Sigma,
nquad = 40)
}
exp.means
[1] 0.5419 0.4315 0.3795 0.4614 0.6418 0.6144 0.5069 0.2610 0.6259Population-Level Variances and Covariances
We can compute the population level variance for any item using
\[\begin{equation} \sigma_{Y}^{2}=E[Y]-E[Y]^2.\tag{3} \end{equation}\]
This equation holds as \(Y\) is binary and can take only two values (0, 1). This is straightforward since we know how to compute \(E[Y]\). For instance, the population variance for Item 1 would be 0.2482.
mu <- c(0,0,0)
Sigma <- diag(3)
i = 1 # item number
e <- e.m4pl(a1 = ipar[i,]$a1,
a2 = ipar[i,]$a2,
a3 = ipar[i,]$a3,
d = ipar[i,]$d,
g = ipar[i,]$g,
u = ipar[i,]$u,
M = mu, S = Sigma, nquad = 40)
e - e^2
[1] 0.2482To compute the covariance for an item pair given the set of item parameters, we first need to write a new function for the joint probability to integrate over.
\[\begin{equation} E[Y_iY_j]=\int_{\theta_1}\int_{\theta_2}\int_{\theta_3}P(Y_i=1)P(Y_j=1)f(\theta_1,\theta_2,\theta_3)d_{\theta_1}d_{\theta_2}d_{\theta_3}\tag{4} \end{equation}\]
Below is what we get after tweaking a little bit of the above code written for \(E[Y]\).
e12.m4pl <- function(a11,a12,a13,d1,g1,u1,a21,a22,a23,d2,g2,u2,M,S,nquad){
p12.m4pl <- function(theta,a11,a12,a13,d1,g1,u1,a21,a22,a23,d2,g2,u2,M,S){
z1 = a11*theta[,1] + a12*theta[,2] + a13*theta[,3] + d1
z2 = a21*theta[,1] + a22*theta[,2] + a23*theta[,3] + d2
p1 = g1 + (u1-g1)*(1/(1+exp(-z1)))
p2 = g2 + (u2-g2)*(1/(1+exp(-z2)))
f = dmvnorm(theta,mean=M,sigma=S)
p1*p2*f
}
nw <- createNIGrid(dim=3, type="GHe", level=nquad)
quadrature(p12.m4pl,
grid = nw,
a11 = a11,
a12 = a12,
a13 = a13,
d1 = d1,
g1 = g1,
u1 = u1,
a21 = a21,
a22 = a22,
a23 = a23,
d2 = d2,
g2 = g2,
u2 = u2,
M = M,
S = S)
}
Now, we can compute the expected covariance for any item pair using
\[\begin{equation} \sigma_{Y_1Y_2}=E[Y_1Y_2]-E[Y_1]E[Y_2].\tag{5} \end{equation}\]
For instance, the population covariance between Item 1 and Item 2 would be 0.0414.
mu <- c(0,0,0)
Sigma <- diag(3)
i = 1 # item number for the first item
j = 2 # item number for the second item
e1 <- e.m4pl(a1 = ipar[i,]$a1,
a2 = ipar[i,]$a2,
a3 = ipar[i,]$a3,
d = ipar[i,]$d,
g = ipar[i,]$g,
u = ipar[i,]$u,
M = mu, S = Sigma, nquad = 40)
e2 <- e.m4pl(a1 = ipar[j,]$a1,
a2 = ipar[j,]$a2,
a3 = ipar[j,]$a3,
d = ipar[j,]$d,
g = ipar[j,]$g,
u = ipar[j,]$u,
M = mu, S = Sigma, nquad = 40)
e12 <- e12.m4pl(a11 = ipar[i,]$a1,
a12 = ipar[i,]$a2,
a13 = ipar[i,]$a3,
d1 = ipar[i,]$d,
g1 = ipar[i,]$g,
u1 = ipar[i,]$u,
a21 = ipar[j,]$a1,
a22 = ipar[j,]$a2,
a23 = ipar[j,]$a3,
d2 = ipar[j,]$d,
g2 = ipar[j,]$g,
u2 = ipar[j,]$u,
M = mu,
S = Sigma,
nquad = 40)
e12 - e1*e2
[1] 0.04137Let’s create a wrapper function using these smaller functions to return a population variance-covariance matrix for a given set of item parameters.
pop.vcov <- function(ip, M, S, nq){
VCOV <- matrix(nrow=nrow(ip),ncol=nrow(ip))
for(i in 1:nrow(ip)){
e <- e.m4pl(a1 = ip[i,]$a1,
a2 = ip[i,]$a2,
a3 = ip[i,]$a3,
d = ip[i,]$d,
g = ip[i,]$g,
u = ip[i,]$u,
M = M, S = S, nquad = nq)
VCOV[i,i] <- e - e^2
}
for(i in 1:(nrow(ip)-1)){
for(j in (i+1):nrow(ip)){
e1 <- e.m4pl(a1 = ip[i,]$a1,
a2 = ip[i,]$a2,
a3 = ip[i,]$a3,
d = ip[i,]$d,
g = ip[i,]$g,
u = ip[i,]$u,
M = M, S = S, nquad = nq)
e2 <- e.m4pl(a1 = ip[j,]$a1,
a2 = ip[j,]$a2,
a3 = ip[j,]$a3,
d = ip[j,]$d,
g = ip[j,]$g,
u = ip[j,]$u,
M = M, S = S, nquad = nq)
e12 <- e12.m4pl(a11 = ip[i,]$a1,
a12 = ip[i,]$a2,
a13 = ip[i,]$a3,
d1 = ip[i,]$d,
g1 = ip[i,]$g,
u1 = ip[i,]$u,
a21 = ip[j,]$a1,
a22 = ip[j,]$a2,
a23 = ip[j,]$a3,
d2 = ip[j,]$d,
g2 = ip[j,]$g,
u2 = ip[j,]$u,
M = M,
S = S,
nquad = nq)
VCOV[i,j] <- VCOV[j,i] <- e12 - e1*e2
}
}
return(VCOV=VCOV)
}
We can now compute the population variance-covariance matrix for all these 9 items, and then convert it to a correlation matrix using the cov2cor() function.
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.248 0.041 0.036 0.014 0.008 0.002 0.002 0.003 0.002
[2,] 0.041 0.245 0.047 0.018 0.009 0.002 0.018 0.022 0.015
[3,] 0.036 0.047 0.235 0.017 0.010 0.003 0.013 0.017 0.011
[4,] 0.014 0.018 0.017 0.249 0.037 0.019 0.013 0.016 0.014
[5,] 0.008 0.009 0.010 0.037 0.230 0.018 0.006 0.007 0.008
[6,] 0.002 0.002 0.003 0.019 0.018 0.237 0.003 0.004 0.004
[7,] 0.002 0.018 0.013 0.013 0.006 0.003 0.250 0.054 0.038
[8,] 0.003 0.022 0.017 0.016 0.007 0.004 0.054 0.193 0.042
[9,] 0.002 0.015 0.011 0.014 0.008 0.004 0.038 0.042 0.234 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 1.000 0.168 0.148 0.055 0.034 0.009 0.008 0.013 0.007
[2,] 0.168 1.000 0.197 0.073 0.039 0.010 0.072 0.100 0.062
[3,] 0.148 0.197 1.000 0.072 0.043 0.013 0.056 0.078 0.048
[4,] 0.055 0.073 0.072 1.000 0.155 0.079 0.053 0.072 0.057
[5,] 0.034 0.039 0.043 0.155 1.000 0.077 0.027 0.035 0.035
[6,] 0.009 0.010 0.013 0.079 0.077 1.000 0.013 0.017 0.017
[7,] 0.008 0.072 0.056 0.053 0.027 0.013 1.000 0.244 0.156
[8,] 0.013 0.100 0.078 0.072 0.035 0.017 0.244 1.000 0.199
[9,] 0.007 0.062 0.048 0.057 0.035 0.017 0.156 0.199 1.000[1] 1.6102 1.2309 1.1392 0.9374 0.8530 0.8449 0.8406 0.7982 0.7456plot(eigens)
Partially Compensatory Models
The process above can be applied to any probability model for item responses. The only thing that needs to be adjusted is how we define the probability of correct response given the item and person parameters. Suppose you are working with a two-dimensional 4-PL partially compensatory item response model. In this model, the probability of a correct response can be written as,
\[ \begin{equation} P(Y=1 | \theta_1,\theta_2,a_1,a_2,b_1,b_2,g,u) = g + (u-g) \left(\prod_{l=1}^{2}\frac{1}{1+e^{-a_l(\theta_l-b_l)}}\right)\tag{6}, \end{equation} \]
in which \(\theta_l\), \(a_l\), and \(b_l\) are the person parameter, item discrimination parameter, and item difficulty parameters on the \(l^{th}\) latent dimension, respectively; \(g\) is the item guessing parameter (lower bound), and \(u\) is the slipping parameter (upper bound). Below is a table of item parameters for some hypothetical items.
ipar <- read.csv('data/ipar2.csv')
ipar
a1 a2 b1 b2 g u
1 0.94 1.99 -1.36 -1.65 0.05 0.98
2 1.83 1.38 0.11 -2.63 0.04 0.97
3 1.75 1.25 -2.97 -1.09 0.10 0.95
4 1.10 0.95 -0.29 -0.08 0.09 0.96
5 1.89 1.92 -0.49 0.45 0.07 0.98Download the item parameter file
Below are the auxiliary and wrapper functions for this model to compute the population variance-covariance matrix and correlations.
# Function to compute the marginal expectation
e.pcm4pl <- function(a1,a2,b1,b2,g,u,M,S,nquad){
p.pcm4pl <- function(theta,a1,a2,b1,b2,d,g,u,M,S){
z1 = a1*(theta[,1]-b1)
z2 = a2*(theta[,2]-b2)
p = g + (u-g)*(1/(1+exp(-z1)))*(1/(1+exp(-z2)))
f = dmvnorm(theta,mean=M,sigma=S)
p*f
}
nw <- createNIGrid(dim=2, type="GHe", level=nquad)
quadrature(p.pcm4pl,
grid = nw,
a1=a1,
a2=a2,
b1=b1,
b2=b2,
g=g,
u=u,
M=M,
S=S)
}
# Function to compute the expectation for joint probability
e12.pcm4pl <- function(a11,a12,b11,b12,g1,u1,a21,a22,b21,b22,g2,u2,M,S,nquad){
p12.pcm4pl <- function(theta,a11,a12,b11,b12,g1,u1,a21,a22,b21,b22,g2,u2,M,S){
z11 = a11*(theta[,1]-b11)
z12 = a12*(theta[,2]-b12)
z21 = a21*(theta[,1]-b21)
z22 = a22*(theta[,2]-b22)
p1 = g1 + (u1-g1)*(1/(1+exp(-z11)))*(1/(1+exp(-z12)))
p2 = g2 + (u2-g2)*(1/(1+exp(-z21)))*(1/(1+exp(-z22)))
f = dmvnorm(theta,mean=M,sigma=S)
p1*p2*f
}
nw <- createNIGrid(dim=2, type="GHe", level=nquad)
quadrature(p12.pcm4pl,
grid = nw,
a11 = a11,
a12 = a12,
b11 = b11,
b12 = b12,
g1 = g1,
u1 = u1,
a21 = a21,
a22 = a22,
b21 = b21,
b22 = b22,
g2 = g2,
u2 = u2,
M = M,
S = S)
}
# Wrapper function
pop.vcov <- function(ip, M, S, nq){
VCOV <- matrix(nrow=nrow(ip),ncol=nrow(ip))
for(i in 1:nrow(ip)){
e <- e.pcm4pl(a1 = ip[i,]$a1,
a2 = ip[i,]$a2,
b1 = ip[i,]$b1,
b2 = ip[i,]$b2,
g = ip[i,]$g,
u = ip[i,]$u,
M = M, S = S, nquad = nq)
VCOV[i,i] <- e - e^2
}
for(i in 1:(nrow(ip)-1)){
for(j in (i+1):nrow(ip)){
e1 <- e.pcm4pl(a1 = ip[i,]$a1,
a2 = ip[i,]$a2,
b1 = ip[i,]$b1,
b2 = ip[i,]$b2,
g = ip[i,]$g,
u = ip[i,]$u,
M = M, S = S, nquad = nq)
e2 <- e.pcm4pl(a1 = ip[j,]$a1,
a2 = ip[j,]$a2,
b1 = ip[j,]$b1,
b2 = ip[j,]$b2,
g = ip[j,]$g,
u = ip[j,]$u,
M = M, S = S, nquad = nq)
e12 <- e12.pcm4pl(a11=ipar[i,1],
a12=ipar[i,2],
b11=ipar[i,3],
b12=ipar[i,4],
g1=ipar[i,5],
u1=ipar[i,6],
a21=ipar[j,1],
a22=ipar[j,2],
b21=ipar[j,3],
b22=ipar[j,4],
g2=ipar[j,5],
u2=ipar[j,6],
M=mu,S=Sigma,nquad=10)
VCOV[i,j] <- VCOV[j,i] <- e12 - e1*e2
}
}
return(VCOV=VCOV)
}
We can now compute the population level variance-covariance and correlations for these 5 items using the wrapper function.
[,1] [,2] [,3] [,4] [,5]
[1,] 0.221 0.037 0.021 0.023 0.027
[2,] 0.037 0.248 0.009 0.029 0.029
[3,] 0.021 0.009 0.203 0.018 0.027
[4,] 0.023 0.029 0.018 0.225 0.029
[5,] 0.027 0.029 0.027 0.029 0.204 [,1] [,2] [,3] [,4] [,5]
[1,] 1.000 0.158 0.100 0.101 0.125
[2,] 0.158 1.000 0.041 0.123 0.131
[3,] 0.100 0.041 1.000 0.083 0.132
[4,] 0.101 0.123 0.083 1.000 0.134
[5,] 0.125 0.131 0.132 0.134 1.000[1] 1.4562 0.9700 0.9054 0.8517 0.8167plot(eigens)
