Overview

This document demonstrates an IRT analysis using Mplus where the outcome variables (items) are ordered categorical variables. The example uses Apgar data from the CPP (collaborative perinatal project) and are publicly available (ICPSR). We demonstrate how to use an R package for running Mplus from R (MplusAutomation). In this example, use the default Mplus estimator for ordinal dependent variables (WLSMV estimator, delta parameterization) which implements a multivariate probit model. We generate boundary response functions, expected item response functions, test response functions, and item characteristic curves or category response functions.

Model and estimator

This set of examples uses the default estimator options in Mplus, which for ordinal dependent variables is weighted least squares with mean and variance standardization (WLSMV), a probit link, and so-called delta parameterization. This is not my favorite estimator for categorical dependent variables using Mplus, for reasons I describe elsewhere. But there is nothing wrong with this estimator. The reason why I don’t like it is the parameter estimates are on a scale that is different from the 2PL model, and this makes it harder to compare results or switch from one estimator (maximum likelihood, Bayes) and be able to compare results.

Predicted probabilities with Mplus WLSMV delta

When we run an IRT-like model with Mplus (a single factor confirmatory factor analysis model with categorical dependent variables), the model estimates include measurement slopes (\(\lambda\)), item thresholds (\(\tau\)), and a variance parameter for the underlying latent trait (\(\psi\)).

Binary dependent variables

In the case of binary dependent variables, the probability of observing a response to item \(u_j\) in category 1 vs category 0 given an estimate of the underlying latent trait (\(\eta\)) is, using the parameters from a Mplus WLSMV/delta probit model:

\[P(u_{ij}=1|\eta_i) = \Phi^{-1}\left(\frac{-\tau_{1j} + \lambda_j\eta_i}{\sqrt{1-\lambda_j^2}}\right)\]

Where \(j\) is indexing items and subscript \(i\) identifies variables that vary over individual people; \(u_{ij}\) is an item response to binary test item \(j\) by person \(i\), (\(u\) is observed as 0 or 1), \(\eta_i\) is a person \(i\)’s level of a latent trait or latent ability, \(\lambda_j\) is the measurement slope for item \(j\), \(\tau_{1j}\) is threshold 1 for item \(j\), and \(\Phi^{-1}\) is a link function, in this case a normal probability or probit link function.

Outcome variables with more than 2 categories: introducing the latent response variable framework

In the case of ordered categorical outcomes, the function is a little bit more complicated. In the rest of this subsection I will be paraphrasing from Mplus Web Note 4.

We will consider the case of an outcome variable that is observed in ordered categories. The variable has \(C\) categories, and we’ll denote each category with letter \(c\) (and \(c = 0, 1, 2, ..., C-1\)). In the Mplus WLS approach to these data, we use the latent response variable formulation to model the relationships among the observed categorical dependent variables \(u\). A latent response variable, that we will represent as \(y^*\) (read that as “y-star”), is an unobserved, continuously distributed, and presumably normally distributed variable that we assume exists. Think of a \(y_j^*\) as a propensity of a person to respond in a higher category on the observed dependent variable item \(u_j\). This latent response variable \(y^*\) is not the same thing as the common factor latent variable (\(\eta\)): \(\eta\) represents the shared variance of multiple \(y^*\) variables and the tendency to respond in a higher category across multiple latent response variables. It’s also not the same thing as \(\epsilon_{ij}\), which is also latent and also related to the probability that a person has a higher level on the latent response variable (and therefore responds in a higher category on the observed categorical dependent variable): \(\epsilon\) is independent from \(\eta\). Here’s an equation for the latent response variable framework:

\[\begin{align} y^*_i & = \nu + \lambda \eta_i + \epsilon_i \\ E(y^*) & = \mu^* = \nu + \lambda \alpha \\ V(y^*) & = \sigma^* = \lambda^2\psi + \theta \\ \end{align}\]

