using Pumas
using DataFrames
using DataFramesMeta
using Random
using LinearAlgebra
using CairoMakie
using AlgebraOfGraphics
using SummaryTables
using PumasUtilities
using Pumas: OrderedLogit
Random.seed!(123)
Item Response Theory
Item response theory (IRT) describes how a set of discrete questionnaire or clinical-test responses (“items”) depends on an unobserved latent trait (such as disease severity). Each item has its own characteristic curve, driven by a difficulty (cut-off) and a discrimination, and the latent trait can vary over time and across subjects via a nonlinear mixed-effects (NLME) progression model. This tutorial walks through the different ways of defining and fitting an IRT model in Pumas.
1 Learning objectives
By the end of this tutorial you will be able to:
- Describe the bounded ordinal IRT model with item-specific discriminations.
- Simulate ordinal longitudinal responses from known parameters.
- Fit the IRT model by marginal likelihood maximization.
- Obtain asymptotic confidence intervals of the population parameters.
- Add weakly informative priors on the population parameters and fit the Bayesian model with maximum a-posteriori (MAP) estimation or Markov chain Monte Carlo (MCMC).
First, we load the necessary packages.
2 What is IRT?
IRT models describe the progression of multiple ordinal outcomes (e.g., questionnaire items, clinical tests) as a function of a latent trait (e.g., disease severity). The different outcomes are called items and the latent trait is often called ability. Each item has its own response curve, often called an item characteristic curve (ICC), which maps the latent trait to a probability of a certain response.
In a clinical setting, we might have a few tests that measure different levels of disease severity at different time points. Each test has its own characteristics (difficulty, discrimination) that determine how likely a patient with a certain disease severity is to respond positively to that test. IRT allows us to model these relationships and estimate both the latent trait (disease severity) and the item parameters (difficulty, discrimination) simultaneously from observed responses. For instance, for test A, a patient with a disease severity of 0.5 might have a 70% chance of responding positively (positive is bad), while the same patient might only have a 30% chance of responding positively to test B. As the disease severity increases, the probability of responding positively to both tests will increase, but at different rates depending on the item/test parameters.
Note
The disease severity (or ability) is a latent trait that we cannot directly observe, but we can infer it from the observed (discrete) responses to the items.
Tip
A mixture of discrete and continuous outcomes can also be modeled in a joint IRT framework, driven by a shared latent disease severity variable. For example, we might have a continuous biomarker that is measured alongside the discrete item responses, and we can model both types of outcomes together to get a more comprehensive understanding of the disease progression.
Warning
Simple IRT models assume that the latent trait (ability) is unidimensional and that all the items measure the same underlying trait. A violation of this assumption can lead to biased parameter estimates and incorrect inferences. If increasing the underlying disease severity leads to a higher score on some items but a lower score on others, then the unidimensionality assumption is likely violated and a more complex model (e.g., multidimensional IRT) may be needed. This can be checked during exploratory data analysis by looking at the correlation structure of the items’ responses. If the items are all positively correlated, this is a good sign that they are measuring the same underlying trait. On the other hand, if some items are negatively correlated with others, this may indicate that they are measuring different traits.
3 Bounded ordinal score models
With the conceptual picture in place, we now fix the notation that will carry us through the rest of the tutorial. Consider a study with \(N\) subjects responding to \(J\) items, where each item has a bounded ordinal response. The model below is the graded response model (GRM) of Samejima (1969), the polytomous extension of standard IRT, with item-specific discriminations. We use the following notation:
- \(i \in \{1, \dots, N\}\): subject index
- \(j \in \{1, \dots, J\}\): item index
- \(K_j\): number of response categories for item \(j\), so the response \(y_{i,j}\) takes values in \(\{1, 2, \dots, K_j\}\)
- \(k \in \{1, 2, \dots, K_j\}\): category index
- \(y_{i,j}\): observed ordinal response of subject \(i\) to item \(j\)
- \(f_i \in \mathbb{R}\): latent ability of subject \(i\)
- \(a_j > 0\): discrimination parameter of item \(j\) (controls the slope of the response curve)
- \(b_{j,k}\): cut-off (difficulty) for item \(j\) at category threshold \(k\), with \(b_{j,1} < b_{j,2} < \dots < b_{j,K_j - 1}\) ordered (\(K_j - 1\) cut-offs in total)
Important
Pumas’ OrderedLogit distribution assumes the response takes values in \(\{1, 2, \dots, K_j\}\), i.e., consecutive integers starting at 1. If your observed scores are ordinal but encoded differently (e.g., \(\{0, 1, 2, 3\}\), \(\{-2, -1, 0, 1, 2\}\), non-consecutive integers like \(\{0, 1, 2, 4\}\), fractional values such as \(\{0.0, 0.125, 0.25, \dots, 3.0\}\), or textual labels like "low" / "medium" / "high"), map them to \(1, 2, \dots, K_j\) before passing them to read_pumas or simobs. The mapping must preserve the ordering, and the same mapping should be used consistently for all subjects and time points so that a given integer always denotes the same category.
3.1 Category probabilities
The probability of each ordinal category is then defined as the difference between two cumulative probabilities. The cumulative probabilities are defined by the logistic (inverse logit) function as follows:
\[ \begin{aligned} P(y_{i,j} = k) & = P(y_{i,j} \leq k) - P(y_{i,j} \leq k - 1) \quad \forall k \in 1 \dots K_j \\ P(y_{i,j} \leq k) & = \begin{cases} 0 & k = 0 \\ 1 & k = K_j \\ \mathrm{logit}^{-1}\!\left(a_j \cdot (b_{j,k} - f_i)\right) & 1 \leq k \leq K_j - 1 \end{cases} \end{aligned} \]
Each cut-off \(b_{j,k}\) governs the cumulative probability \(P(y_{i,j} \leq k)\) of being in category \(k\) or below, with \(k = 1, \dots, K_j - 1\) indexing the \(K_j - 1\) cut-offs. The discrimination \(a_j > 0\) controls how sharply the response probability changes with ability: items with larger \(a_j\) have sharper transitions between categories (a small change in ability moves a lot of probability mass), while items with smaller \(a_j\) have flatter ICCs.
Note
The OrderedLogit distribution in Pumas. OrderedLogit(η) takes a vector \(\eta \in \mathbb{R}^{K-1}\) of (strictly increasing) values and defines a discrete distribution over \(\{1, 2, \dots, K\}\) via
\[ P(y \leq k) = \mathrm{logit}^{-1}(\eta_k), \qquad k = 1, \dots, K - 1, \]
with the boundary cases \(P(y \leq 0) = 0\) and \(P(y \leq K) = 1\). Each \(\eta_k\) is therefore the log-odds of being in category \(k\) or lower. Per-category probabilities are obtained by differencing:
\[ P(y = k) = P(y \leq k) - P(y \leq k - 1). \]
In the IRT model above, we set \(\eta_k = a_j (b_{j,k} - f_i)\), so larger ability \(f_i\) pushes \(\eta\) downward and shifts probability mass toward higher categories. The ordering \(b_{j,1} < \dots < b_{j,K_j - 1}\) implies \(\eta_1 < \dots < \eta_{K_j - 1}\), which is exactly the monotonicity that OrderedLogit requires.
For a more detailed treatment of ordinal regression with OrderedLogit in Pumas (including exploratory data analysis, covariate effects, and category-probability diagnostics), see the dedicated Ordinal Regression tutorial, which the IRT model here builds upon.
3.2 Example: visualizing ICCs
To build intuition, the figure below shows item characteristic curves (ICCs) for a single ordinal item with \(K = 4\) categories (\(y \in \{1, 2, 3, 4\}\) and \(K - 1 = 3\) cut-offs \(b_{j,1}, b_{j,2}, b_{j,3}\)). The left panel shows the cumulative probabilities \(P(y_{i,j} \leq k)\) as a function of ability \(f_i\), and the right panel shows the per-category probabilities \(P(y_{i,j} = k)\).
logistic(x) = 1 / (1 + exp(-x))
fi = -4:0.05:4
# Single item: 4 categories, 3 cut-offs
a_ex = 1.5
b_ex = [-1.0, 0.0, 1.2]
# Cumulative probabilities P(y ≤ k) for k = 1, 2, 3
cum_probs = [logistic.(a_ex .* (b .- fi)) for b in b_ex]
# Category probabilities P(y = k) for k = 1, 2, 3, 4
p_y1 = cum_probs[1]
p_y2 = cum_probs[2] .- cum_probs[1]
p_y3 = cum_probs[3] .- cum_probs[2]
p_y4 = 1 .- cum_probs[3]
icc_cum_df = DataFrame(
ability = repeat(collect(fi), 3),
probability = vcat(cum_probs...),
threshold = repeat(["P(y ≤ 1)", "P(y ≤ 2)", "P(y ≤ 3)"]; inner = length(fi)),
)
icc_cat_df = DataFrame(
ability = repeat(collect(fi), 4),
probability = vcat(p_y1, p_y2, p_y3, p_y4),
category = repeat(["y = 1", "y = 2", "y = 3", "y = 4"]; inner = length(fi)),
)
fig = Figure(; size = (700, 300))
ag1 = draw!(
fig[1, 1],
data(icc_cum_df) *
mapping(:ability => "ability (f)", :probability, color = :threshold => nonnumeric) *
visual(Lines; linewidth = 3);
axis = (; title = "Cumulative ICCs"),
)
legend!(fig[1, 2], ag1)
ag2 = draw!(
fig[1, 3],
data(icc_cat_df) *
mapping(:ability => "ability (f)", :probability, color = :category => nonnumeric) *
visual(Lines; linewidth = 3);
axis = (; title = "Per-category ICCs"),
)
legend!(fig[1, 4], ag2)
fig
The cumulative curves are ordered logistic functions of ability, each centred at one cut-off \(b_{j,k}\) and sharing the slope \(a_j\). The per-category probabilities sum to 1 at every value of \(f_i\); the middle categories appear as bumps between adjacent cut-offs, while the extreme categories (\(y = 1\) and \(y = K_j\)) are monotone in ability: \(P(y = 1)\) decreasing and \(P(y = K_j)\) increasing as \(f_i\) grows.
4 A Pumas ordinal IRT model with a linear latent ability
Having grounded the cross-sectional picture in a plot, we now turn it into a longitudinal Pumas model. Concretely, we build a graded response model with \(J = 2\) items (score1, score2), each with \(K = 4\) response categories, and each subject measured at multiple time points. In this first example, the latent ability is taken to be linear in time:
\[ \begin{aligned} f_i(t) & = f_{0,i} + s_i \cdot t \\ f_{0,i} & \sim \mathrm{Normal}\!\left(\mu_f,\; \omega_f^2\right) \\ s_i & \sim \mathrm{Normal}\!\left(\mu_s,\; \omega_s^2\right) \end{aligned} \]
Both the intercept \(f_{0,i}\) and the slope \(s_i\) of this linear trajectory carry fixed and random effects: the population means \(\mu_f, \mu_s\) describe the typical disease trajectory, while the per-subject draws \(f_{0,i}, s_i\) around those means capture inter-individual variability. This makes the model a nonlinear mixed effects (NLME) model. Concretely, the population-level parameters are:
- \(\mu_f, \mu_s\): fixed-effect (population) mean intercept and progression slope of ability
- \(\omega_f^2, \omega_s^2\): variances of the random effects on the intercept and slope (so \(\omega_f, \omega_s\) are SDs)
- \(f_{0,i}, s_i\): subject-specific random effects (the actual baseline ability and slope for subject \(i\))
The item-level parameters per item \(j \in \{1, 2\}\) are:
- \(b_j \in \mathbb{R}^{K_j - 1}\): ordered per-category cut-offs of item \(j\), with \(b_{j,1} < b_{j,2} < \dots < b_{j,K_j - 1}\). The length of \(b_j\) is the number of cut-offs, one fewer than the number of categories. Each item has its own cut-off parameter, accommodating items with different numbers of categories.
- \(a_j > 0\): discrimination of item \(j\) (positive scalar).
The OrderedLogit threshold input for item \(j\) at time \(t\) is \(\eta_{j,k} = a_j \cdot (b_{j,k} - f_i(t))\) for each cut-off \(k = 1, \dots, K_j - 1\), matching the bounded ordinal model defined earlier: \(\mathrm{logit}^{-1}(\eta_{j,k}) = P(y_{i,j} \leq k)\).
We stick with the linear form throughout the rest of this tutorial because it keeps the focus on the IRT-specific parts of the model. More generally, any disease progression model can play the role of \(f_i(t)\), such as exponential, EMAX, indirect-response, or mechanistic ODE-based forms. Switching between progressions is a localized change to @pre (and optionally @dynamics) and its associated parameters; the rest of the model (item cut-offs, discriminations, OrderedLogit observation) is unaffected. We come back to this in the Extensions section.
4.1 Enforcing ordering via an unconstrained parameterization
The ordered cut-off vector \(b_j\) lives in a constrained region of \(\mathbb{R}^{K_j - 1}\), namely the strictly increasing sequences. To let the optimizer work in a simple, unconstrained space, we introduce an unconstrained surrogate \(\theta^b_j \in \mathbb{R}^{K_j - 1}\) (called θb1 and θb2 in the Pumas @param block) and obtain the ordered cut-offs \(b_j\) from \(\theta^b_j\) via a cumulative log-exp transform:
\[ \begin{aligned} b_{j,1} & = \theta^b_{j,1} & & \text{(first cut-off is unconstrained on the real line)} \\ \theta^b_{j,k} & = \log\!\left(b_{j,k} - b_{j,k-1}\right), \quad k > 1 & & \text{(log of the gap between consecutive cut-offs)} \\ b_{j,k} & = b_{j,k-1} + e^{\theta^b_{j,k}}, \quad k > 1 & & \text{(reconstruct ordered cut-offs)} \end{aligned} \]
So when item \(j\)’s first cut-off is estimated:
- \(\theta^b_{j,1}\) is the location of the lowest cut-off on the ability scale, free to take any real value.
- \(\theta^b_{j,k}\) for \(k > 1\) is the logarithm of the gap between consecutive cut-offs. Because \(e^{\theta^b_{j,k}} > 0\), the reconstructed \(b_{j,k}\) sequence is strictly increasing for any value of \(\theta^b_j\), so the ordering constraint is satisfied automatically. The \(\theta\) prefix distinguishes them from the ordered cut-offs \(b_j\).
When item \(j\)’s first cut-off is instead fixed (which we do for item 1 below, to anchor the latent scale), the lowest cut-off drops out of \(\theta^b_j\) entirely: \(\theta^b_j\) then consists only of log-gaps, has length \(K_j - 2\), and the reconstruction prepends the fixed value before applying the same cumulative log-exp transform. So the meaning of \(\theta^b_{j,1}\) differs between the anchor item (a log-gap) and the rest (a location).
The reparameterization is implemented as a pair of plain Julia functions used both inside the model and when post-processing the fit:
# Map an unconstrained vector θ ∈ R^m to an ordered vector b ∈ R^m:
# b[1] = θ[1], b[k] = b[k-1] + exp(θ[k]) for k > 1.
function to_ordered(θ)
b = similar(θ)
b[1] = θ[1]
for k = 2:length(θ)
b[k] = b[k-1] + exp(θ[k])
end
return b
end
# Inverse: map an ordered vector b ∈ R^m to its unconstrained θ ∈ R^m:
# θ[1] = b[1], θ[k] = log(b[k] - b[k-1]) for k > 1.
function to_raw(b)
θ = similar(b)
θ[1] = b[1]
for k = 2:length(b)
θ[k] = log(b[k] - b[k-1])
end
return θ
end4.2 Pinning a reference cut-off for location identifiability
The model is invariant under a joint shift of the latent ability \(f\) and all cut-offs \(b_{j,k}\) by the same constant, so without an external constraint, the location of \(f\) (and hence of every \(b_{j,k}\)) is unidentified. Only the relative positions of the cut-offs across items carry information; the common origin is a convention. To make the parameters recoverable, we therefore pick any one cut-off and pin it to a fixed value (a “reference cut-off”), and let everything else vary freely. The choice of which cut-off is arbitrary (it does not have to be the lowest, the highest, or come from any particular item). For concreteness, we use
\[ b_{1,1} \;=\; 0 \]
throughout this tutorial. In the Pumas model, this is hard-coded directly inside @pre rather than declared in @param, so \(b_{1,1}\) is not a fitted parameter at all. Concretely, \(\theta^b_1\) now holds only the \(K_1 - 2\) log-gaps of item 1, and the ordered cut-offs are reconstructed by prepending \(0\) before applying the cumulative transform.
4.3 Pinning a reference discrimination for scale identifiability
Location is not the only invariance. The model is also invariant under a joint rescaling: for any \(c > 0\), mapping \(f_i \to c f_i\), \(b_{j,k} \to c b_{j,k}\), and \(a_j \to a_j / c\) leaves \(\eta_{j,k} = a_j(b_{j,k} - f_i)\) unchanged. Fixing \(b_{1,1} = 0\) does not resolve this; it only kills the shift. We therefore additionally pin the discrimination of item 1,
\[ a_1 \;=\; 1, \]
also hard-coded in @pre. Only the remaining discriminations \(a_2, \dots, a_J\) are estimated, and they are interpreted relative to item 1: \(a_j = 1.5\) means item \(j\)’s response curve is 1.5× as steep as item 1’s, not 1.5 on some absolute scale.
4.4 The Pumas model
The Pumas model for the IRT model with a linear latent ability trajectory is defined as follows.
n_items = 2
# Per-item number of categories; can differ across items.
# Each item j has n_cats[j] categories and n_cats[j] - 1 cut-offs.
n_cats = [4, 4]
m = n_cats .- 1
irt_ordinal_model = @model begin
@param begin
# Ability fixed effects (population intercept and slope)
μf ∈ RealDomain()
μs ∈ RealDomain()
# Ability random-effect SDs (subject-level variability in intercept and slope)
ωf ∈ RealDomain(; lower = 0.05)
ωs ∈ RealDomain(; lower = 1e-3)
# Per-item log-gap vectors. Item 1's first cut-off is hard-coded to 0
# in @pre as the location reference (not estimated), so θb1 contributes
# only m[1] - 1 free entries (the log-gaps). Item 2 contributes m[2]
# entries (first cut-off + log-gaps).
θb1 ∈ VectorDomain(m[1] - 1)
θb2 ∈ VectorDomain(m[2])
# Discriminations for items 2..n_items, on the positive scale. Item 1's
# discrimination is hard-coded to 1 in @pre as the scale reference, so
# a_free has length n_items - 1.
a_free ∈ VectorDomain(n_items - 1; lower = 0.0)
end
@random begin
# Subject-specific intercept and slope drawn directly from their
# population distributions, parameterized by mean and SD
f0 ~ Normal(μf, ωf)
slope ~ Normal(μs, ωs)
end
@pre begin
# Latent ability trajectory for the subject
f = f0 + slope * t
# Reconstruct ordered cut-offs b_j. For item 1, prepend the reference
# cut-off (hard-coded to 0) before the log-gaps θb1.
b1_ordered = to_ordered(vcat(0.0, θb1))
b2_ordered = to_ordered(θb2)
# Full discrimination vector with a[1] = 1 fixed as the scale reference.
a = vcat(1.0, a_free)
# OrderedLogit thresholds: η_{j,k} = a_j (b_{j,k} - f)
lle1 = a[1] .* (b1_ordered .- f)
lle2 = a[2] .* (b2_ordered .- f)
end
@derived begin
score1 ~ @. OrderedLogit(lle1)
score2 ~ @. OrderedLogit(lle2)
end
endPumasModel
Parameters: μf, μs, ωf, ωs, θb1, θb2, a_free
Random effects: f0, slope
Covariates:
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: score1, score2
Observed: score1, score25 Simulate data
With the model defined, we generate a synthetic dataset that we can later fit it on, so that we can compare the recovered estimates against parameter values we actually know. We simulate ordinal responses for 300 subjects measured at 5 time points (e.g., weeks 0, 4, 8, 12, 24) from fixed ordinal parameter values using Pumas’ simobs. First we build a template population (one Subject per id), then call simobs at the true parameter values.
n_subjects = 300
obs_times = [0.0, 4.0, 8.0, 12.0, 24.0]
# Target ordered cut-offs per item (length m[j] = n_cats[j] - 1).
# Item 1 anchors the latent scale: its first cut-off is fixed at 0.
b1_ordered_true = [0.0, 1.5, 2.5]
b2_ordered_true = [-1.0, 0.3, 1.2]
# Drop the (fixed-at-0) first cut-off of item 1; keep only the log-gaps,
# matching the @param decomposition where θb1 holds the log-gaps alone.
# a_free holds only the free discriminations (items 2..n_items); item 1's
# discrimination is fixed to 1 inside @pre as the scale reference.
true_params = (
μf = -0.3,
ωf = 0.8,
μs = 0.04,
ωs = 0.02,
θb1 = to_raw(b1_ordered_true)[2:end],
θb2 = to_raw(b2_ordered_true),
a_free = [0.75],
)
template_pop = map(1:n_subjects) do i
Subject(; id = i, time = obs_times)
end
sim_obs = simobs(irt_ordinal_model, template_pop, true_params)
irt_df = DataFrame(sim_obs)
first(irt_df, 10)10×13 DataFrame
| Row | id | time | score1 | score2 | evid | f0 | slope | f | b1_ordered | b2_ordered | a | lle1 | lle2 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| String | Float64 | Int64? | Int64? | Int64? | Float64? | Float64? | Float64? | Array…? | Array…? | Array…? | Array…? | Array…? | |
| 1 | 1 | 0.0 | 2 | 4 | 0 | 0.482736 | 0.00734241 | 0.482736 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-0.482736, 1.01726, 2.01726] | [-1.11205, -0.137052, 0.537948] |
| 2 | 1 | 4.0 | 2 | 1 | 0 | 0.482736 | 0.00734241 | 0.512106 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-0.512106, 0.987894, 1.98789] | [-1.13408, -0.159079, 0.515921] |
| 3 | 1 | 8.0 | 2 | 2 | 0 | 0.482736 | 0.00734241 | 0.541475 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-0.541475, 0.958525, 1.95852] | [-1.15611, -0.181106, 0.493894] |
| 4 | 1 | 12.0 | 1 | 2 | 0 | 0.482736 | 0.00734241 | 0.570845 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-0.570845, 0.929155, 1.92916] | [-1.17813, -0.203134, 0.471866] |
| 5 | 1 | 24.0 | 2 | 4 | 0 | 0.482736 | 0.00734241 | 0.658954 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-0.658954, 0.841046, 1.84105] | [-1.24422, -0.269215, 0.405785] |
| 6 | 2 | 0.0 | 1 | 4 | 0 | 1.56184 | 0.0543695 | 1.56184 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-1.56184, -0.0618438, 0.938156] | [-1.92138, -0.946383, -0.271383] |
| 7 | 2 | 4.0 | 3 | 4 | 0 | 1.56184 | 0.0543695 | 1.77932 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-1.77932, -0.279322, 0.720678] | [-2.08449, -1.10949, -0.434491] |
| 8 | 2 | 8.0 | 1 | 4 | 0 | 1.56184 | 0.0543695 | 1.9968 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-1.9968, -0.4968, 0.5032] | [-2.2476, -1.2726, -0.5976] |
| 9 | 2 | 12.0 | 2 | 4 | 0 | 1.56184 | 0.0543695 | 2.21428 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-2.21428, -0.714278, 0.285722] | [-2.41071, -1.43571, -0.760708] |
| 10 | 2 | 24.0 | 4 | 3 | 0 | 1.56184 | 0.0543695 | 2.86671 | [0.0, 1.5, 2.5] | [-1.0, 0.3, 1.2] | [1.0, 0.75] | [-2.86671, -1.36671, -0.366712] | [-2.90003, -1.92503, -1.25003] |
5.1 Checking category coverage
Before fitting, it is worth checking that every response category for every score is observed at every time point; otherwise some cut-offs will be poorly identified by the likelihood:
@chain irt_df begin
groupby([:time, :score1])
combine(nrow => :n)
simple_table
end| Time | score1 | n |
| 0 | 1 | 197 |
| 0 | 2 | 62 |
| 0 | 3 | 26 |
| 0 | 4 | 15 |
| 4 | 1 | 172 |
| 4 | 2 | 87 |
| 4 | 3 | 26 |
| 4 | 4 | 15 |
| 8 | 1 | 159 |
| 8 | 2 | 88 |
| 8 | 3 | 31 |
| 8 | 4 | 22 |
| 12 | 1 | 140 |
| 12 | 2 | 100 |
| 12 | 3 | 39 |
| 12 | 4 | 21 |
| 24 | 1 | 142 |
| 24 | 2 | 69 |
| 24 | 3 | 43 |
| 24 | 4 | 46 |
@chain irt_df begin
groupby([:time, :score2])
combine(nrow => :n)
simple_table
end| Time | score2 | n |
| 0 | 1 | 111 |
| 0 | 2 | 75 |
| 0 | 3 | 43 |
| 0 | 4 | 71 |
| 4 | 1 | 113 |
| 4 | 2 | 68 |
| 4 | 3 | 30 |
| 4 | 4 | 89 |
| 8 | 1 | 95 |
| 8 | 2 | 65 |
| 8 | 3 | 48 |
| 8 | 4 | 92 |
| 12 | 1 | 101 |
| 12 | 2 | 61 |
| 12 | 3 | 47 |
| 12 | 4 | 91 |
| 24 | 1 | 71 |
| 24 | 2 | 54 |
| 24 | 3 | 49 |
| 24 | 4 | 126 |
All four categories appear with a reasonable count at every time point for both items, so we proceed.
5.2 Reading into a Pumas population
We convert the simulated DataFrame into a Pumas population, ready for fitting:
irt_pop = read_pumas(
irt_df;
id = :id,
time = :time,
observations = [:score1, :score2],
event_data = false,
)Population
Subjects: 300
Observations: score1, score26 Fit the IRT model with LaplaceI
With the simulated population defined, we now estimate the population parameters using the LaplaceI() algorithm. We deliberately initialize the optimizer away from the true values so that recovery is not trivial:
# Initial unconstrained θb vectors, split out for clarity:
# - θb1_init holds only the log-gaps of item 1 (length m[1] - 1); the first
# cut-off is hard-coded to 0 inside @pre and is not estimated.
# - θb2_init holds the first cut-off and log-gaps of item 2 (length m[2]).
θb1_init = fill(-0.5, m[1] - 1)
θb2_init = vcat(0.0, fill(-0.5, m[2] - 1))
init_irt = (
μf = 0.5,
ωf = 0.4,
μs = 0.2,
ωs = 0.1,
θb1 = θb1_init,
θb2 = θb2_init,
a_free = fill(0.5, n_items - 1),
)(μf = 0.5, ωf = 0.4, μs = 0.2, ωs = 0.1, θb1 = [-0.5, -0.5], θb2 = [0.0, -0.5, -0.5], a_free = [0.5])fpm_irt_laplace = fit(irt_ordinal_model, irt_pop, init_irt, LaplaceI())[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 4.800428e+03 4.779925e+03 * time: 0.03737187385559082 1 4.331239e+03 2.074564e+03 * time: 6.526118040084839 2 4.164350e+03 3.266075e+02 * time: 7.230453968048096 3 4.074510e+03 4.890505e+02 * time: 7.866087913513184 4 3.855667e+03 1.975355e+03 * time: 8.47200894355774 5 3.758155e+03 1.080697e+03 * time: 9.148177862167358 6 3.709927e+03 4.542381e+02 * time: 9.740301847457886 7 3.695732e+03 3.184177e+02 * time: 10.290653944015503 8 3.660219e+03 1.091091e+02 * time: 10.84486985206604 9 3.655934e+03 9.540397e+01 * time: 11.415019035339355 10 3.630042e+03 1.664781e+02 * time: 11.971274852752686 11 3.618930e+03 7.081763e+01 * time: 12.509055852890015 12 3.614430e+03 5.305270e+01 * time: 13.06850290298462 13 3.613672e+03 3.028912e+01 * time: 13.625885009765625 14 3.613007e+03 2.132361e+01 * time: 14.171612024307251 15 3.611778e+03 5.387535e+01 * time: 14.762471914291382 16 3.610270e+03 1.577074e+01 * time: 15.307481050491333 17 3.609388e+03 2.147710e+01 * time: 15.899247884750366 18 3.609051e+03 1.736887e+01 * time: 16.47326397895813 19 3.608733e+03 1.774374e+01 * time: 17.050946950912476 20 3.608112e+03 1.990767e+01 * time: 17.583266019821167 21 3.607103e+03 2.192472e+01 * time: 18.138062953948975 22 3.606085e+03 2.264348e+01 * time: 18.71022605895996 23 3.605450e+03 1.753239e+01 * time: 19.315476894378662 24 3.605051e+03 1.535280e+01 * time: 19.84327507019043 25 3.604787e+03 1.588642e+01 * time: 20.384187936782837 26 3.604144e+03 1.741700e+01 * time: 20.955333948135376 27 3.603174e+03 1.633688e+01 * time: 21.545450925827026 28 3.602358e+03 1.079944e+01 * time: 22.093586921691895 29 3.601993e+03 1.052361e+01 * time: 22.646564960479736 30 3.601874e+03 1.168567e+01 * time: 23.263288021087646 31 3.601821e+03 1.131774e+01 * time: 23.913433074951172 32 3.601703e+03 9.676022e+00 * time: 24.483294010162354 33 3.601507e+03 8.356853e+00 * time: 25.05222988128662 34 3.601288e+03 1.015789e+01 * time: 25.615659952163696 35 3.601165e+03 8.779162e+00 * time: 26.238352060317993 36 3.601111e+03 6.432333e+00 * time: 26.753882884979248 37 3.601075e+03 4.746521e+00 * time: 27.305309057235718 38 3.601013e+03 6.194610e+00 * time: 27.848477840423584 39 3.600894e+03 7.217212e+00 * time: 28.38909602165222 40 3.600739e+03 8.244509e+00 * time: 28.92282485961914 41 3.600612e+03 6.055542e+00 * time: 29.46276307106018 42 3.600561e+03 2.519780e+00 * time: 30.02649688720703 43 3.600554e+03 2.595838e+00 * time: 30.616481065750122 44 3.600551e+03 2.580745e+00 * time: 31.149869918823242 45 3.600540e+03 2.577949e+00 * time: 31.768063068389893 46 3.600515e+03 4.724372e+00 * time: 32.350982904434204 47 3.600455e+03 7.989893e+00 * time: 32.91207003593445 48 3.600354e+03 1.098940e+01 * time: 33.45792889595032 49 3.600257e+03 1.003705e+01 * time: 34.005555868148804 50 3.600196e+03 3.905297e+00 * time: 34.79412007331848 51 3.600186e+03 2.647766e-01 * time: 35.80467104911804 52 3.600185e+03 8.021899e-02 * time: 36.54373598098755 53 3.600185e+03 7.332845e-03 * time: 37.2733850479126 54 3.600185e+03 6.224305e-03 * time: 37.89951491355896 55 3.600185e+03 7.417728e-04 * time: 38.4845929145813
FittedPumasModel
Dynamical system type: No dynamical model
Number of subjects: 300
Observation records: Active Missing
score1: 1500 0
score2: 1500 0
Total: 3000 0
Number of parameters: Constant Optimized
0 10
Likelihood approximation: LaplaceI
Likelihood optimizer: BFGS
Termination Reason: GradientNorm
Log-likelihood value: -3600.1853
--------------------
Estimate
--------------------
μf -0.60264
μs 0.04482
ωf 0.8417
ωs 0.038788
θb1₁ 0.44853
θb1₂ 0.14895
θb2₁ -1.4797
θb2₂ 0.47622
θb2₃ 0.11072
a_free₁ 0.60231
--------------------Recall the true parameters used to simulate the data, to compare against the estimates above:
true_params(μf = -0.3, ωf = 0.8, μs = 0.04, ωs = 0.02, θb1 = [0.4054651081081644, 0.0], θb2 = [-1.0, 0.26236426446749106, -0.1053605156578264], a_free = [0.75])Despite the deliberately off starting point, the optimizer converges to estimates closer (but not identical) to the parameters used to simulate the data. Next, we check whether the confidence intervals of the estimates contain the true values.
7 Confidence intervals
Next, we want to get confidence intervals (CIs) for the population parameters and check that they contain the true values we used to simulate the data. Asymptotic confidence intervals for the population parameters can be obtained using infer on the LaplaceI fitted Pumas model.
infer_irt = infer(fpm_irt_laplace)
ci_df = coeftable(infer_irt)
simple_table(ci_df)[ Info: Calculating: variance-covariance matrix. [ Info: Done.
| parameter | description | constant | estimate | se | relative_se | ci_lower | ci_upper |
| μf | μf | false | -0.603 | 0.0957 | 0.159 | -0.79 | -0.415 |
| μs | μs | false | 0.0448 | 0.00601 | 0.134 | 0.033 | 0.0566 |
| ωf | ωf | false | 0.842 | 0.0983 | 0.117 | 0.649 | 1.03 |
| ωs | ωs | false | 0.0388 | 0.00864 | 0.223 | 0.0218 | 0.0557 |
| θb1₁ | θb1₁ | false | 0.449 | 0.0477 | 0.106 | 0.355 | 0.542 |
| θb1₂ | θb1₂ | false | 0.149 | 0.0818 | 0.549 | -0.0114 | 0.309 |
| θb2₁ | θb2₁ | false | -1.48 | 0.223 | 0.15 | -1.92 | -1.04 |
| θb2₂ | θb2₂ | false | 0.476 | 0.144 | 0.302 | 0.194 | 0.758 |
| θb2₃ | θb2₃ | false | 0.111 | 0.15 | 1.35 | -0.183 | 0.404 |
| a_free₁ | a_free₁ | false | 0.602 | 0.0844 | 0.14 | 0.437 | 0.768 |
We then append true_value and init_value columns by flattening true_params and init_irt in the same row order that coeftable uses. coeftable lists the population scalars as μf, μs, ωf, ωs (all means first, then all SDs), followed by the vector parameters θb1, θb2, a_free:
ci_df.true_value = vcat(
[true_params.μf, true_params.μs, true_params.ωf, true_params.ωs],
true_params.θb1,
true_params.θb2,
true_params.a_free,
)
ci_df.init_value = vcat(
[init_irt.μf, init_irt.μs, init_irt.ωf, init_irt.ωs],
init_irt.θb1,
init_irt.θb2,
init_irt.a_free,
)
simple_table(ci_df)| parameter | description | constant | estimate | se | relative_se | ci_lower | ci_upper | true_value | init_value |
| μf | μf | false | -0.603 | 0.0957 | 0.159 | -0.79 | -0.415 | -0.3 | 0.5 |
| μs | μs | false | 0.0448 | 0.00601 | 0.134 | 0.033 | 0.0566 | 0.04 | 0.2 |
| ωf | ωf | false | 0.842 | 0.0983 | 0.117 | 0.649 | 1.03 | 0.8 | 0.4 |
| ωs | ωs | false | 0.0388 | 0.00864 | 0.223 | 0.0218 | 0.0557 | 0.02 | 0.1 |
| θb1₁ | θb1₁ | false | 0.449 | 0.0477 | 0.106 | 0.355 | 0.542 | 0.405 | -0.5 |
| θb1₂ | θb1₂ | false | 0.149 | 0.0818 | 0.549 | -0.0114 | 0.309 | 0 | -0.5 |
| θb2₁ | θb2₁ | false | -1.48 | 0.223 | 0.15 | -1.92 | -1.04 | -1 | 0 |
| θb2₂ | θb2₂ | false | 0.476 | 0.144 | 0.302 | 0.194 | 0.758 | 0.262 | -0.5 |
| θb2₃ | θb2₃ | false | 0.111 | 0.15 | 1.35 | -0.183 | 0.404 | -0.105 | -0.5 |
| a_free₁ | a_free₁ | false | 0.602 | 0.0844 | 0.14 | 0.437 | 0.768 | 0.75 | 0.5 |
Each row now shows the LaplaceI point estimate, its confidence interval, the true value used to simulate the data, and the initial value handed to the optimizer, making it easy to check both that every true parameter falls inside its CI and that the fit moved away from its starting point.
8 Predicted item characteristic curves
We opened the tutorial by plotting ICCs from hypothetical parameters; now we plot the estimated ICCs and overlay the true curves. We reconstruct each item’s cut-offs and discrimination from the fit (via to_ordered, prepending the pinned 0 for item 1) and reuse the per-category probability construction from the introductory example.
# Reconstruct estimated item parameters from the LaplaceI fit
est = coef(fpm_irt_laplace)
b1_hat = to_ordered(vcat(0.0, est.θb1))
b2_hat = to_ordered(est.θb2)
a_hat = vcat(1.0, est.a_free)
# Per-category probabilities P(y = k) over an ability grid for one item,
# given its discrimination `a` and ordered cut-offs `b`.
function category_probs(a, b, grid)
cum = [logistic.(a .* (bk .- grid)) for bk in b] # P(y ≤ k), k = 1 … K-1
K = length(b) + 1
return [k == 1 ? cum[1] : k == K ? 1 .- cum[end] : cum[k] .- cum[k-1] for k = 1:K]
end
# Tidy long-format ICC table for one item and one parameter source.
function icc_long(a, b, grid, item, source)
probs = category_probs(a, b, grid)
K = length(probs)
DataFrame(
ability = repeat(collect(grid), K),
probability = reduce(vcat, probs),
category = repeat(["y = $k" for k = 1:K]; inner = length(grid)),
item = "item $item",
source = source,
)
end
icc_compare = reduce(
vcat,
[
icc_long(a_hat[1], b1_hat, fi, 1, "estimated"),
icc_long(a_hat[2], b2_hat, fi, 2, "estimated"),
icc_long(1.0, b1_ordered_true, fi, 1, "true"),
icc_long(true_params.a_free[1], b2_ordered_true, fi, 2, "true"),
],
)
plt =
data(icc_compare) *
mapping(
:ability => "ability (f)",
:probability => "P(y = k)",
color = :category => nonnumeric,
linestyle = :source => nonnumeric => "source",
col = :item => nonnumeric,
) *
visual(Lines; linewidth = 2)
draw(plt; figure = (; size = (700, 300)))
The estimated curves (solid) sit close to the true curves (dashed) for both items, confirming visually what the confidence-interval table reported. Item 2’s flatter transitions reflect its smaller estimated discrimination.
9 Model diagnostics and qualification
With real data there is no ground truth to compare against, so we rely on diagnostics that ask whether the fitted model is consistent with the observed responses.
9.1 Information criteria
aic and bic are most useful for comparing competing models (e.g. a linear vs. an EMAX latent trajectory) fitted to the same data:
(; AIC = aic(fpm_irt_laplace), BIC = bic(fpm_irt_laplace))(AIC = 7220.3706800337, BIC = 7280.434355710203)9.2 Visual predictive check
The key qualification tool is the visual predictive check (VPC), which overlays the observed category proportions over time on the distribution simulated from the fitted model. Pumas builds a categorical VPC with one panel per response category, so we plot one item at a time. The panels default to first-appearance order, so we pass levels to force them to read \(1 \to K\). Pumas names the per-category columns <observation>_<value>, which we reconstruct from the sorted category values:
# Build sorted level names matching Pumas' `<observation>_<value>` columns.
pop_df = DataFrame(irt_pop)
ordered_levels(score) =
[Symbol("$(score)_$(v)") for v in sort(unique(skipmissing(pop_df[!, score])))]ordered_levels (generic function with 1 method)irt_vpc_1 = vpc(fpm_irt_laplace; observations = [:score1])
vpc_plot(irt_vpc_1; levels = ordered_levels(:score1), figure = (; size = (600, 1200)))[ Info: Discrete VPC [ Info: Automatically choosing bandwidth for smooth discrete VPC. [ Info: Automatic bandwidth choice complete with bandwidth: 24.0.
irt_vpc_2 = vpc(fpm_irt_laplace; observations = [:score2])
vpc_plot(irt_vpc_2; levels = ordered_levels(:score2), figure = (; size = (600, 1200)))[ Info: Discrete VPC [ Info: Automatically choosing bandwidth for smooth discrete VPC. [ Info: Automatic bandwidth choice complete with bandwidth: 24.0.
Each panel shows one response category’s observed proportion over time (for one item) against the simulated prediction band; bands that envelop the observed proportions across categories and items indicate an adequate fit.
10 Priors as regularization when data is sparse
With 300 subjects, every category was observed enough times that pure maximum likelihood worked fine. But what happens when we don’t have that much data?
The fit we just ran uses marginal maximum likelihood (LaplaceI()): each parameter is declared in @param with ∈ RealDomain(...) or ∈ VectorDomain(...), which specifies only its support, with no probabilistic information. With the \(N = 300\) subjects and 5 observations per item per subject we just used, every category of every item is observed enough times that the likelihood alone identifies all parameters well, and MLE works fine.
The story changes when some categories are barely observed: for example, a clinical scale where the worst category is reached by only a handful of subjects, or an item where one threshold is rarely crossed. In that regime, the corresponding cut-off has very little likelihood information, leading to numerical instability and non-estimability of some parameters.
10.1 A prior-aware Pumas model
The standard remedy is to regularize via priors, i.e., to fit by maximum a posteriori (MAP) instead of MLE. In Pumas, this means redeclaring each parameter with the ~ syntax instead of ∈, which both fixes its support and attaches a prior distribution:
irt_ordinal_model_map = @model begin
@param begin
# Ability fixed effects: weakly informative Gaussians centred at 0.
μf ~ Normal(0, 2)
μs ~ Normal(0, 2)
# Random-effect SDs: half-Cauchy on the positive line. Heavy tails
# allow large values; the small lower bound keeps the optimizer
# from collapsing the SD to exactly zero.
ωf ~ truncated(Cauchy(0, 1); lower = 0.05)
ωs ~ truncated(Cauchy(0, 1); lower = 1e-3)
# Item 1: θb1 contains only log-gaps (first cut-off is fixed to 0).
# N(0, 1) on each log-gap implies a prior 95% interval of roughly
# (e^{-2}, e^{2}) ≈ (0.14, 7.4) for the gap between consecutive
# cut-offs. This is broad relative to typical category spacings but not
# absurd.
θb1 ~ MvNormal(zeros(m[1] - 1), I(m[1] - 1))
# Item 2: θb2[1] is the lowest cut-off (a location on the latent
# scale) and gets a wider N(0, 2) prior; θb2[2:end] are log-gaps
# and get the same N(0, 1) prior as for item 1.
θb2 ~ MvNormal(zeros(m[2]), Diagonal(vcat(4.0, ones(m[2] - 1))))
# Free item discriminations (items 2..n_items): N(0, 2) on log(a_j).
# Item 1's discrimination is fixed to 1 as the scale reference.
a_free ~ MvLogNormal(zeros(n_items - 1), 4 * I(n_items - 1))
end
@random begin
f0 ~ Normal(μf, ωf)
slope ~ Normal(μs, ωs)
end
@pre begin
f = f0 + slope * t
b1_ordered = to_ordered(vcat(0.0, θb1))
b2_ordered = to_ordered(θb2)
a = vcat(1.0, a_free)
lle1 = a[1] .* (b1_ordered .- f)
lle2 = a[2] .* (b2_ordered .- f)
end
@derived begin
score1 ~ @. OrderedLogit(lle1)
score2 ~ @. OrderedLogit(lle2)
end
endPumasModel
Parameters: μf, μs, ωf, ωs, θb1, θb2, a_free
Random effects: f0, slope
Covariates:
Dynamical system variables:
Dynamical system type: No dynamical model
Derived: score1, score2
Observed: score1, score2Only the @param block has changed; everything downstream is identical. The prior choices above are deliberately weakly informative, and they treat locations on the latent scale and log-gaps between cut-offs differently. A wide prior on the log scale corresponds to an absurdly wide prior on the gap itself, which is rarely what we want:
- \(\mathrm{Normal}(0, 2)\) on the ability fixed effects (\(\mu_f, \mu_s\)) and the lowest cut-off of item 2 (\(\theta^b_{2,1}\), which lives on the ability scale). Wide enough not to dominate informative data, but tight enough to keep the latent scale anchored.
- \(\mathrm{Normal}(0, 1)\) on every log-gap component (\(\theta^b_{1,:}\) and \(\theta^b_{2,2:}\)). This places the prior 95% interval for each consecutive-cut-off gap at \((e^{-2}, e^{2}) \approx (0.14, 7.4)\). Bumping this up to \(\mathrm{Normal}(0, 2)\) (the prior we use on locations) would imply gaps anywhere in \((e^{-4}, e^{4}) \approx (0.02, 55)\), which is far too diffuse and gives the optimizer / sampler very little to hold on to in sparse-data regimes.
- \(\mathrm{LogNormal}(0, 2)\) on the free item discriminations \(a_{2:J}\) (item 1’s discrimination is fixed to 1 as the scale reference, so it has no prior). Equivalent to \(\mathrm{Normal}(0, 2)\) on \(\log a_j\), putting the bulk of the mass on \(a_j \in (e^{-4}, e^{4})\) on the positive scale.
- Half-Cauchy \(\mathrm{Cauchy}(0, 1)\) truncated to a small positive lower bound on the random-effect SDs: heavy-tailed, permits very small values (near-zero between-subject variability) while still allowing large ones.
10.2 Fitting MAP and comparing to MLE
We then use MAP(LaplaceI()) to fit the model, which adds the log-priors to the Laplace-approximated log-likelihood as a regularizer. To make the regularization actually matter, we downsample the simulated population to its first 50 subjects (irt_pop[1:50]). At that size, some categories of each item are only observed a handful of times, so the likelihood alone leaves several cut-offs and discriminations weakly identified:
fpm_irt_map = fit(irt_ordinal_model_map, irt_pop[1:50], init_irt, MAP(LaplaceI()))[ Info: Checking the initial parameter values. [ Info: The initial negative log-likelihood and its gradient are finite. Check passed. Iter Function value Gradient norm 0 8.188248e+02 1.610446e+03 * time: 3.0994415283203125e-5 1 7.222346e+02 1.282099e+02 * time: 4.3816978931427 2 7.072352e+02 9.952945e+01 * time: 4.498085975646973 3 6.674699e+02 5.595134e+01 * time: 4.617419958114624 4 6.574952e+02 1.277119e+02 * time: 4.736804008483887 5 6.483717e+02 6.208605e+01 * time: 4.860898017883301 6 6.355318e+02 4.758679e+01 * time: 4.985477924346924 7 6.253840e+02 4.853089e+01 * time: 5.113236904144287 8 6.235519e+02 3.084219e+01 * time: 5.235969066619873 9 6.200235e+02 2.329154e+01 * time: 5.353359937667847 10 6.187342e+02 1.943061e+01 * time: 5.9160308837890625 11 6.178400e+02 1.844845e+01 * time: 6.016913890838623 12 6.171348e+02 1.889776e+01 * time: 6.1155219078063965 13 6.160006e+02 2.340980e+01 * time: 6.213999032974243 14 6.148306e+02 1.536122e+01 * time: 6.330228090286255 15 6.139827e+02 2.905593e+01 * time: 6.45032000541687 16 6.135959e+02 1.627549e+01 * time: 6.569782972335815 17 6.134894e+02 8.537191e+00 * time: 6.687136888504028 18 6.133777e+02 5.553852e+00 * time: 6.804409980773926 19 6.130012e+02 5.536909e+00 * time: 7.076196908950806 20 6.124912e+02 6.524996e+00 * time: 7.181751012802124 21 6.119685e+02 5.037053e+00 * time: 7.286895990371704 22 6.117146e+02 3.802199e+00 * time: 7.396717071533203 23 6.116419e+02 3.899200e+00 * time: 7.513093948364258 24 6.116002e+02 3.941805e+00 * time: 7.63162088394165 25 6.115004e+02 4.091526e+00 * time: 7.750792026519775 26 6.112839e+02 4.273605e+00 * time: 7.996869087219238 27 6.109031e+02 5.570665e+00 * time: 8.10174298286438 28 6.105850e+02 8.064348e+00 * time: 8.20870304107666 29 6.102714e+02 9.233028e+00 * time: 8.313905954360962 30 6.101787e+02 2.244962e+00 * time: 8.430595874786377 31 6.101709e+02 2.268945e+00 * time: 8.547758102416992 32 6.101416e+02 2.302496e+00 * time: 8.664813995361328 33 6.101108e+02 2.180672e+00 * time: 8.806199073791504 34 6.100670e+02 1.700142e+00 * time: 8.892472982406616 35 6.100385e+02 1.188981e+00 * time: 8.974653005599976 36 6.100237e+02 2.170499e+00 * time: 9.06362795829773 37 6.100148e+02 2.177115e+00 * time: 9.161412000656128 38 6.100083e+02 1.186545e+00 * time: 9.295752048492432 39 6.100030e+02 1.047516e+00 * time: 9.382628917694092 40 6.099957e+02 1.880258e+00 * time: 9.467108964920044 41 6.099813e+02 3.866894e+00 * time: 9.554481983184814 42 6.099512e+02 6.042576e+00 * time: 9.69982099533081 43 6.098984e+02 7.265061e+00 * time: 9.788635015487671 44 6.098341e+02 6.341669e+00 * time: 9.872602939605713 45 6.097912e+02 3.792921e+00 * time: 9.95626187324524 46 6.097807e+02 8.147522e-01 * time: 10.112936973571777 47 6.097788e+02 7.833767e-01 * time: 10.213021993637085 48 6.097757e+02 2.415934e+00 * time: 10.300143957138062 49 6.097686e+02 4.666322e+00 * time: 10.38418197631836 50 6.097530e+02 7.293401e+00 * time: 10.518183946609497 51 6.097257e+02 8.942912e+00 * time: 10.603687047958374 52 6.097002e+02 6.912695e+00 * time: 10.687766075134277 53 6.096896e+02 2.285373e+00 * time: 10.770373106002808 54 6.096875e+02 1.023681e-01 * time: 10.903662919998169 55 6.096875e+02 1.600101e-01 * time: 10.98941707611084 56 6.096875e+02 4.417520e-02 * time: 11.074667930603027 57 6.096875e+02 3.950546e-03 * time: 11.15807294845581 58 6.096875e+02 6.106902e-05 * time: 11.287482976913452
FittedPumasModel
Dynamical system type: No dynamical model
Number of subjects: 50
Observation records: Active Missing
score1: 250 0
score2: 250 0
Total: 500 0
Number of parameters: Constant Optimized
0 10
Estimation method: MAP (LaplaceI)
--------------------
Estimate
--------------------
μf -0.65986
μs 0.051183
ωf 0.89282
ωs 0.045609
θb1₁ 0.52516
θb1₂ 0.035213
θb2₁ -1.5777
θb2₂ 0.6691
θb2₃ 0.30676
a_free₁ 0.50861
--------------------We can compare the MAP estimates (50 subjects, priors active) to the earlier LaplaceI estimates (full 300 subjects, MLE) by appending them to ci_df:
map_est = coef(fpm_irt_map)
ci_df.map_50 = vcat(
[map_est.μf, map_est.μs, map_est.ωf, map_est.ωs],
map_est.θb1,
map_est.θb2,
map_est.a_free,
)
simple_table(ci_df)| parameter | description | constant | estimate | se | relative_se | ci_lower | ci_upper | true_value | init_value | map_50 |
| μf | μf | false | -0.603 | 0.0957 | 0.159 | -0.79 | -0.415 | -0.3 | 0.5 | -0.66 |
| μs | μs | false | 0.0448 | 0.00601 | 0.134 | 0.033 | 0.0566 | 0.04 | 0.2 | 0.0512 |
| ωf | ωf | false | 0.842 | 0.0983 | 0.117 | 0.649 | 1.03 | 0.8 | 0.4 | 0.893 |
| ωs | ωs | false | 0.0388 | 0.00864 | 0.223 | 0.0218 | 0.0557 | 0.02 | 0.1 | 0.0456 |
| θb1₁ | θb1₁ | false | 0.449 | 0.0477 | 0.106 | 0.355 | 0.542 | 0.405 | -0.5 | 0.525 |
| θb1₂ | θb1₂ | false | 0.149 | 0.0818 | 0.549 | -0.0114 | 0.309 | 0 | -0.5 | 0.0352 |
| θb2₁ | θb2₁ | false | -1.48 | 0.223 | 0.15 | -1.92 | -1.04 | -1 | 0 | -1.58 |
| θb2₂ | θb2₂ | false | 0.476 | 0.144 | 0.302 | 0.194 | 0.758 | 0.262 | -0.5 | 0.669 |
| θb2₃ | θb2₃ | false | 0.111 | 0.15 | 1.35 | -0.183 | 0.404 | -0.105 | -0.5 | 0.307 |
| a_free₁ | a_free₁ | false | 0.602 | 0.0844 | 0.14 | 0.437 | 0.768 | 0.75 | 0.5 | 0.509 |
The MAP estimates on 50 subjects should sit close to the true values for the well-identified parameters and shrink toward the prior mode for the weakly identified ones.
10.3 Marginal log-likelihood diagnostic
A useful diagnostic is to evaluate the marginal log-likelihood on the same 50-subject dataset at both the MAP estimate and the true parameters:
ll_map = loglikelihood(irt_ordinal_model_map, irt_pop[1:50], map_est, LaplaceI())
ll_true = loglikelihood(irt_ordinal_model_map, irt_pop[1:50], true_params, LaplaceI())
(; ll_map, ll_true, Δll = ll_map - ll_true)(ll_map = -597.9485464513559, ll_true = -601.2558270522167, Δll = 3.307280600860736)In this case, the log-likelihood of the MAP estimate is similar to that of the true parameters on the 50-subject dataset, which implies a good fit to the observed data. Notice that the log-likelihood of the MAP estimate is actually higher than that of the true parameters on this dataset. This is because the true parameters, although they generated the data, are not guaranteed to be the best fit to any particular subset of it (especially a small one) due to sampling variability.
11 Markov chain Monte Carlo (MCMC)
To fully quantify the uncertainty in the parameters of the model in the presence of sparse data, a full Bayesian fit is recommended. Pumas has a No-U-Turn Sampler (NUTS)-based MCMC sampler. We can simply use the BayesMCMC algorithm with the same model definition as for MAP. Here, we start the sampler from the MAP estimate:
fpm_irt_bayes = fit(
irt_ordinal_model_map,
irt_pop[1:50],
coef(fpm_irt_map), # use the MAP estimate as the MCMC initial value
BayesMCMC(; nsamples = 2_000, nadapts = 1_000, nchains = 4),
)
fpm_irt_bayes_samples = discard(fpm_irt_bayes; burnin = 1_000)┌ Warning: Only a single thread available: sampling with a single thread. └ @ Pumas ~/run/_work/PumasTutorials.jl/PumasTutorials.jl/custom_julia_depot/packages/Pumas/btqOJ/src/estimation/bayes/nuts.jl:239 [ Info: Checking the initial parameter values. [ Info: The initial log probability and its gradient are finite. Check passed. [ Info: Checking the initial parameter values. [ Info: The initial log probability and its gradient are finite. Check passed. [ Info: Checking the initial parameter values. [ Info: The initial log probability and its gradient are finite. Check passed. [ Info: Checking the initial parameter values. [ Info: The initial log probability and its gradient are finite. Check passed.
Chains MCMC chain (1000×10×4 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 4
Samples per chain = 1000
Wall duration = 153.93 seconds
Compute duration = 152.82 seconds
parameters = μf, μs, ωf, ωs, θb1₁, θb1₂, θb2₁, θb2₂, θb2₃, a_free₁
Use `describe(chains)` for summary statistics and quantiles.11.1 Convergence diagnostics
The returned object is a collection of MCMC chains rather than a single point estimate. The standard summaries are obtained with describe (mean, std, ESS, R̂, quantiles, etc.):
describe(fpm_irt_bayes_samples)Chains MCMC chain (1000×10×4 Array{Float64, 3}):
Iterations = 1:1:1000
Number of chains = 4
Samples per chain = 1000
Wall duration = 153.93 seconds
Compute duration = 152.82 seconds
parameters = μf, μs, ωf, ωs, θb1₁, θb1₂, θb2₁, θb2₂, θb2₃, a_free₁
Summary Statistics
parameters mean std mcse ess_bulk ess_tail rhat ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
μf -0.6888 0.2504 0.0049 2639.1410 2742.7259 1.0016 ⋯
μs 0.0530 0.0180 0.0004 2616.8937 2619.5156 1.0032 ⋯
ωf 1.0104 0.2495 0.0069 1308.4742 1811.8742 1.0023 ⋯
ωs 0.0555 0.0237 0.0010 596.2368 751.2739 1.0074 ⋯
θb1₁ 0.5484 0.1204 0.0020 3541.1312 2648.9642 1.0012 ⋯
θb1₂ 0.0427 0.1984 0.0027 5651.1759 3086.1739 1.0041 ⋯
θb2₁ -1.8487 0.5915 0.0158 1577.1445 1709.7115 1.0030 ⋯
θb2₂ 0.8085 0.3068 0.0084 1353.7266 1837.3992 1.0038 ⋯
θb2₃ 0.4368 0.3255 0.0085 1470.5661 1746.9652 1.0027 ⋯
a_free₁ 0.4593 0.1471 0.0039 1371.0046 1593.8607 1.0027 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
μf -1.1991 -0.8550 -0.6844 -0.5217 -0.2171
μs 0.0186 0.0409 0.0527 0.0649 0.0887
ωf 0.5631 0.8419 0.9964 1.1684 1.5397
ωs 0.0073 0.0399 0.0560 0.0714 0.1026
θb1₁ 0.3046 0.4693 0.5485 0.6307 0.7839
θb1₂ -0.3615 -0.0895 0.0506 0.1818 0.4042
θb2₁ -3.2920 -2.1680 -1.7623 -1.4338 -0.9213
θb2₂ 0.2396 0.6039 0.7993 1.0062 1.4186
θb2₃ -0.1997 0.2253 0.4292 0.6434 1.1168
a_free₁ 0.2296 0.3589 0.4404 0.5384 0.7966
Before reading anything off the posterior, glance at the rhat and ess_* columns. As a rule of thumb, rhat should be very close to 1.0 (certainly below 1.01) for every parameter, and the effective sample sizes (both bulk and tail) should be at least a few hundred per parameter. If rhat is much above 1.01 or ESS is much smaller than the nominal sample count, the chains have not mixed and the credible intervals below should not be trusted; rerun with more nsamples / nadapts, or revisit the model parameterization and priors.
11.2 Credible intervals
We can get 95% credible intervals for the parameters from the quantiles of the marginal posteriors:
quantile_df = DataFrame(quantile(fpm_irt_bayes_samples))
quantile_df.true_value = vcat(
[true_params.μf, true_params.μs, true_params.ωf, true_params.ωs],
true_params.θb1,
true_params.θb2,
true_params.a_free,
)
simple_table(quantile_df)| parameters | 2.5% | 25.0% | 50.0% | 75.0% | 97.5% | true_value |
| :μf | -1.2 | -0.855 | -0.684 | -0.522 | -0.217 | -0.3 |
| :μs | 0.0186 | 0.0409 | 0.0527 | 0.0649 | 0.0887 | 0.04 |
| :ωf | 0.563 | 0.842 | 0.996 | 1.17 | 1.54 | 0.8 |
| :ωs | 0.00733 | 0.0399 | 0.056 | 0.0714 | 0.103 | 0.02 |
| :θb1₁ | 0.305 | 0.469 | 0.548 | 0.631 | 0.784 | 0.405 |
| :θb1₂ | -0.362 | -0.0895 | 0.0506 | 0.182 | 0.404 | 0 |
| :θb2₁ | -3.29 | -2.17 | -1.76 | -1.43 | -0.921 | -1 |
| :θb2₂ | 0.24 | 0.604 | 0.799 | 1.01 | 1.42 | 0.262 |
| :θb2₃ | -0.2 | 0.225 | 0.429 | 0.643 | 1.12 | -0.105 |
| :a_free₁ | 0.23 | 0.359 | 0.44 | 0.538 | 0.797 | 0.75 |
If the true parameters are not mostly contained in their 95% credible intervals, this could imply that the priors are not weak enough to let the data speak for itself. For more on Bayesian estimation in Pumas, see the Bayesian tutorials in Pumas.
12 Extensions
The model we have built and fitted is intentionally minimal: a linear scalar latent ability \(f_i(t)\) with OrderedLogit observations on two items, no covariates, and no treatment effect. There are at least four ways to extend the IRT model presented here. Each extension requires only a localized change to the model, with the rest of the workflow (simobs, read_pumas, fit(..., LaplaceI()) / MAP(LaplaceI()) / BayesMCMC) unchanged.
12.1 Different disease progression function forms
The linear model \(f_i(t) = f_{0,i} + s_i \, t\) assumes the latent ability changes at a constant rate. In many clinical settings, the disease severity saturates (it changes quickly at first and then approaches a plateau), and a natural parametric alternative is the EMAX form
\[ f_i(t) = f_{0,i} + E_{\max,i} \cdot \frac{t}{EC_{50} + t}, \]
where \(f_{0,i}\) is the subject-specific baseline ability, \(E_{\max,i}\) is the maximum change from baseline (so the asymptote as \(t \to \infty\) is \(f_{0,i} + E_{\max,i}\)), and \(EC_{50} > 0\) is the time at which half of the maximum change is reached. The progression rate at \(t = 0\) is \(E_{\max,i} / EC_{50}\), and the curve approaches its plateau monotonically. Other parametric forms (exponential, Weibull, indirect-response, mechanistic ODE-based, etc.) slot in the same way. Each is a localized change to the @pre block, optionally the @dynamics block, and their associated parameters in @param / @random; the rest of the model (item cut-offs, discriminations, OrderedLogit observation, fitting with LaplaceI) is unaffected.
Sometimes no fixed parametric form is plausible; for example, the trajectory may be non-monotonic and its structure is difficult to specify a priori. In such cases, the structural part of \(f_i(t)\) can be replaced by a neural network:
\[ f_i(t) = \mathrm{NN}_{\boldsymbol{\phi}}\!\left(t,\; \boldsymbol{\eta}_i,\; \boldsymbol{x}_i\right), \]
where \(\boldsymbol{\eta}_i\) are subject-specific random effects, \(\boldsymbol{x}_i\) are covariates, and \(\boldsymbol{\phi}\) are the (population-level) network weights, turning the model to a DeepNLME model. This lets a neural component be embedded directly in @pre (using DeepPumas) alongside the usual @param / @random / @derived machinery. From the rest of the IRT model’s point of view, nothing changes (f is still just a scalar function of time and individual effects), so the same OrderedLogit observation block, simobs, read_pumas, and fit(..., LaplaceI()) workflow apply.
12.2 Covariate and treatment effects
The latent ability in our model depends only on time and subject-specific random effects. In practice, \(f_i(t)\) is often also driven by baseline covariates (such as age, sex, disease subtype, or a baseline biomarker) and by treatment (active vs. placebo, dose, regimen). These can be defined either directly in @pre or using an ODE in the @dynamics block, e.g., using a pharmacokinetic model to drive the disease progression via an exposure-response relationship. The key point is that the IRT-specific parts of the model (item cut-offs, discriminations, OrderedLogit observation) are unaffected by the addition of covariates and treatment effects, so the same workflow applies. If complex covariate relationships are expected, a neural network can be used to flexibly model the covariate effects on the latent trajectory.
12.3 Bounded integer model
The extensions above all leave the OrderedLogit likelihood untouched. The next extension is to replace that likelihood. An alternative is the bounded integer (BI) model, in which the observed score \(y_{i,j} \in \{1, \dots, K_j\}\) arises by discretizing a latent continuous variable on the score scale through a Gaussian CDF, rather than a logistic function as in the OrderedLogit.
12.4 Multi-dimensional IRT
The fourth extension touches the latent state itself. The model above assumes a single scalar latent ability \(f_i(t)\). When the items are driven by multiple, partially independent latent dimensions (for instance, motor function and cognition, or fatigue and pain), the unidimensional assumption (cf. the warning at the top of this tutorial) is violated and a multi-dimensional IRT model is appropriate. The latent state becomes a vector \(\boldsymbol{f}_i(t) \in \mathbb{R}^D\), and each item’s OrderedLogit threshold reads off a linear combination of the latent dimensions via item-specific loadings \(\boldsymbol{w}_j \in \mathbb{R}^D\):
\[ \eta_{j,k} \;=\; a_j \left(b_{j,k} - \boldsymbol{w}_j^\top \boldsymbol{f}_i(t)\right). \]
The multi-dimensional latent trajectory \(\boldsymbol{f}_i(t)\) can be built in two complementary ways:
- Traditional disease progression models. Each component of \(\boldsymbol{f}_i(t)\) is given its own structural form (linear, EMAX, exponential, ODE-based, …), with its own fixed and random effects. This is the natural choice when each latent dimension corresponds to a clinically interpretable construct with a well-understood time course.
- DeepNLME (using DeepPumas). \(\boldsymbol{f}_i(t) = \mathrm{NN}_{\boldsymbol{\phi}}(t, \boldsymbol{\eta}_i, \boldsymbol{x}_i)\) is the output of a neural network with \(D\) outputs. This is well suited to high-dimensional or multi-modal trajectories where no single parametric form fits well across patients, or to settings with many covariates whose interactions with time are non-trivial.
In both cases, the rest of the IRT model (item cut-offs, discriminations, OrderedLogit observation) is unchanged from the unidimensional case.