EPSY 906/CLDP 948 Example 5 (Updated 2022)

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:

• a1: the factor loading/item discrimination/slope: the change in the log-odds of \(P(Y_{si}=1)\) per one unit change in \(\theta_s\) (approximate change in the 3/4PL models)
• d: the item intercept the expected value of the log-odds of \(P(Y_{si}=1)\) for \(\theta_s = 0\) (approximate expectation in the 3/4PL models)
• g: the lower asymptote or guessing parameter (for 3PL this will be non-zero)
• u: the upper asymptote (for 4PL or, in mirt, a type of a 3PL, this will be non-zero)

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

We can see that our test information peaks around a theta of -.5, meaning scores near -.5 will be the most reliable.

Finally, we can use the fscores() function to get the estimated trait scores. Note, there are several types of scores available. The standard used for score reporting is method = "EAP", which are scores that use the expected value of the posterior distribution of the score. For doing secondary analyses, multiple “plausible” scores should be used with the option plausible.draws = # where # is the number of scores to draw. We plot the density of the test scores following estimating them with fscores():

theta = fscores(object = model2PLmirt, method = "EAP")
hist(theta)

The plot shows a number of people with scores at the maximum value – but very few around the most reliable portion of the test, -.5

Here, we will plot the scores along with the standard error of the scores to show the relationship:

theta2 = fscores(object = model2PLmirt, method = "EAP", full.scores.SE = TRUE)
plot(x = theta2[complete.cases(iadlData),1], y = theta2[complete.cases(iadlData),2], xlab = expression(theta), ylab = "SE", main = "Complete Data Theta vs. SE(Theta)")

For the cases with complete data, this is the best the SE gets. When plotting all the data, you can see the impact of missing data on SE: fewer observations means a higher SE for the trait.

plot(x = theta2[,1], y = theta2[,2], xlab = expression(theta), ylab = "SE", main = "All Data Theta vs. SE(Theta)")

We can also plot the item difficulty values to get a sense of the the scale is telling us about the trait:

itemDif = unlist(lapply(X = coef2PL, FUN = function(x) return(x[5])))
plot(x = 1:7, y = itemDif[1:7], type = "l", xlab = "Item", ylab = "Item Difficulty (b)")

Here are the items, again:

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%

Note that no items are available to measure above-average abilities well! The item difficulty for most items covers values of Theta between -1.0 to -0.5.

Estimation in lavaan: Limited Information Only

lavaan only provides limited information estimates of IRT/IFA models, which are parallel (but not necessarily equal to) Mplus’ WLSMV methods. This is a limitation of lavaan, so if ML versions of estimates are needed, Mplus will have to be used. Alternatively, the mirt package in R estimates IRT models (and many other types), but has rather limited capabilities for SEM with the models used and does not include model functions for continous observed variables and provides very few SEM fit statistics.

We will compare the one-parameter vs. two-parameter models for Binary Responses using WLSMV Probit model. Beginning with the two-parameter model.

Two-Parameter Model