where \(\nu\) is the intercept, \(\lambda\) is the measurement slope or factor loading, \(\eta_i\) is the common latent variable or trait or ability of interest, and \(\epsilon_i\) is a residual or error term. This is an ordinary linear regression model, other than the fact that all of the variables in the model latent. The expectation of \(y^*\) is determined by the mean of the common latent variable \(\eta\), which is represented by \(\alpha\) in the equation above. The variance of \(y^*\) is a function of the factor loading, the variance of the latent variable \(\eta\) represented with \(\psi\), and the residual variance of the latent response variable (i.e., the variance of the residual \(\epsilon\)), represented with \(\theta\).

Now that we have a model for the distribution of a continuous latent response variable (\(y^* \sim N(\mu^*,\sigma^*)\)), we need to map that latent response variable according to the categories we see in the observed categorical response variable \(u\). This is accomplished by identifying thresholds (\(\tau\)) on the latent response variable such that if we cut the distribution at the thresholds, the proportion of people so classified would match the observed proportion in the observed categories. Because we use the thresholds as cut points, we will have one fewer thresholds than there are response categories. (It takes one threshold to cut a continuous distribution into two groups; it takes two thresholds to cut a continuous distribution into three groups, etc.). For an individual participant, their most likely response category to item \(u\) is determined by how their estimated \(y^*\) value underlying that \(u\) compares to the estimated thresholds \(\tau\):

\[ u = c, \; \text{if} \; \tau_c < y^* \le \tau_{c+1} \]

for categories \(c = 0, 1, 2, \ldots, C-1\), and where \(\tau_0 = - \infty\) and \(\tau_C = + \infty\). From this we can represent a conditional probability expression

\[P(u_{ij} \ge c|\eta_i) = \Phi^{-1}\left(\frac{-\tau_{cj} + \nu + \lambda_j\eta_i}{\sqrt{\theta}}\right)\] where, remember, \(\theta\) is the variance of the residual \(\epsilon\).

Identifying assumptions that are involved in the Mplus WLSMV estimator, probit link, and delta parameterization (the Mplus default with categorical dependent variables) are as follows. First, the probit model informs that \(\Phi\) is a normal function. It is typical to standardize the model such that the expectation of the latent variable \(\eta\) is 0, so \(\alpha = 0\), and therefore \(\mu^* = \nu\). It is also common to standardize such that \(\nu = 0\), and \(\sigma^* = 1\). With these standardization assumptions, the variance of the residual is given as a remainder

\[ \theta = 1 - \lambda^2\psi\]

or \(\theta = 1-\lambda^2\) if \(\psi = 1\). So if \(\nu = 0\) and \(\psi = 1\), the conditional probability expression reduces to:

\[P(u_{ij} \ge c|\eta_i) = \Phi^{-1}\left(\frac{-\tau_{cj} + \lambda_j\eta_i}{\sqrt{1-\lambda^2}}\right)\] .

The reason why we divide \((-\tau + \lambda\eta)\) by \(\sqrt{1-\lambda^2}\) when the model is multivariate probit under the delta parameterization is explained in Mplus Web Note 4. The value \((1-\lambda^2\psi)\) (where \(\psi\) is the common latent variable, \(\eta\), variance or residual variance; I have left it out of the equation above because in our example, \(\psi=1\)) is the residual variance of the latent response variable (the continuous \(y^*\) variable underlying the observed categorical \(y\) variable). Dividing WLSMV delta parameterization estimates in \(\tau\) and \(\lambda\) by the standard deviation of \(y^*\) (i.e., the square root of the variance of \(y^*\)) is necessary to use the standard form of the normal probability CDF (pnorm). In the theta parameterization, the theta’s are parameters (not solved by \(1-\lambda^2\psi\)) and therefore the delta parameters are constrained (to 1 in the single group case as a default). The delta parameters are the inverse variance terms for the latent response variable (\(\Delta = 1/\sqrt{\sigma^*}\)).

Trace line functions

The function \(P(u_{ij} \ge c|\eta_i)\) is of interest, and can be used to prepare boundary response functions. These are plots that show the probability of responding in category \(c\) or higher as a function of the underlying latent trait (\(\eta\)). Each item will have \(C-1\) trace lines in a boundary response function, one for each threshold.

