Example data: Resampled with replacement data set of 635 older adults (age 80-100) self-reporting on 7 items assessing the Instrumental Activities of Daily Living (IADL) as follows:

  1. Housework (cleaning and laundry): 1=64%
  2. Bedmaking: 1=84%
  3. Cooking: 1=77%
  4. Everyday shopping: 1=66%
  5. Getting to places outside of walking distance: 1=65%
  6. Handling banking and other business: 1=73%
  7. Using the telephone 1=94%

Two versions of a response format were available:

Binary -> 0 = “needs help”, 1 = “does not need help”

Categorical -> 0 = “can’t do it”, 1=”big problems”, 2=”some problems”, 3=”no problems”

Higher scores indicate greater function. We will look at each response format in turn.

Package Installation and Loading

if (require(lavaan) == FALSE){
  install.packages("lavaan")
}
Loading required package: lavaan
This is lavaan 0.6-10
lavaan is FREE software! Please report any bugs.
library(lavaan)

if (require(mirt) == FALSE){
  install.packages("mirt")
}
Loading required package: mirt
Loading required package: stats4
Loading required package: lattice
library(mirt)

if (require(ggplot2) == FALSE){
  install.packages("ggplot2")
}
Loading required package: ggplot2
library(ggplot2)

Data Import into R

The data are in a text file named adl.dat orignally used in Mplus (so no column names were included at the top of the file). The file contains more items than we will use, so we select only items above from the whole file.


#read in data file (Mplus file format so having to label columns)
# adlData = read.table(file = "adlPermute.dat", header = FALSE, na.strings = ".", col.names = c("case", paste0("dpa", 1:14), paste0("dia", 1:7), paste0("cpa", 1:14), paste0("cia", 1:7)))

adlData = read.csv(file = "adlPermute.csv")

#select Situations items and PersonID variables
iadlDataInit = adlData[c(paste0("dia", 1:7))]

#remove cases with all items missing
removeCases = which(apply(X = iadlDataInit, MARGIN = 1, FUN = function (x){ length(which(is.na(x)))}) == 7)

iadlData = iadlDataInit[-removeCases,]

Estimation with Marginal Maximum Likelihood

We will introduce the mirt package as a method for estimating IRT models. Overall, the package is very good, but typically is used for scaling purposes (measurement rather than use of latent variables in additional model equations). We use the package to demonstrate estimating IRT models using marginal maximum likelihood. If you wish to use the latent trait estimates in secondary analyses (which you would otherwise use SEM for simultaneously), there are additional steps to take to ensure the error associated with each score is carried over to the subsequent analysese

When all the items of the model are the same type, the mirt syntax is very short. The mirt() function is used to provide estimates, with options model=1 for all items measuring the same trait and itemtype="2PL" for the two-parameter logistic model ("Rasch" is used for the 1PL shorthand). The "Rasch" designation estimates a model where the loadings are all set to one and the factor/latent trait variance is estimated) – which is an equivalent model to the one estimated below but we seek to keep the latent trait standardized. We will estimate both simultaneously here:

mirt1PLsyntax = "
IADL = 1-7
CONSTRAIN = (1-7, a1)
COV = 1
"

model1PLmirt = mirt(data = iadlData, model = mirt1PLsyntax)

