rerun_dgp <- FALSE
refit_model <- FALSE
library(pks)
library(rstan)
options(mc.cores = 4)
set.seed(1617)
## For plotting with Rgraphviz
myattrs <- list(
graph = list(rankdir = "BT", bgcolor = "white"),
edge = list(arrowsize = NULL),
node = list(fillcolor = "lightblue", fontsize = 10,
shape = "rectangle", fixedsize = FALSE)
)Knowledge structure theory (KST) seeks to provide procedures for the effective diagnosis of a student’s knowledge state in a certain domain, such as algebra, geometry, physics, or statistics (Doignon and Falmagne 1999; Falmagne et al. 1990; Heller and Stefanutti 2024).
The basic notions in KST are
the knowledge domain: a finite set \(Q\) of (dichotomous) problems (also called items);
the knowledge state: the subset \(K \subseteq Q\) of problems in the domain \(Q\) a student masters;
the knowledge structure: a collection \(\mathcal{K}\) of knowledge states that contains at least the empty set \(\emptyset\) and the set \(Q\).
Only a subset of all possible response patterns \(\mathcal{R}\) are also knowledge states in \(\mathcal{K}\).
Lakshminarayan (1995) presented 959 undergraduate students with the five problems in high school geometry shown in Figure 1.
With the knowledge domain \(Q = \{a,b,c,d,e\}\) consisting of five items, there are \(2^{|Q|} = 32\) possible response patterns ranging from no problem solved (\(00000\)) to all problems solved (\(11111\)). Not all of them are plausible knowledge states because solving a more difficult problem may require that a simpler problem is solved first. Doignon and Falmagne (1999) introduced a knowledge structure that is used here for illustration; it consists of \(|\mathcal{K}| = 9\) knowledge states (Figure 2).
data(DoignonFalmagne7, package = "pks") # provides knowledge structure
K <- DoignonFalmagne7$K
rownames(K) <- as.pattern(K, useNames = TRUE)
rownames(K)[rownames(K) == "abcde"] <- "Q"
plot(relations::as.relation(is.subset(K)), attrs = myattrs, main = "")
The knowledge structure in Figure 2 is special in two ways: First, it is closed under union and thus called a knowledge space. Second, it is well-graded: it can be decomposed into learning paths along which a student can traverse – by learning how to solve one problem at a time – from a state of ignorance \(\emptyset\) to full mastery \(Q\). Consequently, \(\mathcal{K}\) is downgradable: For any nonempty state \(K\) there is some problem \(q\) such that \(K \smallsetminus \{q\} \in \mathcal{K}\). The problems a student in state \(K\) is ready to learn at each step, \[ K^\mathcal{O} = \{q \notin K \mid K \cup \{q\} \in \mathcal{K}\}, \] is called the outer fringe of \(K\) (see Figure 3).
rel <- relations::as.relation(is.subset(K))
mycol <- setNames(rep("white", 9),
unlist(relations::relation_domain(rel)[[1]]))
mycol["{}"] <-
mycol["b"] <-
mycol["ab"] <-
mycol["abc"] <-
mycol["abcd"] <-
mycol["Q"] <- "lightblue"
plot(rel,
nodeAttrs = list(fillcolor = mycol),
attrs = myattrs, main = "")
mycol <- setNames(rep("white", 9),
unlist(relations::relation_domain(rel)[[1]]))
mycol["ab"] <- "lightblue"
mycol["abd"] <-
mycol["abc"] <- "lightgreen"
plot(rel,
nodeAttrs = list(fillcolor = mycol),
attrs = myattrs, main = "")
In order to formulate a probabilistic model, a distinction is introduced between knowledge state and response pattern (Falmagne and Doignon 1988a, 1988b). A knowledge state \(K\) is a latent class, and the response pattern \(R\) is a manifest indicator of the knowledge state. As a consequence, this gives rise to two types of response errors,
a careless error, which occurs if the response is incorrect although the problem is contained in \(K\);
a lucky guess, which occurs if the response is correct although the problem is not contained in \(K\).
The latent knowledge states themselves are considered to be random variables with a probability distribution. Formally, a probabilistic knowledge structure consists of
a knowledge structure \(\mathcal{K}\) on a knowledge domain \(Q\) of problems;
the conditional probabilities \(P(R \mid K)\) to observe response pattern \(R\) given state \(K\);
the marginal distribution \(P(K)\) on the knowledge states \(K \in \mathcal{K}\).
The probability of the response pattern \(R \in \mathcal{R} = 2^Q\) then is predicted by \[ P(R) = \sum_{K \in \mathcal{K}} P(R \mid K) \cdot P(K), \] and the response frequencies \(y = (N_R)_{R \in \mathcal{R}}\) follow a multinomial likelihood: \[ p(y \mid \theta) = \prod_{R \in \mathcal{R}} P(R)^{N_R}. \]
Two prevalent models that make specific assumptions about the conditional and marginal probabilities are the basic local independence model and the simple learning model.
For the conditional distribution \(P(R \mid K)\), the BLIM assumes that given the knowledge state \(K\) of a student, responses are stochastically independent over problems (local independence). This implies that the response to each problem \(q\) only depends on the probabilities of a careless error \(\beta_q\) and a lucky guess \(\eta_q\). Accordingly, the probability of the response pattern \(R\) given the knowledge state \(K\) becomes \[ P(R \mid K) = \prod_{q \in K \setminus R} \beta_q \cdot \prod_{q \in K \cap R} (1 - \beta_q) \cdot \prod_{q \in R \setminus K} \eta_q \cdot \prod_{q \in \bar{R} \cap \bar{K}} (1 - \eta_q). \] For example, the probability that a student in state \(\{a,b\}\) solves the problems \(b\) and \(d\) is \[ P(01010 \mid \{a,b\}) = \beta_a \cdot (1 - \beta_b) \cdot (1 - \eta_c) \cdot \eta_d \cdot (1 - \eta_e). \] For the marginal distribution \(P(K)\), the BLIM does not imply any restrictions; for each knowledge state \(K\) it introduces a parameter \(\pi_K\) with \(\pi_K \in (0, 1)\) and \(\sum_K \pi_K = 1\). Thus, in total it has \(2 \cdot |Q| + |\mathcal{K}| - 1\) free parameters.
Just like the BLIM, the SLM assumes local independence for the conditional distribution of responses given a knowledge state \(P(R \mid K)\). On top of that, the SLM establishes restrictions on the marginal distribution of knowledge states \(P(K)\). It requires \(\mathcal{K}\) to be a well-graded knowledge space, that is a knowledge space that is downgradable. For each problem, it introduces a solvability parameter \(g_q\), and the probability of a knowledge state \(K\) becomes \[ P(K) = \prod_{q \in K} g_q \! \prod_{q \in K^\mathcal{O}} \! (1 - g_q). \] This is the probability of mastering all problems in \(K\) while not mastering any problems accessible from \(K\). Accessible here means: problems that are in the outer fringe \(K^\mathcal{O}\) of \(K\). For example, with the well-graded knowledge space in Figure 3, the probability of being in state \(\{a,b\}\) is \[ P(\{a,b\}) = g_a \cdot g_b \cdot (1 - g_c) \cdot (1 - g_d). \] In total, the SLM has \(3 \cdot |Q|\) free parameters.
Using the simulate function from the pks package, we generate response patterns for 2000 students using the knowledge structure in Figure 2 and the BLIM as data generating process (dgp). We assume a fixed distribution over the knowledge states and set all careless error and guessing probabilities to 0.1.
dgp <- function(K) {
nitem <- ncol(K)
m0 <- list(
P.K = setNames(c(0.15, 0.1, 0.1, 0.08, 0.1, 0.12, 0.1, 0.1, 0.15),
as.pattern(K)),
beta = rep(0.1, nitem),
eta = rep(0.1, nitem),
K = K,
ntotal = 2000
)
class(m0) <- "blim"
y <- simulate(m0, method = "multinomial", zeropad = TRUE)
attr(y, "dgp") <- m0
y
}N.R <- dgp(K)
dput(N.R, file = "data/data_dgp_blim.R")N.R <- dget("data/data_dgp_blim.R")
datlist <- list(
R = as.binmat(N.R), # distinct response patterns
K = K, # knowledge states
Ko = getKFringe(K), # outer fringes
y = c(N.R), # by-pattern multinomial count
ycat = match(rep(names(N.R), N.R), names(N.R)), # by-subj response category
nstate = nrow(K),
nitem = ncol(K),
N = sum(N.R)
)
datlist$npattern <- nrow(datlist$R)
datlist$y |> print()00000 10000 01000 11000 00100 10100 01100 11100 00010 10010 01010 11010 00110
207 141 152 173 23 36 37 179 27 33 33 160 8
10110 01110 11110 00001 10001 01001 11001 00101 10101 01101 11101 00011 10011
22 17 166 18 25 14 33 3 17 19 145 5 10
01011 11011 00111 10111 01111 11111
8 47 0 15 24 203 datlist$ycat |> table()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
207 141 152 173 23 36 37 179 27 33 33 160 8 22 17 166 18 25 14 33
21 22 23 24 25 26 27 28 30 31 32
3 17 19 145 5 10 8 47 15 24 203 The list of data that is passed on to Stan includes two versions of the generated responses: y is a named vector that contains the multinomial counts \(N_R\), some may be zero, but its length always is npattern, which is the number of all possible distinct response patterns; ycat is a vector that contains the response category of each student as a position index in y, it must be a number between 1 and npattern. Whenever some response patterns were not observed, the response frequency vector y should be completed with zeros, as is achieved by the zeropad = TRUE setting above. The models below do not account for missingness of single item responses (see Chiusole et al. 2015).
In Stan, the BLIM can be written as follows:
functions {
matrix getPRK(vector beta, vector eta, matrix K, matrix R) { // P(R|K)
return exp(
(1 - R) * diag_pre_multiply(log( beta), K' ) +
R * diag_pre_multiply(log1m(beta), K' ) +
R * diag_pre_multiply(log( eta), (1 - K)') +
(1 - R) * diag_pre_multiply(log1m( eta), (1 - K)')
);
}
}
data {
int<lower=1> N; // number of students
int<lower=1> nstate; // number of knowledge states
int<lower=1> nitem; // number of problems
int<lower=1> npattern; // number of distinct patterns
array[npattern] int<lower=0> y; // response frequencies
array[N] int<lower=1, upper=npattern> ycat; // response categories
matrix<lower=0, upper=1>[nstate, nitem] K; // knowledge structure matrix
matrix<lower=0, upper=1>[npattern, nitem] R; // response pattern matrix
}
parameters {
vector<lower=0, upper=0.5>[nitem] beta; // error probabilities
vector<lower=0, upper=0.5>[nitem] eta; // guess probabilities
simplex[nstate] pi; // state probabilities
}
model {
// Priors
beta ~ beta(2, 9);
eta ~ beta(2, 9);
pi ~ dirichlet(rep_row_vector(1.0, nstate));
// Likelihood
y ~ multinomial(getPRK(beta, eta, K, R) * pi);
}
generated quantities {
vector[npattern] PR = getPRK(beta, eta, K, R) * pi; // P(R)
array[npattern] int yrep = multinomial_rng(PR, N);
vector[N] log_lik;
for (i in 1:N) {
log_lik[i] = categorical_lpmf(ycat[i] | PR);
}
}The conditional distribution \(P(R \mid K)\) is implemented as a separate getPRK function, which is called in the model block to build the multinomial likelihood. The weakly informative priors for the \(\beta_q\) and \(\eta_q\) parameters are beta distributions assuming low error rates. In addition, the error rates are restricted to fall below 0.5; this is to avoid nonsensical behavior, such as a correct response from guessing being more likely than from mastering the problem. The prior for the \(\pi_K\) parameters is a diffuse Dirichlet distribution, as both high and low probabilities for the knowledge states are plausible a priori.
In summary, the Bayesian BLIM can be written as: \[\begin{align*} \beta_q &\sim beta(2, 9), \quad \beta_q \in (0, 0.5)\\ \eta_q &\sim beta(2, 9), \quad \eta_q \in (0, 0.5)\\ \pi_K &\sim dirichlet(\mathbf{\alpha}), \quad \mathbf{\alpha} = (1, \dots, 1)\\ y &\sim multinomial(P(R)), \quad y = (N_R)_{R \in \mathcal{R}}, \quad \sum_{R \in \mathcal{R}} N_R = N\\ P(R) &= \sum_{K \in \mathcal{K}} P(R \mid K) \cdot P(K)\\ P(R \mid K) &= \prod_{q \in K \setminus R} \beta_q \cdot \prod_{q \in K \cap R} (1 - \beta_q) \cdot \prod_{q \in R \setminus K} \eta_q \cdot \prod_{q \in \bar{R} \cap \bar{K}} (1 - \eta_q)\\ P(K) &= \pi_K, \quad \pi_K \in (0, 1), \quad \sum_{K \in \mathcal{K}} \pi_K = 1, \end{align*}\] where \(R\) denotes the response pattern, \(N_R\) its frequency, and \(K\) a knowledge state.
Stan fits the model to the generated responses and draws posterior samples for its parameters:
## Code takes time to run; instead use saved output
blim1 <- stan("stan/blim.stan", data = datlist, seed = 1634)
saveRDS(blim1, "saved_fit/blim1.rds")blim1 <- readRDS(file = "saved_fit/blim1.rds")
print(blim1, pars = c("beta", "eta", "pi"), probs = c(.025, .975))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 97.5% n_eff Rhat
beta[1] 0.10 0 0.01 0.07 0.12 3383 1
beta[2] 0.09 0 0.01 0.07 0.12 3482 1
beta[3] 0.14 0 0.03 0.09 0.19 3556 1
beta[4] 0.20 0 0.10 0.03 0.38 1340 1
beta[5] 0.20 0 0.10 0.03 0.40 1397 1
eta[1] 0.19 0 0.10 0.03 0.39 1262 1
eta[2] 0.19 0 0.10 0.03 0.39 1852 1
eta[3] 0.11 0 0.02 0.08 0.15 2926 1
eta[4] 0.12 0 0.02 0.09 0.16 3509 1
eta[5] 0.10 0 0.01 0.07 0.13 3274 1
pi[1] 0.20 0 0.05 0.13 0.31 1459 1
pi[2] 0.09 0 0.04 0.01 0.15 1104 1
pi[3] 0.08 0 0.03 0.01 0.15 1447 1
pi[4] 0.05 0 0.03 0.00 0.10 1226 1
pi[5] 0.09 0 0.02 0.04 0.14 1454 1
pi[6] 0.12 0 0.02 0.08 0.17 1676 1
pi[7] 0.09 0 0.04 0.01 0.15 1352 1
pi[8] 0.08 0 0.04 0.01 0.15 1291 1
pi[9] 0.20 0 0.05 0.13 0.31 1307 1
Samples were drawn using NUTS(diag_e) at Mon Sep 21 11:59:30 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).The convergence diagnostics and effective sample sizes look OK for all parameters. Figure 4 shows the posterior distribution of the BLIM parameters and their true values.
posteriorBlimPlot <- function(samples, npar = 5 + 5 + 9,
N.R = NULL, dgp = FALSE) {
with(samples,
boxplot(cbind(beta, eta, pi)[, npar:1], horizontal = TRUE,
border = "gray", axes = FALSE, xlab = "Parameter value")
)
points(rev(coef(attr(N.R, "dgp"))), 1:npar, pch = 4, col = "darkblue")
axis(1)
axis(2, at = 1:npar, rev(c(
expression(beta[a]),
expression(beta[b]),
expression(beta[c]),
expression(beta[d]),
expression(beta[e]),
expression(eta[a]),
expression(eta[b]),
expression(eta[c]),
expression(eta[d]),
expression(eta[e]),
expression(pi["{}"]),
expression(pi[a]),
expression(pi[b]),
expression(pi[ab]),
expression(pi[abc]),
expression(pi[abd]),
expression(pi[abcd]),
expression(pi[abce]),
expression(pi[Q])
)), las = 1)
legend("topright", if(dgp) "BLIM (dgp)" else "BLIM", bty = "n")
box()
}
par(mai = c(.6, .5, .02, .03),
mgp = c(2, .7, 0))
posteriorBlimPlot(extract(blim1), N.R = N.R, dgp = TRUE)
legend("bottomright", inset = c(0, -.02),
legend = c("true", "estimated"), pch = c(4, 16),
col = c("darkblue", "gray"), bty = "n", ncol = 2, xpd = NA)
Parameter recovery is not perfect, but certain tradeoffs between parameters due to nonidentifiability were expected. As a consequence of the knowledge structure being a well-graded knowledge space, some parameters of the BLIM are not identifiable (Heller 2017; Stefanutti 2012).
c(
space = pks:::is.knowledgespace(K),
downgradable = is.downgradable(K),
forward = names(which(is.forward.graded(K))),
backward = names(which(is.backward.graded(K)))
) space downgradable forward1 forward2 backward1 backward2
"TRUE" "TRUE" "a" "b" "d" "e" The knowledge structure (Figure 2) is forward graded in problems \(a\) and \(b\) (they can be added to any state and still yield a state) and backward graded in \(d\) and \(e\) (they can be removed from any state and still yield a state). Thus, the error and guess probabilities of these problems trade off with some of the state probabilities. The BLIM’s partial nonidentifiability, however, does not compromise its predictive performance as shown below (Figure 5).
In Stan, the SLM can be written like this:
functions {
matrix getPRK(vector beta, vector eta, matrix K, matrix R) {
return exp(
(1 - R) * diag_pre_multiply(log( beta), K' ) +
R * diag_pre_multiply(log1m(beta), K' ) +
R * diag_pre_multiply(log( eta), (1 - K)') +
(1 - R) * diag_pre_multiply(log1m( eta), (1 - K)')
);
}
vector getPK(row_vector g, matrix K, matrix Ko, int nstate) { // P(K)
vector[nstate] log_pi;
for (i in 1:nstate) {
log_pi[i] = sum(K[i, ] .* log(g)) + sum(Ko[i, ] .* log1m(g));
}
return exp(log_pi);
}
}
data {
int<lower=1> N;
int<lower=1> nstate;
int<lower=1> nitem;
int<lower=1> npattern;
array[npattern] int<lower=0> y;
array[N] int<lower=1, upper=npattern> ycat;
matrix<lower=0, upper=1>[nstate, nitem] K;
matrix<lower=0, upper=1>[nstate, nitem] Ko; // outer fringes
matrix<lower=0, upper=1>[npattern, nitem] R;
}
parameters {
vector<lower=0, upper=0.5>[nitem] beta;
vector<lower=0, upper=0.5>[nitem] eta;
row_vector<lower=0, upper=1>[nitem] g;
}
model {
beta ~ beta(2, 9);
eta ~ beta(2, 9);
g ~ beta(5, 1);
y ~ multinomial(getPRK(beta, eta, K, R) * getPK(g, K, Ko, nstate));
}
generated quantities {
vector[npattern] PR = getPRK(beta, eta, K, R) * getPK(g, K, Ko, nstate);
array[npattern] int yrep = multinomial_rng(PR, N);
vector[N] log_lik;
for (i in 1:N) {
log_lik[i] = categorical_lpmf(ycat[i] | PR);
}
}For the marginal state distribution \(P(K)\), the code includes the getPK function that is called in the model block. The solvability parameters \(g_q\) have a beta prior. Its parameters should be chosen to produce an a priori plausible marginal state distribution. The model-implied marginal distribution evaluated at the prior mean can be checked via:
getSlmPK(g = rep(5/6, 5), K = K, Ko = getKFringe(K)) {} a b ab abc abd abcd
0.02777778 0.13888889 0.13888889 0.01929012 0.01607510 0.09645062 0.08037551
abce Q
0.08037551 0.40187757 In summary, the Bayesian SLM can be written as: \[\begin{align*} \beta_q &\sim beta(2, 9), \quad \beta_q \in (0, 0.5)\\ \eta_q &\sim beta(2, 9), \quad \eta_q \in (0, 0.5)\\ g_q &\sim beta(5, 1)\\ y &\sim multinomial(P(R)), \quad y = (N_R)_{R \in \mathcal{R}}, \quad \sum_{R \in \mathcal{R}} N_R = N\\ P(R) &= \sum_{K \in \mathcal{K}} P(R \mid K) \cdot P(K)\\ P(R \mid K) &= \prod_{q \in K \setminus R} \beta_q \cdot \prod_{q \in K \cap R} (1 - \beta_q) \cdot \prod_{q \in R \setminus K} \eta_q \cdot \prod_{q \in \bar{R} \cap \bar{K}} (1 - \eta_q)\\ P(K) &= \prod_{q \in K} g_q \! \prod_{q \in K^\mathcal{O}} \! (1 - g_q), \end{align*}\] where \(R\) denotes the response pattern, \(N_R\) its frequency, \(K\) a knowledge state, and \(K^\mathcal{O}\) its outer fringe.
Fitting in Stan yields:
slm1 <- readRDS(file = "saved_fit/slm1.rds")
print(slm1, pars = c("beta", "eta", "g"), probs = c(.025, .975))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 97.5% n_eff Rhat
beta[1] 0.13 0 0.02 0.10 0.16 3588 1
beta[2] 0.12 0 0.02 0.09 0.15 2941 1
beta[3] 0.14 0 0.03 0.09 0.19 2945 1
beta[4] 0.27 0 0.11 0.05 0.45 2369 1
beta[5] 0.27 0 0.12 0.05 0.46 2360 1
eta[1] 0.11 0 0.05 0.02 0.23 3005 1
eta[2] 0.11 0 0.05 0.02 0.22 3100 1
eta[3] 0.12 0 0.02 0.08 0.15 3208 1
eta[4] 0.12 0 0.02 0.09 0.16 3898 1
eta[5] 0.10 0 0.01 0.07 0.12 3378 1
g[1] 0.78 0 0.02 0.73 0.82 2677 1
g[2] 0.77 0 0.02 0.72 0.81 2530 1
g[3] 0.73 0 0.03 0.67 0.80 2896 1
g[4] 0.73 0 0.14 0.51 0.98 2119 1
g[5] 0.71 0 0.13 0.50 0.97 2005 1
Samples were drawn using NUTS(diag_e) at Mon Sep 21 12:01:20 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Also for the SLM, the convergence diagnostics and effective sample sizes look OK for all parameters.
The MNM serves as a reference for the two substantive models. It can be written like this:
data {
int<lower=1> N;
int<lower=1> npattern;
array[npattern] int<lower=0> y;
array[N] int<lower=1, upper=npattern> ycat;
}
parameters {
simplex[npattern] theta;
}
model {
theta ~ dirichlet(rep_row_vector(1.0, npattern));
y ~ multinomial(theta);
}
generated quantities {
array[npattern] int yrep = multinomial_rng(theta, N);
vector[N] log_lik;
for (i in 1:N) {
log_lik[i] = categorical_lpmf(ycat[i] | theta);
}
}The MNM has one parameter for each distinct response pattern, so in total \(2^{|Q|} - 1\) free parameters. We obtain the following abbreviated posterior summary:
mnm1 <- readRDS(file = "saved_fit/mnm1.rds")
print(mnm1, pars = c("theta[1]", "theta[2]", "theta[32]"),
probs = c(.025, .975))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 97.5% n_eff Rhat
theta[1] 0.10 0 0.01 0.09 0.12 7228 1
theta[2] 0.07 0 0.01 0.06 0.08 5320 1
theta[32] 0.10 0 0.01 0.09 0.11 6056 1
Samples were drawn using NUTS(diag_e) at Mon Sep 21 12:02:42 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Again convergence diagnostics and effective sample sizes look good.
Figure 5 shows the BLIM-generated and predicted frequencies for each distinct response pattern. The BLIM achieves a close fit between observations and predictions, not much worse than the MNM. The SLM’s predictions are further off; this is especially visible for patterns with higher frequency.
postPredPlot <- function(blim, slm, mnm, N.R, dgp = 1) {
par(mfrow = c(1, 3), omi = c(.4, .25, 0, 0), mai = c(.1, .1, .3, .02),
mgp = c(2, .7, 0))
blim$yrep |>
boxplot(border = "gray", horizontal = TRUE, axes = FALSE,
main = if(dgp == 1) "BLIM (dgp)" else "BLIM")
points(N.R, 1:32, pch = 4, col = "darkblue")
axis(1)
axis(2, 1:32, names(N.R), las = 1, tick = FALSE, line = -1,
outer = TRUE, cex.axis = 0.82)
box()
#
slm$yrep |>
boxplot(border = "gray", horizontal = TRUE, axes = FALSE,
main = if(dgp == 2) "SLM (dgp)" else "SLM")
points(N.R, 1:32, pch = 4, col = "darkblue")
axis(1)
axis(2, 1:32, labels = FALSE)
box()
#
mnm$yrep |>
boxplot(border = "gray", horizontal = TRUE, axes = FALSE, main = "MNM")
points(N.R, 1:32, pch = 4, col = "darkblue")
axis(1)
axis(2, 1:32, labels = FALSE)
box()
mtext(paste0("Frequency (N = ", sum(N.R), ")"),
side = 1, line = 1.5, outer = TRUE, cex = .8)
legend("bottomright", inset = c(0, -.12), pch = c(4, 16),
legend = c(if(dgp > 0) "generated" else "observed", "predicted"),
col = c("darkblue", "gray"), bty = "n", ncol = 2, xpd = NA)
}
postPredPlot(extract(blim1), extract(slm1), extract(mnm1), N.R)
When each model computes the observation-wise log-likelihood in the generated quantities block, approximate leave-one-out cross-validation (LOO-CV, Vehtari et al. 2017) is available to compare the models with respect to their ability to predict out-of-sample observations.
l1 <- loo(blim1) |> print()
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -5853.5 38.6
p_loo 13.6 0.2
looic 11707.1 77.2
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.4]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.l2 <- loo(slm1) |> print()
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -5935.5 36.6
p_loo 11.6 0.2
looic 11871.0 73.2
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.3]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.l3 <- loo(mnm1) |> print()
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -5861.2 38.4
p_loo 29.1 1.0
looic 11722.4 76.9
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [1.0, 1.8]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.loo::loo_compare(list(BLIM = l1, SLM = l2, MNM = l3)) model elpd_diff se_diff p_worse diag_diff diag_elpd
BLIM 0.0 0.0 NA
MNM -7.6 3.7 0.98
SLM -82.0 13.0 1.00 For all three models the approximation to LOO-CV is OK. The model most successful at predicting new response frequencies is the BLIM, followed by the MNM. The SLM is the least successful of the three. Remember that all three models were estimated based on the BLIM-generated data.
We generate a new set of response frequencies from a rudimentary SLM. The solvability parameters \(g_q\) are all set to 0.8; this gives rise to a distribution over the knowledge states that is different to the one above, where the BLIM was the data generator. Otherwise, the settings are the same as above.
dgp <- function(K) {
nitem <- ncol(K)
m0 <- list(
beta = rep(0.1, nitem),
eta = rep(0.1, nitem),
g = rep(0.8, nitem),
K = K,
ntotal = 2000
)
m0$P.K <- getSlmPK(g = m0$g, K = K, Ko = getKFringe(K))
class(m0) <- c("slm", "blim")
y <- simulate(m0, method = "multinomial", zeropad = TRUE)
attr(y, "dgp") <- m0
y
}N.R <- dgp(K)
dput(N.R, file = "data/data_dgp_slm.R")N.R <- dget("data/data_dgp_slm.R")
datlist$y <- c(N.R) |> print()00000 10000 01000 11000 00100 10100 01100 11100 00010 10010 01010 11010 00110
83 207 190 86 8 27 22 69 17 40 45 133 5
10110 01110 11110 00001 10001 01001 11001 00101 10101 01101 11101 00011 10011
19 17 178 11 16 24 39 5 20 14 149 0 5
01011 11011 00111 10111 01111 11111
12 46 4 33 54 422 datlist$ycat <- match(rep(names(N.R), N.R), names(N.R))
datlist$N <- sum(N.R)
datlist$ycat |> table()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
83 207 190 86 8 27 22 69 17 40 45 133 5 19 17 178 11 16 24 39
21 22 23 24 26 27 28 29 30 31 32
5 20 14 149 5 12 46 4 33 54 422 slm2 <- readRDS(file = "saved_fit/slm2.rds")
print(slm2, pars = c("beta", "eta", "g"), probs = c(.025, .975))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 97.5% n_eff Rhat
beta[1] 0.10 0 0.01 0.08 0.12 3611 1
beta[2] 0.08 0 0.01 0.07 0.10 3717 1
beta[3] 0.10 0 0.02 0.07 0.13 3392 1
beta[4] 0.17 0 0.06 0.04 0.27 2607 1
beta[5] 0.17 0 0.07 0.03 0.28 2490 1
eta[1] 0.11 0 0.05 0.02 0.22 2101 1
eta[2] 0.11 0 0.05 0.02 0.21 1980 1
eta[3] 0.09 0 0.01 0.06 0.12 3975 1
eta[4] 0.14 0 0.02 0.11 0.17 4060 1
eta[5] 0.10 0 0.01 0.07 0.12 4090 1
g[1] 0.81 0 0.02 0.77 0.84 2257 1
g[2] 0.79 0 0.02 0.76 0.82 2294 1
g[3] 0.83 0 0.02 0.80 0.87 3319 1
g[4] 0.86 0 0.08 0.72 0.99 2565 1
g[5] 0.86 0 0.08 0.71 0.99 2520 1
Samples were drawn using NUTS(diag_e) at Mon Sep 21 12:04:27 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Figure 6 shows the posterior distribution of SLM parameters and their true values. Again, parameter recovery is not perfect. There is a visible tradeoff between the solvability and careless-error parameters for problems \(d\) and \(e\): Higher solvability goes along with more careless errors and vice versa.
posteriorSlmPlot <- function(samples, npar = 5 + 5 + 5,
N.R = NULL, dgp = FALSE) {
with(samples,
boxplot(cbind(beta, eta, g)[, npar:1], horizontal = TRUE,
border = "gray", axes = FALSE, xlab = "Parameter value")
)
points(rev(coef(attr(N.R, "dgp"))), 1:npar, pch = 4, col = "darkblue")
axis(1)
axis(2, at = 1:npar, rev(c(
expression(beta[a]),
expression(beta[b]),
expression(beta[c]),
expression(beta[d]),
expression(beta[e]),
expression(eta[a]),
expression(eta[b]),
expression(eta[c]),
expression(eta[d]),
expression(eta[e]),
expression(g[a]),
expression(g[b]),
expression(g[c]),
expression(g[d]),
expression(g[e])
)), las = 1)
legend("topright", if(dgp) "SLM (dgp)" else "SLM", bty = "n")
box()
}
par(mai = c(.6, .5, .02, .03),
mgp = c(2, .7, 0))
posteriorSlmPlot(extract(slm2), N.R = N.R, dgp = TRUE)
legend("bottomright", inset = c(0, -.02),
legend = c("true", "estimated"), pch = c(4, 16),
col = c("darkblue", "gray"), bty = "n", ncol = 2, xpd = NA)
In addition to the SLM, the BLIM and the MNM were fitted to the new data (not shown, but see code supplement). Figure 7 shows the SLM-generated and predicted frequencies for each distinct response pattern. All three models achieve a close fit between observations and predictions.
postPredPlot(extract(blim2), extract(slm2), extract(mnm2), N.R, dgp = 2)
l1 <- loo(blim2) |> print()
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -5619.4 45.7
p_loo 13.6 0.2
looic 11238.9 91.4
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.8]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.l2 <- loo(slm2) |> print()
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -5617.6 45.8
p_loo 11.7 0.2
looic 11235.2 91.6
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.9, 1.5]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.l3 <- loo(mnm2) |> print()
Computed from 4000 by 2000 log-likelihood matrix.
Estimate SE
elpd_loo -5622.6 45.6
p_loo 29.3 1.1
looic 11245.2 91.2
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [1.2, 2.4]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.loo::loo_compare(list(BLIM = l1, SLM = l2, MNM = l3)) model elpd_diff se_diff p_worse diag_diff diag_elpd
SLM 0.0 0.0 NA
BLIM -1.9 0.5 1.00 |elpd_diff| < 4
MNM -5.0 5.0 0.84
Diagnostic flags present.
See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
or https://mc-stan.org/loo/reference/loo-glossary.html.Now that the SLM has generated the responses that were used for parameter estimation, it also turns out to be the model with the highest predictive accuracy, closely followed by the BLIM. The MNM is worst of the three. Considering the uncertainty as indicated by the se_diff and p_worse values, the differences are not large.
Lakshminarayan (1995, 1996) presented the high school geometry problems shown in Figure 1 to 959 undergraduate students before and after a lesson. The data were split into a training and a held-out set; only the 800 responses in the training data set are reported. In this application, the posttest responses (Figure 8) will be analyzed. They are referred to by Doignon and Falmagne (1999, chap. 8).
data(hsgeometry, package = "pks")
dotchart(angles$N.R[, "posttest"], pch = 4,
main = "High school geometry problems",
xlab = "Response frequency (pre- [o] and posttest [x])")
points(angles$N.R[, "pretest"], 1:32)
mtext("Lakshminarayan (1995)", side = 3, line = 0.5)
N.R <- angles$N.R[, "posttest"] |>
print()00000 00001 00010 00011 00100 00101 00110 00111 01000 01001 01010 01011 01100
12 0 1 0 3 0 4 0 14 1 9 0 24
01101 01110 01111 10000 10001 10010 10011 10100 10101 10110 10111 11000 11001
2 15 2 5 0 2 0 7 0 7 3 23 5
11010 11011 11100 11101 11110 11111
31 10 160 22 288 150 
The knowledge structure proposed by Lakshminarayan (1995) is displayed in Figure 9.
K <- angles$K
rownames(K) <- as.pattern(K, useNames = TRUE)
rownames(K)[rownames(K) == "abcde"] <- "Q"
plot(relations::as.relation(is.subset(K)), attrs = myattrs, main = "")
With an open-answer format like in this problem set, it is plausible to assume that guessing the correct response is impossible. Therefore, all guessing parameters \(\eta_q\) are set to zero. Accordingly, the getPRK function has to be adapted, so a no-guessing BLIM can be written in Stan like this:
functions {
matrix getPRK(vector beta, matrix eta_filter, matrix K, matrix R) {
return exp(
(1 - R) * diag_pre_multiply(log( beta), K') +
R * diag_pre_multiply(log1m(beta), K')
) .* eta_filter;
}
}
data {
int<lower=1> N;
int<lower=1> nstate;
int<lower=1> nitem;
int<lower=1> npattern;
array[npattern] int<lower=0> y;
array[N] int<lower=1, upper=npattern> ycat;
matrix<lower=0, upper=1>[nstate, nitem] K;
matrix<lower=0, upper=1>[npattern, nitem] R;
}
transformed data { // or do it in R and pass as data: R %*% t(1 - K) == 0
matrix<lower=0>[npattern, nstate] eta_filter = R * (1 - K)';
for (i in 1:npattern) { // if item in R but not in K:
for (j in 1:nstate) { // set P(R|K) to zero
eta_filter[i, j] = eta_filter[i, j] > 0 ? 0 : 1;
}
}
}
parameters {
vector<lower=0, upper=0.5>[nitem] beta;
simplex[nstate] pi;
}
model {
beta ~ beta(2, 9);
pi ~ dirichlet(rep_row_vector(1.0, nstate));
y ~ multinomial(getPRK(beta, eta_filter, K, R) * pi);
}
generated quantities {
vector[npattern] PR = getPRK(beta, eta_filter, K, R) * pi;
array[npattern] int yrep = multinomial_rng(PR, N);
vector[N] log_lik;
for (i in 1:N) {
log_lik[i] = categorical_lpmf(ycat[i] | PR);
}
}The above code assumes that when response patterns are missing, their frequencies are set to zero such that the frequency vector is always complete. An alternative strategy is to leave out missing response patterns: R only contains the observed patterns and y their frequencies; npattern \(<
2^{|Q|}\) is their number. Then, the likelihood in the model block should avoid the multinomial distribution statement (which requires a simplex as an argument). Instead, a custom likelihood can be coded like this:
transformed parameters {
vector[npattern] logPR = log(getPRK(beta, eta_filter, K, R) * pi);
}
model {
beta ~ beta(2, 9);
pi ~ dirichlet(rep_row_vector(1.0, nstate));
target += sum(y .* logPR); // custom likelihood
}We pass the data to Stan and obtain posterior samples for the no-guessing versions of the BLIM and the SLM.
datlist <- list(
R = as.binmat(N.R),
K = K,
Ko = getKFringe(K),
y = N.R,
ycat = match(rep(names(N.R), N.R), names(N.R)),
nstate = nrow(K),
nitem = ncol(K),
N = sum(N.R)
)
datlist$npattern <- nrow(datlist$R)blim3 <- readRDS(file = "saved_fit/blim3.rds")
print(blim3, pars = c("beta", "pi"), probs = c(.025, .975))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 97.5% n_eff Rhat
beta[1] 0.04 0 0.01 0.02 0.07 2543 1
beta[2] 0.04 0 0.01 0.02 0.05 4192 1
beta[3] 0.09 0 0.02 0.06 0.12 2537 1
beta[4] 0.16 0 0.03 0.11 0.21 2490 1
beta[5] 0.19 0 0.11 0.03 0.44 1970 1
pi[1] 0.02 0 0.00 0.01 0.03 4989 1
pi[2] 0.01 0 0.00 0.00 0.01 6662 1
pi[3] 0.02 0 0.00 0.01 0.03 5542 1
pi[4] 0.01 0 0.01 0.00 0.02 3648 1
pi[5] 0.00 0 0.00 0.00 0.01 4331 1
pi[6] 0.02 0 0.01 0.01 0.04 4118 1
pi[7] 0.16 0 0.02 0.11 0.20 2582 1
pi[8] 0.01 0 0.01 0.00 0.04 2307 1
pi[9] 0.02 0 0.01 0.00 0.04 2464 1
pi[10] 0.42 0 0.06 0.29 0.52 1826 1
pi[11] 0.31 0 0.05 0.24 0.43 1791 1
Samples were drawn using NUTS(diag_e) at Tue Sep 15 15:22:18 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).slm3 <- readRDS(file = "saved_fit/slm3.rds")
print(slm3, pars = c("beta", "g"), probs = c(.025, .975))Inference for Stan model: anon_model.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 97.5% n_eff Rhat
beta[1] 0.03 0 0.01 0.02 0.06 3141 1
beta[2] 0.04 0 0.01 0.02 0.05 4284 1
beta[3] 0.09 0 0.02 0.06 0.12 2709 1
beta[4] 0.17 0 0.03 0.12 0.23 2345 1
beta[5] 0.26 0 0.12 0.05 0.48 1965 1
g[1] 0.92 0 0.01 0.89 0.95 3309 1
g[2] 0.98 0 0.01 0.97 0.99 4410 1
g[3] 0.96 0 0.02 0.93 0.99 2762 1
g[4] 0.80 0 0.03 0.75 0.86 2392 1
g[5] 0.48 0 0.09 0.35 0.68 1719 1
Samples were drawn using NUTS(diag_e) at Tue Sep 15 15:23:43 2026.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Figure 10 shows the posterior distribution of the BLIM and SLM parameters.
posteriorBlimPlot <- function(samples, npar = 5 + 11,
N.R = NULL, dgp = FALSE) {
with(samples,
boxplot(cbind(beta, pi)[, npar:1], horizontal = TRUE,
xlim = c(1, 16), ylim = 0:1, border = "gray", axes = FALSE)
)
axis(1)
axis(2, at = 1:npar, rev(c(
expression(beta[a]),
expression(beta[b]),
expression(beta[c]),
expression(beta[d]),
expression(beta[e]),
expression(pi["{}"]),
expression(pi[a]),
expression(pi[b]),
expression(pi[ab]),
expression(pi[ad]),
expression(pi[bc]),
expression(pi[abc]),
expression(pi[abd]),
expression(pi[bcd]),
expression(pi[abcd]),
expression(pi[Q])
)), las = 1)
legend("topright", if(dgp) "BLIM (dgp)" else "BLIM", bty = "n")
box()
}
posteriorSlmPlot <- function(samples, npar = 5 + 5,
N.R = NULL, dgp = FALSE) {
with(samples,
boxplot(cbind(beta, g)[, npar:1], at = 1:npar + 6, horizontal = TRUE,
xlim = c(1, 16), ylim = 0:1, border = "gray", axes = FALSE)
)
axis(1)
axis(2, at = 1:npar + 6, rev(c(
expression(beta[a]),
expression(beta[b]),
expression(beta[c]),
expression(beta[d]),
expression(beta[e]),
expression(g[a]),
expression(g[b]),
expression(g[c]),
expression(g[d]),
expression(g[e])
)), las = 1)
legend("topright", if(dgp) "SLM (dgp)" else "SLM", bty = "n")
box()
}
par(mfrow = 1:2, mai = c(.4, .5, .02, .03), omi = c(.3, 0, 0, 0),
mgp = c(2, .7, 0))
posteriorBlimPlot(extract(blim3))
posteriorSlmPlot(extract(slm3))
mtext("Parameter value", side = 1, outer = TRUE)
Figure 11 shows observed response frequencies and their predictions from the BLIM, SLM, and MNM. The BLIM predictions are almost as close to the observations as the MNM predictions. The SLM is further off, especially for response patterns with high frequency.
postPredPlot(extract(blim3), extract(slm3), extract(mnm3), N.R, dgp = 0)
l1 <- loo(blim3) |> print()
Computed from 4000 by 800 log-likelihood matrix.
Estimate SE
elpd_loo -1654.4 37.1
p_loo 11.9 0.8
looic 3308.7 74.3
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.7]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.l2 <- loo(slm3) |> print()
Computed from 4000 by 800 log-likelihood matrix.
Estimate SE
elpd_loo -1683.7 38.1
p_loo 8.7 0.5
looic 3367.5 76.3
------
MCSE of elpd_loo is 0.0.
MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.5]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.l3 <- loo(mnm3) |> print()
Computed from 4000 by 800 log-likelihood matrix.
Estimate SE
elpd_loo -1650.9 35.8
p_loo 20.9 1.8
looic 3301.8 71.5
------
MCSE of elpd_loo is 0.1.
MCSE and ESS estimates assume MCMC draws (r_eff in [1.4, 1.9]).
All Pareto k estimates are good (k < 0.7).
See help('pareto-k-diagnostic') for details.loo::loo_compare(list(BLIM = l1, SLM = l2, MNM = l3)) model elpd_diff se_diff p_worse diag_diff diag_elpd
MNM 0.0 0.0 NA
BLIM -3.5 7.0 0.69 |elpd_diff| < 4
SLM -32.9 10.8 1.00
Diagnostic flags present.
See ?`loo-glossary` (sections `diag_diff` and `diag_elpd`)
or https://mc-stan.org/loo/reference/loo-glossary.html.According to LOO-CV, the BLIM is not much worse than the saturated multinomial model at predicting replication data. The uncertainty reflected by the standard error of this difference is high. The SLM is noticeably worse than the other two models.
In this case study, we introduced Bayesian versions of two prevalent models in knowledge structure theory, the basic local independence model and the simple learning model. We showed how to code these models in Stan, obtain posterior estimates, run diagnostic checks, and compare them using cross validation. The flexibility of the Bayesian framework makes it attractive for implementing other KST models in the future.
The pks package version 0.8-0 (Heller and Wickelmaier 2013) was used to generate artificial responses, for data handling, for providing the high school geometry data, and for plotting knowledge structures leveraging facilities from relations (Meyer and Hornik 2026) and Rgraphviz (Hansen et al. 2025). Stan was interfaced from R using rstan version 2.32.7 (Stan Development Team 2025). The loo package version 2.10.1 (Vehtari et al. 2026) was used for running LOO-CV model comparison.
Florian Wickelmaier, Department of Psychology, University of Tübingen, Schleichstr. 4, 72076 Tübingen, Germany, wickelmaier@web.de. Portions of this material were presented at the Psychoco International Workshop on Psychometric Computing, February 5-6, 2026, Padua, Italy.↩︎