Item characteristic curves or more descriptively, category response functions, represent the probability or responding in category \(c\) as a function of the underlying latent trait \(\eta\). In these plots there will be one trace line for each response category. If \(c = 0, 1, 2, \ldots, C-1\), and \(P(u_{ij} \ge 1|\eta_i)\) is the probability of being in category \(c = 1\) or higher, then the probability of being in category \(c=0\) is \(1-P(u_{ij} \ge 1|\eta_i)\) (i.e., the probability of not being in category 1 or higher). The probability of responding in category 1 is the probability of responding in category 1 or higher (\(P(u_{ij} \ge 1|\eta_i)\)) minus the probability of responding in category 2 or higher (\(P(u_{ij} \ge 2|\eta_i)\)).

\[P(u_{ij} = c|\eta_i) = \Phi^{-1}\left(\frac{-\tau_{cj} + \lambda_j\eta_i}{\sqrt{1-\lambda^2}}\right) - \Phi^{-1}\left(\frac{-\tau_{(c+1)j} + \lambda_j\eta_i}{\sqrt{1-\lambda^2}}\right)\] .

Finally, I often also like to generate expected item response functions, which plot the expected value of the item response (on a scale from 0 to \(C\)) as a function of the underlying latent variable \(\eta\). These functions are simply the sums of the boundary response functions across all thresholds for an item. The test response function is the sum of all expected item response functions for all items in a test, and plots the expected sum score on the test as a function of \(\eta\).

Applied Example

Load packages and libraries

#install.packages("lavaan")
#install.packages("MplusAutomation")
#install.packages("car") 
#install.packages("latex2exp")
#install.packages("gmodels")
#install.packages("haven")
#install.packages("tidyverse")
library(haven) # for reading dta file
library(tidyverse) # for everything
library(MplusAutomation) # for IRT
library(car) # friendly recoding
library(latex2exp) # for including math in ggplot2 labels
library(gmodels) # for a variant of tables (CrossTables)

Read data

apgar <- haven::read_dta("https://s3.amazonaws.com/quantsci/QSP08/apgar.dta")
head(apgar)[,c("cps_46","cps_47","cps_48","cps_49","cps_50","cps_51")]
## # A tibble: 6 × 6
##             cps_46             cps_47           cps_48   cps_49  cps_50   cps_51
##          <dbl+lbl>          <dbl+lbl>        <dbl+lbl> <dbl+lb> <dbl+l> <dbl+lb>
## 1 2 [100 or more]  2 [Crying Lustily] 2 [Well flexed]  2 [Cry]  1 [Blu…  9      
## 2 2 [100 or more]  2 [Crying Lustily] 2 [Well flexed]  2 [Cry]  2 [Ent… 10      
## 3 2 [100 or more]  2 [Crying Lustily] 2 [Well flexed]  2 [Cry]  1 [Blu…  9      
## 4 2 [100 or more]  2 [Crying Lustily] 2 [Well flexed]  2 [Cry]  2 [Ent… 10      
## 5 2 [100 or more]  2 [Crying Lustily] 2 [Well flexed]  2 [Cry]  1 [Blu…  9      
## 6 9 [Not Reported] 9 [Not Reported]   9 [Not Reported] 9 [Not … 9 [Not… 99 [Not…

Make variables

# pull apgar items out into a separate data frame
apgar2 <- apgar[,c("cps_46","cps_47","cps_48","cps_49","cps_50","cps_51","g1_id","g2_id","famid")]
# handle the missing values
apgar2$u1 <- car::recode(apgar2$cps_46,"9=NA;99=NA")
apgar2$u2 <- car::recode(apgar2$cps_47,"9=NA;99=NA")
apgar2$u3 <- car::recode(apgar2$cps_48,"9=NA;99=NA")
apgar2$u4 <- car::recode(apgar2$cps_49,"9=NA;99=NA")
apgar2$u5 <- car::recode(apgar2$cps_50,"9=NA;99=NA")
apgar2$apgartot <- car::recode(apgar2$cps_51,"99=NA;999=NA")
# table(apgar2$cps_51 , useNA="always")
# table(apgar2$apgartot , useNA="always")