Iteration: 1, Log-Lik: -1895.884, Max-Change: 1.91039
Iteration: 2, Log-Lik: -1564.411, Max-Change: 1.20995
Iteration: 3, Log-Lik: -1487.737, Max-Change: 0.59131
Iteration: 4, Log-Lik: -1463.809, Max-Change: 0.36754
Iteration: 5, Log-Lik: -1454.168, Max-Change: 0.25857
Iteration: 6, Log-Lik: -1449.751, Max-Change: 0.18492
Iteration: 7, Log-Lik: -1445.761, Max-Change: 0.06525
Iteration: 8, Log-Lik: -1445.320, Max-Change: 0.04955
Iteration: 9, Log-Lik: -1445.014, Max-Change: 0.03715
Iteration: 10, Log-Lik: -1444.484, Max-Change: 0.01260
Iteration: 11, Log-Lik: -1444.374, Max-Change: 0.01241
Iteration: 12, Log-Lik: -1444.282, Max-Change: 0.01211
Iteration: 13, Log-Lik: -1443.920, Max-Change: 0.00903
Iteration: 14, Log-Lik: -1443.900, Max-Change: 0.00794
Iteration: 15, Log-Lik: -1443.882, Max-Change: 0.00713
Iteration: 16, Log-Lik: -1443.811, Max-Change: 0.00350
Iteration: 17, Log-Lik: -1443.807, Max-Change: 0.00313
Iteration: 18, Log-Lik: -1443.803, Max-Change: 0.00288
Iteration: 19, Log-Lik: -1443.789, Max-Change: 0.00364
Iteration: 20, Log-Lik: -1443.788, Max-Change: 0.00155
Iteration: 21, Log-Lik: -1443.787, Max-Change: 0.00119
Iteration: 22, Log-Lik: -1443.786, Max-Change: 0.00247
Iteration: 23, Log-Lik: -1443.785, Max-Change: 0.00112
Iteration: 24, Log-Lik: -1443.785, Max-Change: 0.00086
Iteration: 25, Log-Lik: -1443.784, Max-Change: 0.00065
Iteration: 26, Log-Lik: -1443.784, Max-Change: 0.00066
Iteration: 27, Log-Lik: -1443.784, Max-Change: 0.00051
Iteration: 28, Log-Lik: -1443.783, Max-Change: 0.00084
Iteration: 29, Log-Lik: -1443.783, Max-Change: 0.00055
Iteration: 30, Log-Lik: -1443.783, Max-Change: 0.00046
Iteration: 31, Log-Lik: -1443.783, Max-Change: 0.00068
Iteration: 32, Log-Lik: -1443.783, Max-Change: 0.00051
Iteration: 33, Log-Lik: -1443.783, Max-Change: 0.00044
Iteration: 34, Log-Lik: -1443.783, Max-Change: 0.00037
Iteration: 35, Log-Lik: -1443.783, Max-Change: 0.00030
Iteration: 36, Log-Lik: -1443.783, Max-Change: 0.00025
Iteration: 37, Log-Lik: -1443.783, Max-Change: 0.00010
Iteration: 38, Log-Lik: -1443.783, Max-Change: 0.00009
model2PLmirt = mirt(data = iadlData, model = 1, itemtype = "2PL")

Iteration: 1, Log-Lik: -1895.884, Max-Change: 1.46672
Iteration: 2, Log-Lik: -1558.097, Max-Change: 1.27474
Iteration: 3, Log-Lik: -1477.872, Max-Change: 1.01881
Iteration: 4, Log-Lik: -1452.137, Max-Change: 0.81113
Iteration: 5, Log-Lik: -1441.831, Max-Change: 0.67794
Iteration: 6, Log-Lik: -1437.020, Max-Change: 0.54253
Iteration: 7, Log-Lik: -1434.593, Max-Change: 0.43074
Iteration: 8, Log-Lik: -1433.277, Max-Change: 0.34108
Iteration: 9, Log-Lik: -1432.514, Max-Change: 0.27723
Iteration: 10, Log-Lik: -1431.338, Max-Change: 0.10121
Iteration: 11, Log-Lik: -1431.215, Max-Change: 0.07949
Iteration: 12, Log-Lik: -1431.113, Max-Change: 0.06283
Iteration: 13, Log-Lik: -1430.856, Max-Change: 0.01857
Iteration: 14, Log-Lik: -1430.807, Max-Change: 0.01384
Iteration: 15, Log-Lik: -1430.765, Max-Change: 0.01339
Iteration: 16, Log-Lik: -1430.619, Max-Change: 0.01480
Iteration: 17, Log-Lik: -1430.603, Max-Change: 0.01345
Iteration: 18, Log-Lik: -1430.589, Max-Change: 0.01313
Iteration: 19, Log-Lik: -1430.528, Max-Change: 0.00775
Iteration: 20, Log-Lik: -1430.523, Max-Change: 0.00658
Iteration: 21, Log-Lik: -1430.519, Max-Change: 0.00721
Iteration: 22, Log-Lik: -1430.514, Max-Change: 0.00645
Iteration: 23, Log-Lik: -1430.512, Max-Change: 0.00610
Iteration: 24, Log-Lik: -1430.509, Max-Change: 0.00550
Iteration: 25, Log-Lik: -1430.503, Max-Change: 0.00461
Iteration: 26, Log-Lik: -1430.502, Max-Change: 0.00362
Iteration: 27, Log-Lik: -1430.501, Max-Change: 0.00248
Iteration: 28, Log-Lik: -1430.500, Max-Change: 0.00176
Iteration: 29, Log-Lik: -1430.499, Max-Change: 0.00264
Iteration: 30, Log-Lik: -1430.499, Max-Change: 0.00244
Iteration: 31, Log-Lik: -1430.496, Max-Change: 0.00168
Iteration: 32, Log-Lik: -1430.496, Max-Change: 0.00097
Iteration: 33, Log-Lik: -1430.496, Max-Change: 0.00121
Iteration: 34, Log-Lik: -1430.495, Max-Change: 0.00284
Iteration: 35, Log-Lik: -1430.495, Max-Change: 0.00068
Iteration: 36, Log-Lik: -1430.495, Max-Change: 0.00025
Iteration: 37, Log-Lik: -1430.495, Max-Change: 0.00021
Iteration: 38, Log-Lik: -1430.495, Max-Change: 0.00012
Iteration: 39, Log-Lik: -1430.495, Max-Change: 0.00059
Iteration: 40, Log-Lik: -1430.495, Max-Change: 0.00067
Iteration: 41, Log-Lik: -1430.495, Max-Change: 0.00017
Iteration: 42, Log-Lik: -1430.495, Max-Change: 0.00049
Iteration: 43, Log-Lik: -1430.495, Max-Change: 0.00011
Iteration: 44, Log-Lik: -1430.495, Max-Change: 0.00049
Iteration: 45, Log-Lik: -1430.495, Max-Change: 0.00015
Iteration: 46, Log-Lik: -1430.495, Max-Change: 0.00013
Iteration: 47, Log-Lik: -1430.495, Max-Change: 0.00038
Iteration: 48, Log-Lik: -1430.495, Max-Change: 0.00041
Iteration: 49, Log-Lik: -1430.495, Max-Change: 0.00019
Iteration: 50, Log-Lik: -1430.495, Max-Change: 0.00009