The two-parameter model is called the two-parameter logisitic model if the logit link function is used (2PL), otherwise it is called the two-parameter normal ogive (2PNO). The syntax is largely identical to that use for CFA, however, item intercepts (item ~ 1 in CFA) are now replaced by item thresholds (item | t# where # is the number of the threshold, in order, from 1 through the number of categories on the item minus one).

The normal ogive model provides a model for a categorical variable’s (\(Y\) “underlying” continuous analog).

model2PSyntax = "

# loadings/discrimination parameters:
IADL =~ dia1 + dia2 + dia3 + dia4 + dia5 + dia6 + dia7

# threshholds use the | operator and start at value 1 after t:
dia1 | t1; dia2 | t1; dia3 | t1; dia4 | t1; dia5 | t1; dia6 | t1; dia7 | t1; 

# factor mean:
IADL ~ 0;

# factor variance:
IADL ~~ 1*IADL

dia1 ~~ dia2

"

model2PEstimates = sem(model = model2PSyntax, data = iadlData, ordered = c("dia1", "dia2", "dia3", "dia4", "dia5", "dia6", "dia7"),
                       mimic = "Mplus", estimator = "WLSMV", std.lv = TRUE, parameterization = "theta")
summary(model2PEstimates, fit.measures = TRUE, rsquare = TRUE, standardized = TRUE)

lavaan 0.6-10 ended normally after 92 iterations

  Estimator                                       DWLS
  Optimization method                           NLMINB
  Number of model parameters                        15
                                                      
                                                  Used       Total
  Number of observations                           614         631
                                                                  
Model Test User Model:
                                              Standard      Robust
  Test Statistic                                46.286      61.136
  Degrees of freedom                                13          13
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  0.797
  Shift parameter                                            3.072
       simple second-order correction (WLSMV)                     

Model Test Baseline Model:

  Test statistic                             25742.046   16653.106
  Degrees of freedom                                21          21
  P-value                                        0.000       0.000
  Scaling correction factor                                  1.546

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.999       0.997
  Tucker-Lewis Index (TLI)                       0.998       0.995
                                                                  
  Robust Comparative Fit Index (CFI)                            NA
  Robust Tucker-Lewis Index (TLI)                               NA

Root Mean Square Error of Approximation:

  RMSEA                                          0.065       0.078
  90 Percent confidence interval - lower         0.045       0.059
  90 Percent confidence interval - upper         0.085       0.098
  P-value RMSEA <= 0.05                          0.104       0.009
                                                                  
  Robust RMSEA                                                  NA
  90 Percent confidence interval - lower                        NA
  90 Percent confidence interval - upper                        NA

Standardized Root Mean Square Residual:

  SRMR                                           0.050       0.050

Weighted Root Mean Square Residual:

  WRMR                                           1.286       1.286

Parameter Estimates:

  Standard errors                           Robust.sem
  Information                                 Expected
  Information saturated (h1) model        Unstructured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  IADL =~                                                               
    dia1              3.572    0.566    6.312    0.000    3.572    0.963
    dia2              4.102    0.865    4.744    0.000    4.102    0.972
    dia3              3.036    0.453    6.708    0.000    3.036    0.950
    dia4              3.541    0.574    6.167    0.000    3.541    0.962
    dia5              2.459    0.282    8.723    0.000    2.459    0.926
    dia6              2.061    0.242    8.530    0.000    2.061    0.900
    dia7              1.744    0.294    5.927    0.000    1.744    0.867

Covariances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
 .dia1 ~~                                                               
   .dia2             -0.018    0.358   -0.049    0.961   -0.018   -0.018

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    IADL              0.000                               0.000    0.000
   .dia1              0.000                               0.000    0.000
   .dia2              0.000                               0.000    0.000
   .dia3              0.000                               0.000    0.000
   .dia4              0.000                               0.000    0.000
   .dia5              0.000                               0.000    0.000
   .dia6              0.000                               0.000    0.000
   .dia7              0.000                               0.000    0.000

Thresholds:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    dia1|t1          -1.330    0.270   -4.924    0.000   -1.330   -0.359
    dia2|t1          -3.852    0.788   -4.888    0.000   -3.852   -0.912
    dia3|t1          -2.247    0.337   -6.664    0.000   -2.247   -0.703
    dia4|t1          -1.514    0.292   -5.187    0.000   -1.514   -0.411
    dia5|t1          -0.940    0.165   -5.684    0.000   -0.940   -0.354
    dia6|t1          -1.635    0.190   -8.601    0.000   -1.635   -0.713
    dia7|t1          -2.991    0.378   -7.904    0.000   -2.991   -1.488

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    IADL              1.000                               1.000    1.000
   .dia1              1.000                               1.000    0.073
   .dia2              1.000                               1.000    0.056
   .dia3              1.000                               1.000    0.098
   .dia4              1.000                               1.000    0.074
   .dia5              1.000                               1.000    0.142
   .dia6              1.000                               1.000    0.191
   .dia7              1.000                               1.000    0.247

Scales y*:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    dia1              0.270                               0.270    1.000
    dia2              0.237                               0.237    1.000
    dia3              0.313                               0.313    1.000
    dia4              0.272                               0.272    1.000
    dia5              0.377                               0.377    1.000
    dia6              0.436                               0.436    1.000
    dia7              0.497                               0.497    1.000

R-Square:
                   Estimate
    dia1              0.927
    dia2              0.944
    dia3              0.902
    dia4              0.926
    dia5              0.858
    dia6              0.809
    dia7              0.753

We can also inspect the residual polychoric correlation matrix to investigate model misfit. Note these are the raw residuals, not a normalized or standardized version. Below that, the modification indices are displayed.

resid(model2PEstimates)

$type
[1] "raw"

$cov
     dia1   dia2   dia3   dia4   dia5   dia6   dia7  
dia1  0.000                                          
dia2  0.000  0.000                                   
dia3  0.034  0.007  0.000                            
dia4 -0.036  0.008 -0.065  0.000                     
dia5  0.002  0.007 -0.055  0.021  0.000              
dia6 -0.070 -0.044 -0.050  0.035 -0.006  0.000       
dia7 -0.135 -0.057 -0.015  0.035 -0.104  0.120  0.000

$mean
dia1 dia2 dia3 dia4 dia5 dia6 dia7 
   0    0    0    0    0    0    0 

$th
dia1|t1 dia2|t1 dia3|t1 dia4|t1 dia5|t1 dia6|t1 dia7|t1 
      0       0       0       0       0       0       0 

modificationindices(model2PEstimates, sort. = TRUE)

From the modifiation indices, we see several trends we saw with the M2() residuals in R: that items 5 and 3 need some additional help (as do items 1 and 3 and items 5 and 4).

ICC Plots

lavaanCatItemPlot = function(lavObject, varname, sds = 3){
  output = inspect(object = lavObject, what = "est")
  if (!varname %in% rownames(output$lambda)) stop(paste(varname, "not found in lavaan object"))
  if (dim(output$lambda)[2]>1) stop("plots only given for one factor models")
  
  itemloading = output$lambda[which(rownames(output$lambda) == varname),1]
  itemthresholds = output$tau[grep(pattern = varname, x = rownames(output$tau))]
  
  factorname = colnames(output$lambda)
  factormean = output$alpha[which(rownames(output$alpha) == factorname)]
  factorvar = output$psi[which(rownames(output$psi) == factorname)]
  
  factormin = factormean - 3*sqrt(factorvar)
  factormax = factormean + 3*sqrt(factorvar)
  
  factorX = seq(factormin, factormax, .01)
  itemloc = which(lavObject@Data@ov$name == varname)      
  itemlevels = unlist(strsplit(x = lavObject@Data@ov$lnam[itemloc], split = "\\|"))  
  if (length(itemthresholds)>1){

    
    plotdata = NULL
    plotdata2 = NULL
    itemY = NULL
    itemY2 = NULL
    itemX = NULL
    itemText = NULL
    for (level in 1:length(itemthresholds)){
      
      itemY = pnorm(q = -1*itemthresholds[level] + itemloading*factorX)
      itemY2 = cbind(itemY2, pnorm(q = -1*itemthresholds[level] + itemloading*factorX))
      itemText = paste0("P(", varname, " > ", itemlevels[level], ")")
      itemText2 = paste0("P(", varname, " = ", itemlevels[level], ")")
      plotdata = rbind(plotdata, data.frame(factor = factorX, prob = itemY, plot = itemText))
      
      if (level == 1){
        plotdata2 = data.frame(factor = factorX, plot = itemText2, prob = matrix(1, nrow = dim(itemY2)[1], ncol=1) - itemY2[,level])
      } else if (level == length(itemthresholds)){
        plotdata2 = rbind(plotdata2, data.frame(factor = factorX, plot = itemText2, prob = itemY2[,level-1] - itemY2[,level]))
        plotdata2 = rbind(plotdata2, data.frame(factor = factorX, plot = paste0("P(", varname, " = ", itemlevels[level+1], ")"), prob = itemY2[,level]))                  
      } else {
        plotdata2 = rbind(plotdata2, data.frame(factor = factorX, plot = itemText2, prob = itemY2[,level-1] - itemY2[,level]))
      }
      
    }
    
    names(plotdata) = c(factorname , "Probability", "Cumulative")
    ggplot(data = plotdata, aes_string(x = factorname, y = "Probability", colour = "Cumulative")) + geom_line(size = 2)
    
    names(plotdata2) = c(factorname, "Response", "Probability")
    ggplot(data = plotdata2, aes_string(x = factorname, y = "Probability", colour = "Response")) + geom_line(size = 2)
  } else {
    
    itemY = pnorm(q = -1*itemthresholds[1] + itemloading*factorX)
    itemText2 = paste0("P(", varname, " = ", itemlevels[2], ")")
    plotdata = data.frame(factor = factorX, prob = itemY, plot = itemText2)
    
    names(plotdata) = c(factorname , "Probability", "Response")
    ggplot(data = plotdata, aes_string(x = factorname, y = "Probability", colour = "Response")) + geom_line(size = 2)
    
  }
}

lavaanCatItemPlot(lavObject = model2PEstimates, varname = "dia2", sds = 3)

Conversion to IRT Paramterization (Discrimination/Difficulty)

convertTheta2IRT = function(lavObject){
  if (!lavObject@Options$parameterization == "theta") stop("your model is not estimated with parameterization='theta'")
  
  output = inspect(object = lavObject, what = "est")
  if (ncol(output$lambda)>1) stop("IRT conversion is only valid for one dimensional factor models. Your model has more than one dimension.")
  
  a = output$lambda
  b = output$tau/output$lambda
  return(list(a = a, b=b))
}

convertTheta2IRT(lavObject = model2PEstimates)

$a
      IADL
dia1 3.572
dia2 4.102
dia3 3.036
dia4 3.541
dia5 2.459
dia6 2.061
dia7 1.744

$b
        thrshl
dia1|t1 -0.372
dia2|t1 -0.939
dia3|t1 -0.740
dia4|t1 -0.427
dia5|t1 -0.382
dia6|t1 -0.793
dia7|t1 -1.715

One-Parameter Normal Ogive Model

To estimate the 1PNO model in lavaan, we use the following syntax. The summary, residuals, and modification indices are displayed subsequently.

model1PSyntax = "

# loadings/discrimination parameters:
IADL =~ loading*dia1 + loading*dia2 + loading*dia3 + loading*dia4 + loading*dia5 + loading*dia6 + loading*dia7

# threshholds use the | operator and start at value 1 after t:
dia1 | t1; dia2 | t1; dia3 | t1; dia4 | t1; dia5 | t1; dia6 | t1; dia7 | t1; 

# factor mean:
IADL ~ 0;

# factor variance:
IADL ~~ 1*IADL

"

model1PEstimates = sem(model = model1PSyntax, data = iadlData, ordered = c("dia1", "dia2", "dia3", "dia4", "dia5", "dia6", "dia7"),
                       mimic = "Mplus", estimator = "WLSMV", std.lv = TRUE, parameterization = "theta")
summary(model1PEstimates, fit.measures = TRUE, rsquare = TRUE, standardized = TRUE)

lavaan 0.6-10 ended normally after 26 iterations

  Estimator                                       DWLS
  Optimization method                           NLMINB
  Number of model parameters                        14
  Number of equality constraints                     6
                                                      
                                                  Used       Total
  Number of observations                           614         631
                                                                  
Model Test User Model:
                                              Standard      Robust
  Test Statistic                                83.736      80.544
  Degrees of freedom                                20          20
  P-value (Chi-square)                           0.000       0.000
  Scaling correction factor                                  1.107
  Shift parameter                                            4.921
       simple second-order correction (WLSMV)                     

Model Test Baseline Model:

  Test statistic                             25742.046   16653.106
  Degrees of freedom                                21          21
  P-value                                        0.000       0.000
  Scaling correction factor                                  1.546

User Model versus Baseline Model:

  Comparative Fit Index (CFI)                    0.998       0.996
  Tucker-Lewis Index (TLI)                       0.997       0.996
                                                                  
  Robust Comparative Fit Index (CFI)                            NA
  Robust Tucker-Lewis Index (TLI)                               NA

Root Mean Square Error of Approximation:

  RMSEA                                          0.072       0.070
  90 Percent confidence interval - lower         0.057       0.055
  90 Percent confidence interval - upper         0.088       0.087
  P-value RMSEA <= 0.05                          0.011       0.018
                                                                  
  Robust RMSEA                                                  NA
  90 Percent confidence interval - lower                        NA
  90 Percent confidence interval - upper                        NA

Standardized Root Mean Square Residual:

  SRMR                                           0.071       0.071

Weighted Root Mean Square Residual:

  WRMR                                           1.729       1.729

Parameter Estimates:

  Standard errors                           Robust.sem
  Information                                 Expected
  Information saturated (h1) model        Unstructured

Latent Variables:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
  IADL =~                                                               
    dia1    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947
    dia2    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947
    dia3    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947
    dia4    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947
    dia5    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947
    dia6    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947
    dia7    (ldng)    2.961    0.187   15.849    0.000    2.961    0.947

Intercepts:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    IADL              0.000                               0.000    0.000
   .dia1              0.000                               0.000    0.000
   .dia2              0.000                               0.000    0.000
   .dia3              0.000                               0.000    0.000
   .dia4              0.000                               0.000    0.000
   .dia5              0.000                               0.000    0.000
   .dia6              0.000                               0.000    0.000
   .dia7              0.000                               0.000    0.000

Thresholds:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    dia1|t1          -1.121    0.177   -6.341    0.000   -1.121   -0.359
    dia2|t1          -2.851    0.198  -14.376    0.000   -2.851   -0.912
    dia3|t1          -2.197    0.195  -11.286    0.000   -2.197   -0.703
    dia4|t1          -1.286    0.180   -7.125    0.000   -1.286   -0.411
    dia5|t1          -1.107    0.174   -6.348    0.000   -1.107   -0.354
    dia6|t1          -2.230    0.199  -11.206    0.000   -2.230   -0.713
    dia7|t1          -4.650    0.299  -15.532    0.000   -4.650   -1.488

Variances:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    IADL              1.000                               1.000    1.000
   .dia1              1.000                               1.000    0.102
   .dia2              1.000                               1.000    0.102
   .dia3              1.000                               1.000    0.102
   .dia4              1.000                               1.000    0.102
   .dia5              1.000                               1.000    0.102
   .dia6              1.000                               1.000    0.102
   .dia7              1.000                               1.000    0.102

Scales y*:
                   Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
    dia1              0.320                               0.320    1.000
    dia2              0.320                               0.320    1.000
    dia3              0.320                               0.320    1.000
    dia4              0.320                               0.320    1.000
    dia5              0.320                               0.320    1.000
    dia6              0.320                               0.320    1.000
    dia7              0.320                               0.320    1.000

R-Square:
                   Estimate
    dia1              0.898
    dia2              0.898
    dia3              0.898
    dia4              0.898
    dia5              0.898
    dia6              0.898
    dia7              0.898

resid(object = model1PEstimates)

$type
[1] "raw"

$cov
     dia1   dia2   dia3   dia4   dia5   dia6   dia7  
dia1  0.000                                          
dia2  0.037  0.000                                   
dia3  0.051  0.032  0.000                            
dia4 -0.007  0.046 -0.048  0.000                     
dia5 -0.003  0.010 -0.073  0.015  0.000              
dia6 -0.101 -0.068 -0.094  0.003 -0.070  0.000       
dia7 -0.197 -0.112 -0.089 -0.028 -0.198  0.003  0.000

$mean
dia1 dia2 dia3 dia4 dia5 dia6 dia7 
   0    0    0    0    0    0    0 

$th
dia1|t1 dia2|t1 dia3|t1 dia4|t1 dia5|t1 dia6|t1 dia7|t1 
      0       0       0       0       0       0       0 

modificationindices(model1PEstimates, sort. = TRUE)

To compare models, we use the anova() function, which uses the correct test for us. Here we see the 2PL is preferred to the 1PL, again.

anova(model1PEstimates, model2PEstimates)

Scaled Chi-Squared Difference Test (method = “satorra.2000”)

lavaan NOTE:
    The “Chisq” column contains standard test statistics, not the
    robust test that should be reported per model. A robust difference
    test is a function of two standard (not robust) statistics.
 
                 Df AIC BIC  Chisq Chisq diff Df diff Pr(>Chisq)    
model2PEstimates 13         46.286                                  
model1PEstimates 20         83.737     27.913       7  0.0002281 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1