# Also need to force the variables to be numeric
# otherwise MplusAutomation has trouble
apgar2$u1 <- as.numeric(apgar2$u1)
apgar2$u2 <- as.numeric(apgar2$u2)
apgar2$u3 <- as.numeric(apgar2$u3)
apgar2$u4 <- as.numeric(apgar2$u4)
apgar2$u5 <- as.numeric(apgar2$u5)

Select cases for analysis.

When we used mirt, we had to restrict to listwise complete data. We do not need to do that with Mplus

#apgar2$dropme <- is.na(apgar2$u1) & is.na(apgar2$u2) & is.na(apgar2$u3) & is.na(apgar2$u4) & is.na(apgar2$u5)
#apgar2i <- apgar2[which(apgar2$dropme==0),]
apgar2i <- apgar2

Generate descriptive statitiscs

table(apgar2i$u1 , useNA="always")
## 
##    0    1    2 <NA> 
##    1  111 2151  263
table(apgar2i$u2 , useNA="always")
## 
##    0    1    2 <NA> 
##   55  436 1788  247
table(apgar2i$u3 , useNA="always")
## 
##    0    1    2 <NA> 
##   87  399 1792  248
table(apgar2i$u4 , useNA="always")
## 
##    0    1    2 <NA> 
##   71  267 1941  247
table(apgar2i$u5 , useNA="always")
## 
##    0    1    2 <NA> 
##  359 1687  234  246

Do IRT (using R/MplusAutomation)

# mod1 : name of object that will hold the results of the Mplus model
# mo : Mplus object (the instructions to send to Mplus)


mo <- MplusAutomation::mplusObject(
   VARIABLE = "CATEGORICAL ARE all;" ,
   MODEL = "f by u1-u5*; f@1" ,
   usevariables = c("u1","u2","u3","u4","u5") ,
   rdata = apgar2i
)


mod1 <- MplusAutomation::mplusModeler(mo, modelout = "mod1.inp", run=TRUE)

Estimated coefficients.

In this step we’re going to read the save Mplus output file (mod1.out) and extract the parameter estimates. You could probably start with this step, if for example you had run Mplus interactively separately from R and R/MplusAutomation.

# Mplus model parameter estimates
param <- MplusAutomation::readModels("mod1.out")$parameters$unstandardized %>% as_tibble()

pwd <- param %>%
  mutate(param = tolower(param))

pwd
## # A tibble: 16 × 6
##    paramHeader param    est    se est_se  pval
##    <chr>       <chr>  <dbl> <dbl>  <dbl> <dbl>
##  1 F.BY        u1     0.858 0.022   39.0     0
##  2 F.BY        u2     0.933 0.01    93.1     0
##  3 F.BY        u3     0.919 0.01    92.6     0
##  4 F.BY        u4     0.937 0.01    96.3     0
##  5 F.BY        u5     0.516 0.024   21.5     0
##  6 Thresholds  u1$1  -3.32  0.279  -11.9     0
##  7 Thresholds  u1$2  -1.65  0.045  -37.0     0
##  8 Thresholds  u2$1  -1.98  0.057  -34.9     0
##  9 Thresholds  u2$2  -0.788 0.029  -26.8     0
## 10 Thresholds  u3$1  -1.77  0.048  -36.6     0
## 11 Thresholds  u3$2  -0.795 0.03   -26.9     0
## 12 Thresholds  u4$1  -1.86  0.052  -36.0     0
## 13 Thresholds  u4$2  -1.04  0.032  -32.4     0
## 14 Thresholds  u5$1  -1.00  0.032  -31.7     0
## 15 Thresholds  u5$2   1.27  0.036   35.6     0
## 16 Variances   f      1     0      999     999

In this next step I am going to put the measurement slopes and thresholds in separate lists. I use these later to make my plots.

# lambda (l) parameters as a list
l1 <-  pwd$est[which(pwd$paramHeader=="F.BY" & pwd$param=="u1")]
l2 <-  pwd$est[which(pwd$paramHeader=="F.BY" & pwd$param=="u2")]
l3 <-  pwd$est[which(pwd$paramHeader=="F.BY" & pwd$param=="u3")]
l4 <-  pwd$est[which(pwd$paramHeader=="F.BY" & pwd$param=="u4")]
l5 <-  pwd$est[which(pwd$paramHeader=="F.BY" & pwd$param=="u5")]
l <- c(l1,l2,l3,l4,l5)