Unlike lavaan, mirt does not provide a nice formatting of parameters with the summary statment. Rather, we get parts of estimates through various pieces.

The model log-likelihood and summary information is given by the show() function:

show(model1PLmirt)

Call:
mirt(data = iadlData, model = mirt1PLsyntax)

Full-information item factor analysis with 1 factor(s).
Converged within 1e-04 tolerance after 38 EM iterations.
mirt version: 1.35.1 
M-step optimizer: BFGS 
EM acceleration: Ramsay 
Number of rectangular quadrature: 61
Latent density type: Gaussian 

Log-likelihood = -1443.783
Estimated parameters: 14 
AIC = 2903.566
BIC = 2939.144; SABIC = 2913.745
show(model2PLmirt)

Call:
mirt(data = iadlData, model = 1, itemtype = "2PL")

Full-information item factor analysis with 1 factor(s).
Converged within 1e-04 tolerance after 50 EM iterations.
mirt version: 1.35.1 
M-step optimizer: BFGS 
EM acceleration: Ramsay 
Number of rectangular quadrature: 61
Latent density type: Gaussian 

Log-likelihood = -1430.495
Estimated parameters: 14 
AIC = 2888.99
BIC = 2951.252; SABIC = 2906.804

Also note that the model log-likelihood information does not include a test of the model against an alternative, as does a typical CFA analysis in comparing the model fit of your model to one where all parameters were estimated. This is because the saturated model in IRT is different (for models where all items are binary, it is Multivariate Bernoulli) in that the statistics of interest come in the form of the proportion of people with a given response pattern.

To see estimates, use the coef() function. Here are the estimates for the 1PL model:

coef1PL = coef(model1PLmirt)
coef1PL
$dia1
      a1     d g u
par 4.66 1.803 0 1

$dia2
      a1     d g u
par 4.66 4.541 0 1

$dia3
      a1     d g u
par 4.66 3.509 0 1

$dia4
      a1     d g u
par 4.66 2.055 0 1

$dia5
      a1     d g u
par 4.66 1.725 0 1

$dia6
      a1     d g u
par 4.66 3.477 0 1

$dia7
      a1     d g u
par 4.66 7.292 0 1

$GroupPars
    MEAN_1 COV_11
par      0      1

The coef() function returns an R list of the parameters for each item along with the structural model parameters (the $GroupPars element), which shows the mean and variance of the latent variable. For each item, there are at least four parameters listed:

Note how the item discrimination (the a1 term) is equal for all items – this is done by convention in the 1PL model.

Putting the parameters into equation form, we have a slope/intercept form of the IRT model:

\[P(Y_{si} = 1 | \theta_s) = g_i + (u_i-g_i)\frac{\exp\left(d_i + a1_i \theta_s \right)}{1+\exp\left(d_i + a1_i \theta_s \right)}\]

Another commonly used parameterization of the IRT model is called discrimination/difficulty, given by:

\[P(Y_{si} = 1 | \theta_s) = g_i + (u_i-g_i)\frac{\exp\left(a1_i \left( \theta_s - b_i \right) \right)}{1+\exp\left(a1_i \left( \theta_s - b_i \right) \right)}\]