# t parameters as seperate lists (thresholds)
t11 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u1$1")]
t12 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u2$1")]
t13 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u3$1")]
t14 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u4$1")]
t15 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u5$1")]
t1 <- c(t11,t12,t13,t14,t15)

t21 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u1$2")]
t22 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u2$2")]
t23 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u3$2")]
t24 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u4$2")]
t25 <-  pwd$est[which(pwd$paramHeader=="Thresholds" & pwd$param=="u5$2")]
t2 <- c(t21,t22,t23,t24,t25)

ICCs by hand

These are going to be boundary response functions, since there will be more than one curve per item (one for each threshold) I used color brewer to choose my colors (http://colorbrewer2.org/#type=qualitative&scheme=Set1&n=5)

pwdFunction <- function(x,l,t) { 
   pnorm((-1*t + l*x)/((1-l^2)^.5))
} 
brf <- ggplot(data.frame(x = c(-6, 6)), aes(x = x)) + 
  stat_function(fun = pwdFunction, args = list(l=l[[1]],t=t1[[1]]) , aes(colour = "1-heart rate")) + # red
  stat_function(fun = pwdFunction, args = list(l=l[[1]],t=t2[[1]]) , aes(colour = "1-heart rate")) +
  stat_function(fun = pwdFunction, args = list(l=l[[2]],t=t1[[2]]) , aes(colour = "2-respiration")) + # orange
  stat_function(fun = pwdFunction, args = list(l=l[[2]],t=t2[[2]]) , aes(colour = "2-respiration")) +
  stat_function(fun = pwdFunction, args = list(l=l[[3]],t=t1[[3]]) , aes(colour = "3-muscle tone")) + # green
  stat_function(fun = pwdFunction, args = list(l=l[[3]],t=t2[[3]]) , aes(colour = "3-muscle tone")) +
  stat_function(fun = pwdFunction, args = list(l=l[[4]],t=t1[[4]]) , aes(colour = "4-reflex irr")) + # purple
  stat_function(fun = pwdFunction, args = list(l=l[[4]],t=t2[[4]]) , aes(colour = "4-reflex irr")) +
  stat_function(fun = pwdFunction, args = list(l=l[[5]],t=t1[[5]]) , aes(colour = "5-color")) + # blue
  stat_function(fun = pwdFunction, args = list(l=l[[5]],t=t2[[5]]) , aes(colour = "5-color")) +
  scale_colour_manual("Apgar item", values = c("#e41a1c", "#ff7f00", "#4daf4a", "#984ea3" , "#377eb8")) +
  theme(legend.position="bottom") +
  theme(legend.text = element_text(size=10)) +
  scale_x_continuous(name= TeX("Ability $(\\eta)$")) +
  scale_y_continuous(name= TeX("$P(u \\geq c|\\eta)$"))

brf

Change the order

Instead of assigning the color to the label indicating the content of the item, I assign the color to an alphanumeric label that captures the order I would like to see the lines plotted. Then, I use the labels function in the scale_color_manual option to re-assign the content label to the alphanumeric label. color to

brf2 <- ggplot(data.frame(x = c(-6, 6)), aes(x = x)) + 
  stat_function(fun = pwdFunction, args = list(l=l[[1]],t=t1[[1]]) , aes(color = "item1")) + # heart rate, red = e41a1c
  stat_function(fun = pwdFunction, args = list(l=l[[1]],t=t2[[1]]) , aes(color = "item1")) +
  stat_function(fun = pwdFunction, args = list(l=l[[2]],t=t1[[2]]) , aes(color = "item2")) + # respiration, orange = ff7f00
  stat_function(fun = pwdFunction, args = list(l=l[[2]],t=t2[[2]]) , aes(color = "item2")) +
  stat_function(fun = pwdFunction, args = list(l=l[[3]],t=t1[[3]]) , aes(color = "item3")) + # muscle tone, green = 4daf4a
  stat_function(fun = pwdFunction, args = list(l=l[[3]],t=t2[[3]]) , aes(color = "item3")) +
  stat_function(fun = pwdFunction, args = list(l=l[[4]],t=t1[[4]]) , aes(color = "item4")) + # reflex irr , purple = 984ea3
  stat_function(fun = pwdFunction, args = list(l=l[[4]],t=t2[[4]]) , aes(color = "item4")) +
  stat_function(fun = pwdFunction, args = list(l=l[[5]],t=t1[[5]]) , aes(color = "item5")) + # color , blue = 377eb8
  stat_function(fun = pwdFunction, args = list(l=l[[5]],t=t2[[5]]) , aes(color = "item5")) +
  scale_color_manual(name = "Apgar item", 
    values = c("item1" = "#e41a1c",  "item2" = "#ff7f00", "item3" = "#4daf4a", "item4" = "#984ea3" , "item5" = "#377eb8" ) , 
    labels = c("item1" = "heart rate" , "item2" = "respiration" , "item3" = "muscle tone" , "item4" = "reflex irr" , "item5" = "color")) +
  theme(legend.position="bottom") + 
  theme(legend.text = element_text(size=10)) +
  scale_x_continuous(name= TeX("Ability $(\\eta)$")) +
  scale_y_continuous(name= TeX("$P(u \\geq c|\\eta)$"))


brf2

Plot the Expected Item Response Function

pwdEIRFunction <- function(x,l,t1,t2) { 
  p1 <- pnorm((-1*t1 + l*x)/((1-l^2)^.5))
  p2 <- pnorm((-1*t2 + l*x)/((1-l^2)^.5))
  p1 + p2
} 
eir <- ggplot(data.frame(x = c(-6, 6)), aes(x = x)) + 
  stat_function(fun = pwdEIRFunction, args = list(l=l[[1]],t1=t1[[1]],t2=t2[[1]]) , aes(colour = "1-heart rate")) + # red
  stat_function(fun = pwdEIRFunction, args = list(l=l[[2]],t1=t1[[2]],t2=t2[[2]]) , aes(colour = "2-respiration")) + # orange
  stat_function(fun = pwdEIRFunction, args = list(l=l[[3]],t1=t1[[3]],t2=t2[[3]]) , aes(colour = "3-muscle tone")) + # green
  stat_function(fun = pwdEIRFunction, args = list(l=l[[4]],t1=t1[[4]],t2=t2[[4]]) , aes(colour = "4-reflex irr")) + # purple
  stat_function(fun = pwdEIRFunction, args = list(l=l[[5]],t1=t1[[5]],t2=t2[[5]]) , aes(colour = "5-color")) + # blue
  scale_colour_manual("Apgar item", values = c("#e41a1c", "#ff7f00", "#4daf4a", "#984ea3" , "#377eb8")) +
  theme(legend.position="bottom") +
  theme(legend.text = element_text(size=10)) +
  scale_x_continuous(name= TeX("Ability $(\\eta)$")) +
  scale_y_continuous(name= TeX("$E(u|\\eta)$"))

eir

Plot a test response function

pwdTRFunction <- function(x,
                          l1,t11,t21,
                          l2,t12,t22,
                          l3,t13,t23,
                          l4,t14,t24,
                          l5,t15,t25
                          ) { 
  p11 <- pnorm((-1*t11 + l1*x)/((1-l1^2)^.5))
  p12 <- pnorm((-1*t12 + l2*x)/((1-l2^2)^.5))
  p13 <- pnorm((-1*t13 + l3*x)/((1-l3^2)^.5))
  p14 <- pnorm((-1*t14 + l4*x)/((1-l4^2)^.5))
  p15 <- pnorm((-1*t15 + l5*x)/((1-l5^2)^.5))
  p21 <- pnorm((-1*t21 + l1*x)/((1-l1^2)^.5))
  p22 <- pnorm((-1*t22 + l2*x)/((1-l2^2)^.5))
  p23 <- pnorm((-1*t23 + l3*x)/((1-l3^2)^.5))
  p24 <- pnorm((-1*t24 + l4*x)/((1-l4^2)^.5))
  p25 <- pnorm((-1*t25 + l5*x)/((1-l5^2)^.5))
  p11+p12+p13+p14+p15+p21+p22+p23+p24+p15
} 
trf <- ggplot(data.frame(x = c(-6, 6)), aes(x = x)) + 
  stat_function(fun = pwdTRFunction, args = list(
                                            l1=l[[1]],t11=t1[[1]],t21=t2[[1]],
                                            l2=l[[2]],t12=t1[[2]],t22=t2[[2]],                                                 
                                            l3=l[[3]],t13=t1[[3]],t23=t2[[3]],
                                            l4=l[[4]],t14=t1[[4]],t24=t2[[4]],                                                 
                                            l5=l[[5]],t15=t1[[5]],t25=t2[[5]])) +
  scale_x_continuous(name= TeX("Ability $(\\eta)$")) +
  scale_y_continuous(name= TeX("$E(\\Sigma U|\\eta)$"))

trf

Item characteristic curves

Item characteristic curves for ordered categorical items are displayed with the number of lines equal to the number of response categories. The curves can also be called category response functions. I will only show the item characteristic curve for the first item.

# In the function below (pwd.ordinalicc0), p1 is the probability of responding in category 1 or
# higher. We've defined it before in the boundary response function. But
# with the last line of the function, what the function generates is 
# 1-p1, which is the probably of responding in category 0 (i.e., not
# 1 or higher).

pwd.ordinalicc0 <- function(x,l,t1,t2) { 
  p1 <- pnorm((-1*t1 + l*x)/((1-l^2)^.5))
  1-p1
} 

# In the function below (pwd.ordinalicc1), we compute the probability of responding
# in category 1. This is the probability of responding in category 1 or higher, 
# minus the probability of responding in category 2 or higher.

pwd.ordinalicc1 <- function(x,l,t1,t2) { 
  p1 <- pnorm((-1*t1 + l*x)/((1-l^2)^.5))
  p2 <- pnorm((-1*t2 + l*x)/((1-l^2)^.5))
  p1-p2
} 

# pwd.ordinalicc2 returns the probability of responding in the
# absorbing category, the top category.

pwd.ordinalicc2 <- function(x,l,t1,t2) { 
  pnorm((-1*t2 + l*x)/((1-l^2)^.5))
} 


crf <- ggplot(data.frame(x = c(-6, 0)), aes(x = x)) + 
  stat_function(fun = pwd.ordinalicc0, args = list(l=l[[1]],t1=t1[[1]],t2=t2[[1]]) , aes(colour = "category 0")) + # red
  stat_function(fun = pwd.ordinalicc1, args = list(l=l[[1]],t1=t1[[1]],t2=t2[[1]]) , aes(colour = "category 1")) + # orange
  stat_function(fun = pwd.ordinalicc2, args = list(l=l[[1]],t1=t1[[1]],t2=t2[[1]]) , aes(colour = "category 2")) + # green
  scale_colour_manual("Apgar item 1 Heart Rate", values = c("#e41a1c","#ff7f00", "#4daf4a")) +
  theme(legend.position="bottom") +
  theme(legend.text = element_text(size=10)) +
  scale_x_continuous(name= TeX("Ability $(\\eta)$")) +
  scale_y_continuous(name= TeX("$P(u=c|\\eta)$"))

crf

This figure implies that in order for the first category on heart rate to be the most likely response category, the level of the latent trait for newborn health must be very low: below -4 or so. Since this category is labeled “absent”, that seems pretty severe and makes a lot of sense. Since the most likely response category at latent trait levels greater than -2 is category 2 (labeled a heart rate of 100 or more), this would imply that about 95% of the sample will have this response category observed, based on what we know about the normal distribution and recognizing that we assume the latent trait is normally distributed. We can check:

gmodels::CrossTable(apgar2i$u1)
## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  2263 
## 
##  
##           |         0 |         1 |         2 | 
##           |-----------|-----------|-----------|
##           |         1 |       111 |      2151 | 
##           |     0.000 |     0.049 |     0.951 | 
##           |-----------|-----------|-----------|
## 
## 
## 
##

And feel comfortable that the model is returning interpretable results and our plots agree with what we see in the data.