The two parameterizations are equivalent and one can be found by re-arranging terms of the other. To get the item difficulty from the slope/intercept parameterization:

\[b_i = -\frac{d_i}{a1_i}\]

For our results, we can use the lapply function to add the item difficulties:

getDifficulty = function(itemPar){
  parnames = colnames(itemPar)
  if ("a1" %in% parnames){
    itemPar = c(itemPar, -1*itemPar[2]/itemPar[1])
    names(itemPar) = c(parnames, "b")
    return(itemPar)
  } else {
    return(itemPar)
  }
}

lapply(X = coef1PL, FUN = getDifficulty)
$dia1
        a1          d          g          u          b 
 4.6596434  1.8033173  0.0000000  1.0000000 -0.3870076 

$dia2
        a1          d          g          u          b 
 4.6596434  4.5410128  0.0000000  1.0000000 -0.9745408 

$dia3
        a1          d          g          u          b 
 4.6596434  3.5090781  0.0000000  1.0000000 -0.7530787 

$dia4
        a1          d          g          u          b 
 4.6596434  2.0546432  0.0000000  1.0000000 -0.4409443 

$dia5
        a1          d          g          u          b 
 4.6596434  1.7248642  0.0000000  1.0000000 -0.3701709 

$dia6
        a1          d          g          u          b 
 4.6596434  3.4772527  0.0000000  1.0000000 -0.7462487 

$dia7
       a1         d         g         u         b 
 4.659643  7.292341  0.000000  1.000000 -1.565000 

$GroupPars
    MEAN_1 COV_11
par      0      1
itemPar = coef1PL[[1]]

For the 2PL, we can use a similar method (here condensed to display the item difficulties):

coef2PL = lapply(X = coef(model2PLmirt), FUN = getDifficulty)
coef2PL
$dia1
        a1          d          g          u          b 
 5.3397191  2.0120671  0.0000000  1.0000000 -0.3768114 

$dia2
        a1          d          g          u          b 
 8.7447862  8.0549708  0.0000000  1.0000000 -0.9211169 

$dia3
        a1          d          g          u          b 
 4.7378483  3.5264604  0.0000000  1.0000000 -0.7443169 

$dia4
        a1          d          g          u          b 
 5.9524105  2.5366619  0.0000000  1.0000000 -0.4261571 

$dia5
        a1          d          g          u          b 
 4.5273337  1.6601297  0.0000000  1.0000000 -0.3666904 

$dia6
        a1          d          g          u          b 
 3.3003161  2.5779920  0.0000000  1.0000000 -0.7811349 

$dia7
       a1         d         g         u         b 
 3.219578  5.396144  0.000000  1.000000 -1.676041 

$GroupPars
    MEAN_1 COV_11
par      0      1

As the 1PL is nested within the 2PL, we can use a likelihood ratio test to see which model is preferred. The LRT tests the null hypothesis that all item discriminations are equal against an alternative that not all are equal:

anova(model1PLmirt, model2PLmirt)

Model 1: mirt(data = iadlData, model = mirt1PLsyntax)
Model 2: mirt(data = iadlData, model = 1, itemtype = "2PL")

Here, the test statistic was \(\chi_6 = 18.977\) with a p-value of .004. Therefore, we reject the null hypothesis of equal slopes and conclude the 2PL fits better than the 1PL model.

The LRT, however, assumes both models have a sufficient level of absolute fit to the data. One way to tell is the use of the M2() function, which provides model fit to the 2-way tables (think item-pair covariances). Because our data has some missing responses, we have to use the impute=10 option, imputing 10 values per missing response. Here is the value for the 1PL:

M2(obj = model1PLmirt, na.rm=TRUE)
Sample size after row-wise response data removal: 614
M2(obj = model2PLmirt, na.rm = TRUE)
Sample size after row-wise response data removal: 614

The statistics given from the M2 function are similar to those used in CFA–these show approximate model fit indices such as RMSEA, SRMR, TLI, and CFI. From these, it appears the model fits approximately (CFI and TLI near 1 but relatively poor RMSEA). To find misfitting “residuals” we need complete data and the function M2() and the imputeMissing() functions are not working. So, here is an example with complete data and the 2PL:

model2PLmirtB = mirt(data = iadlData[complete.cases(iadlData),], model = 1, itemtype = "2PL")

Iteration: 1, Log-Lik: -1852.670, Max-Change: 1.46282
Iteration: 2, Log-Lik: -1518.551, Max-Change: 1.24303
Iteration: 3, Log-Lik: -1440.440, Max-Change: 1.01375
Iteration: 4, Log-Lik: -1415.396, Max-Change: 0.77349
Iteration: 5, Log-Lik: -1405.475, Max-Change: 0.64896
Iteration: 6, Log-Lik: -1400.915, Max-Change: 0.51450
Iteration: 7, Log-Lik: -1396.420, Max-Change: 0.24043
Iteration: 8, Log-Lik: -1396.101, Max-Change: 0.19007
Iteration: 9, Log-Lik: -1395.880, Max-Change: 0.15015
Iteration: 10, Log-Lik: -1395.403, Max-Change: 0.04241
Iteration: 11, Log-Lik: -1395.332, Max-Change: 0.03359
Iteration: 12, Log-Lik: -1395.271, Max-Change: 0.02179
Iteration: 13, Log-Lik: -1395.157, Max-Change: 0.01665
Iteration: 14, Log-Lik: -1395.118, Max-Change: 0.01236
Iteration: 15, Log-Lik: -1395.087, Max-Change: 0.01037
Iteration: 16, Log-Lik: -1394.953, Max-Change: 0.01294
Iteration: 17, Log-Lik: -1394.944, Max-Change: 0.01012
Iteration: 18, Log-Lik: -1394.935, Max-Change: 0.01051
Iteration: 19, Log-Lik: -1394.927, Max-Change: 0.00910
Iteration: 20, Log-Lik: -1394.921, Max-Change: 0.00891
Iteration: 21, Log-Lik: -1394.916, Max-Change: 0.00793
Iteration: 22, Log-Lik: -1394.900, Max-Change: 0.00563
Iteration: 23, Log-Lik: -1394.898, Max-Change: 0.00576
Iteration: 24, Log-Lik: -1394.896, Max-Change: 0.00618
Iteration: 25, Log-Lik: -1394.892, Max-Change: 0.00888
Iteration: 26, Log-Lik: -1394.891, Max-Change: 0.00472
Iteration: 27, Log-Lik: -1394.890, Max-Change: 0.00418
Iteration: 28, Log-Lik: -1394.889, Max-Change: 0.00408
Iteration: 29, Log-Lik: -1394.888, Max-Change: 0.00269
Iteration: 30, Log-Lik: -1394.887, Max-Change: 0.00240
Iteration: 31, Log-Lik: -1394.885, Max-Change: 0.00233
Iteration: 32, Log-Lik: -1394.885, Max-Change: 0.00152
Iteration: 33, Log-Lik: -1394.884, Max-Change: 0.00157
Iteration: 34, Log-Lik: -1394.884, Max-Change: 0.00120
Iteration: 35, Log-Lik: -1394.883, Max-Change: 0.00108
Iteration: 36, Log-Lik: -1394.883, Max-Change: 0.00075
Iteration: 37, Log-Lik: -1394.883, Max-Change: 0.00081
Iteration: 38, Log-Lik: -1394.883, Max-Change: 0.00068
Iteration: 39, Log-Lik: -1394.883, Max-Change: 0.00066
Iteration: 40, Log-Lik: -1394.883, Max-Change: 0.00020
Iteration: 41, Log-Lik: -1394.883, Max-Change: 0.00011
Iteration: 42, Log-Lik: -1394.883, Max-Change: 0.00053
Iteration: 43, Log-Lik: -1394.883, Max-Change: 0.00064
Iteration: 44, Log-Lik: -1394.883, Max-Change: 0.00017
Iteration: 45, Log-Lik: -1394.883, Max-Change: 0.00009
M2(obj = model2PLmirtB)

M2(obj = model2PLmirtB, residmat = TRUE)
            dia1          dia2         dia3       dia4         dia5       dia6 dia7
dia1          NA            NA           NA         NA           NA         NA   NA
dia2  0.00426348            NA           NA         NA           NA         NA   NA
dia3  0.06914297  1.907050e-02           NA         NA           NA         NA   NA
dia4 -0.03270413  6.323536e-03 -0.053678261         NA           NA         NA   NA
dia5  0.01037058 -6.319954e-05 -0.046104364 0.02998649           NA         NA   NA
dia6 -0.04041070 -4.516302e-02 -0.027328885 0.06404896  0.007903547         NA   NA
dia7 -0.03203898 -3.850487e-02  0.008642383 0.02416225 -0.027666584 0.08115211   NA

Here we see the biggest descripancy of residual covariances is that for dia5 with dia3 at -.08.

Finally, we can see plots of our model (all shown for the 2PL model). First, the item characteristic curves

plot(model2PLmirt, item=1, type = "trace", theta_lim = c(-3,3))

Next we can see the test information plot:

plot(model2PLmirt, type = "info", theta_lim = c(-3,3))