I would like to submit a case study on probabilistic knowledge structures:
Knowledge structure theory is a field in psychometrics that seeks to
provide procedures for the effective diagnosis of the knowledge state of a
student in a certain domain (such as algebra, geometry, physics, or
statistics).
On the case-study page it just says “please contact us through the Stan
Forums”. Apologies if I’m missing anything. I’d be happy for any feedback.
Hi Florian –
This is indeed how to do it! The next step after you get some feedback would be opening a PR against the stan website to add it to the case study page there; here is an example of what that entails (basically just adding the built html page and then some text to link to it)
Unfortunately, it is difficult to preview the html file you’ve placed in the repository on Codeberg directly. Do you have a website you could host it on temporarily to make viewing easier?
Thanks, @Florian_Wickelmaier, and welcome to the Stan forums.
The writeup looks great.
Should there be arrows on the knowledge structure graph? If they’re implicit, would you mind saying so in the caption? I’m not convinced by this rigid knowledge structure that’s just assumed and the idea that if you know something, you’ll always get the answer right. Does this structure imply that it would be inconsistent with the model to observe someone getting question d correct and getting either a or b wrong?
I see then that the model breaks into “careless error” and “real error” and “lucky guess”. Some comparison to something like an IRT-1PL model that has a continuous difficulty parameter (as opposed to careless/careful) and IRT-3PL, which has a per-item guessing parameter.
If the graph is always subsets of the problems, then it makes sense to define the lattice as being closed under union (it’s not closed under intersection, though). I take it Q is equal to the full set of questions always? I’d be consistent and write either {} or \emptyset everywhere.
What was the motivation for pinning guess and error probabilities at max 0.5 in the BLIM? I would suggest just using <lower=0, upper=1> and letting the prior and data take the value where it should go to be consistent with the data. If you don’t have much data, the uncertainty is likely to extend beyond the 0.5 threshold. When the threshold’s in place, that mass has to go somewhere and it gets pushed down, introducing bias relative to the (0, 1)-constrained version of the model.
I wasn’t sure what “outer fringe” meant in defining `K^0$.
Figure 5 is really nice, but I’m not sure what “predicted” means here. Is there a distinction between the disks (solid circles) and circles (just border)—there’s nothing in the key for circles?
I think it’d help the reader to say a bit more about what LOO is doing so they can interpret the number dump. I’m pretty sure LOO lets you’d dump out the ELPDs and standard errors in a single output. I think that may be what you dump out at the end, but I don’t know what NA means. I think you wound up just below the 95% significance threshold on differences. This is not surprising because LOO isn’t very good at telling binary data apart. If you had more data, they’d probably be significantly different. But then I would expect the more general model to fit very well with more data, too.
What does “LOO CV is OK” mean? And what’s the NA column on the last report? How can
Estimate SE
BLIM: elpd_loo -5853.5 38.6
MNM1: elpd_loo -5861.2 38.4
produce
model elpd_diff se_diff p_worse diag_diff diag_elpd
BLIM 0.0 0.0 NA
MNM -7.6 3.7 0.98
I don’t see how the standard error on the difference can be only 3.7 when the standard errors on the ELPD estimates were 38.
I think you have a couple stray header columns (diag_diff, diag_elpd).
When you say “Again, parameter recovery is not perfect.” I think this is misleading. We don’t expect it to be perfect other than asymptotically as we get more data. If we didn’t get variation within the posterior intervals, we’d know the posterior intervals were too wide. I have the same question about circles and disks in this plot.
Any guesses as to how Lakshminarayan came up with the “knowledge graph” in figure 9? Was it empirically from test results or just a guess based on intuition? That is, if I have my own problem, where does this graph come from? With no guessing, it seems like it’s imposing a strict impossibility on some data responses.
I don’t think continuous responses are unguessable (the answer turns out to be 0 or 1 or 1/2 or sqrt(2) or pi too much of the time for it to be unguessable). I’d be inclined to fit the guessing model with a strong not going to be guessable prior.
I’m not sure what I should be comparing in figure 10. Are the shared parameters intended to be identical between models? It looks like it.
Are there frequentist packages to fit these models and if so, do you get similar results? Are there advantages to the Bayesian approach you want to include in the summary?
If whatever you’re using for layout can handle it, it’d be nice to hide most of the plotting code behind a combo-box that lets you select to show/not show the code. Same for the data munging code. I’d also hide the raw LOO computations and just present the final comparison after the first time (or even the first time as the details don’t matter so much for the point of this case study).
I’d also suggest just following the style guide in the index, which can be applied automatically. As is, there’s a lot of stylized spacing here that’s a bit distracting (like indenting priors and aligning things vertically leaving lots of space inside parens).
In writing math out, constants like function names are conventionally roman (e.g., \textrm{beta} rather than \textit{beta}.
Thanks for using reasonable amounts of precision in reporting.
Rather than a separate code supplement, I’d suggest just putting whatever you think people will need in an appendix. This isn’t to say you shouldn’t link to a GitHub or other repo, just that it’s nice if these things can be self contained. I hate coming to the end of a paper I’ve downloaded and realize I didn’t get the supplementary material.
Whenever there’s a fixed value, like here in the model block,
pi ~ dirichlet(rep_row_vector(1.0, nstate));
It’s better to move the constant up to transformed data so it doesn’t need to get reallocated each time (Stan doesn’t need the .0 like some languages—the C++ compiler figures out what each scalar type should be exactly).
transformed data {
row_vector[nstate] ones_nstate = rep_row_vector(1, nstate);
}
...
pi ~ dirichlet(ones_nstate);
Now having said that, even better is to just remove this sampling statement for pi because it’s a uniform distribution and hence a no-op as the simplex type automatically applies a uniform distribution over values matching its constraints (like all Stan types).
If you need more speed, it’d be faster to recode this:
(1 - R) * diag_pre_multiply(log( beta), K' ) +
R * diag_pre_multiply(log1m(beta), K' )
as
diag_pre_multiply((1 - R) * log(beta) + R * log1m(beta), K');
But whenever you see yourself needing a transpose, the better thing to do is just store K in transposed form so you don’t need the transpose as you go. This is the only place K is used, so just declare the transpose of K (you can call it K_t or K_tr if there would otherwise be confusion with standard naming).
I would suggest 0 as the lower bound for data size. Then you can set 0 size data to do posterior predictive checks.
I think the check eta_filter[i, j] > 0 ? 0 : 1 is easier to understand as !eta_filter[i, j] or eta_filter[i, j] == 0. You could get clever with some floating point operations to vectorize, but I don’t think the obfuscation is worth saving some vertical space.
There’s a download link, which was slow, but it worked.
Thanks Brian and Bob for your feedback!
A preview of the html version is here:
https://fwickel.codeberg.page/BayesianPKS/
I will work on Bob’s comments and report back.
Thanks again @Bob_Carpenter for your thorough review. We have updated the
material in our case study, https://codeberg.org/fwickel/BayesianPKS. The changes can be seen in the diffs. In addition, we address each of your points here.
Should there be arrows on the knowledge structure graph? If they’re
implicit, would you mind saying so in the caption? I’m not convinced by
this rigid knowledge structure that’s just assumed and the idea that if
you know something, you’ll always get the answer right. Does this
structure imply that it would be inconsistent with the model to observe
someone getting question d correct and getting either a or b wrong?
In the KST literature, the arrows usually are left out. The path between two
states represents set inclusion. We clarified this in the caption of Fig. 2.
A knowledge structure like the one displayed in Fig. 2 can be obtained in
different ways. We added a paragraph including references to further
literature.
It is true that a response pattern not included in the knowledge structure
cannot be a knowledge state and must be explained by response errors. For
example, getting only d correct (pattern 00010) is not a knowledge state.
There could be many states underlying this pattern, such as:
The pattern is not inconsistent with the probabilistic models, however, but
inconsistent with the knowledge structure. By allowing for non-zero careless
error and guessing probabilities, both the BLIM and the SLM are able to
predict every possible response pattern (where patterns that coincide with a
knowledge state usually have a higher probability than those that don’t). We
added an explanation and an example to clarify the notions of response
pattern, knowledge state, and responses errors.
I see then that the model breaks into “careless error” and “real error”
and “lucky guess”. Some comparison to something like an IRT-1PL model
that has a continuous difficulty parameter (as opposed to
careless/careful) and IRT-3PL, which has a per-item guessing parameter.
We added a paragraph that puts KST into context of other psychometric
theories: classical test theory, item response theory, and cognitive
diagnostic assessment models; the main asset of KST being that it provides a
discrete non-numerical characterization of the strengths and weaknesses of a
student within a set-theoretic framework.
If the graph is always subsets of the problems, then it makes sense to
define the lattice as being closed under union (it’s not closed under
intersection, though). I take it Q is equal to the full set of questions
always? I’d be consistent and write either {} or \emptyset everywhere.
Knowledge structures in general do not need to be closed under union nor under
intersection. As we say in the introduction, the structure displayed in Fig. 2
is special, as it is closed under union (and called a knowledge space). It is
also closed under intersection (a quasi ordinal knowledge space). Q is the
domain of the knowledge structure, the set of problems under consideration. We
consistently use the emptyset symbol now.
What was the motivation for pinning guess and error probabilities at max
0.5 in the BLIM? I would suggest just using <lower=0, upper=1> and
letting the prior and data take the value where it should go to be
consistent with the data. If you don’t have much data, the uncertainty is
likely to extend beyond the 0.5 threshold. When the threshold’s in place,
that mass has to go somewhere and it gets pushed down, introducing bias
relative to the (0, 1)-constrained version of the model.
We replaced the (0, 0.5) bounds for beta and eta by the less restrictive
beta_q + eta_q < 1. Now it is possible for betas or etas to exceed 0.5, but
not simultaneously for the same problem q. The restriction is necessary to
prevent the case beta_q = eta_q = 0.5 for all q, which implies a uniform
distribution on the response patterns, independent of the underlying knowledge
structure. Such a model would not be identifiable (see Heller, 2017).
I wasn’t sure what “outer fringe” meant in defining `K^0$.
Outer fringe refers to the definition in the introduction. We clarified this
both here and were the definition occurs.
Figure 5 is really nice, but I’m not sure what “predicted” means here. Is
there a distinction between the disks (solid circles) and circles (just
border)—there’s nothing in the key for circles?
We have added a more informative caption here and in all figures where box
plots occur. We clarified that what is displayed is the posterior predictive
distribution. The legend was misleading and has been removed.
I think it’d help the reader to say a bit more about what LOO is doing so
they can interpret the number dump. I’m pretty sure LOO lets you’d dump
out the ELPDs and standard errors in a single output. I think that may be
what you dump out at the end, but I don’t know what NA means. I think you
wound up just below the 95% significance threshold on differences. This
is not surprising because LOO isn’t very good at telling binary data
apart. If you had more data, they’d probably be significantly different.
But then I would expect the more general model to fit very well with more
data, too.
We added a bit more information about LOO-CV and better explain the output. We
agree that it’s better to leave out the number dump and to only display the
summary LOO output (including the individual ELPD values and their standard
errors).
What does “LOO CV is OK” mean? And what’s the NA column on the last
report? How canEstimate SEBLIM: elpd_loo -5853.5 38.6
MNM1: elpd_loo -5861.2 38.4
producemodel elpd_diff se_diff p_worse diag_diff diag_elpd
BLIM 0.0 0.0 NA
MNM -7.6 3.7 0.98I don’t see how the standard error on the difference can be only 3.7 when
the standard errors on the ELPD estimates were 38.
In the loo documentation (via ?loo::loo_compare), it says:
“To compute the standard error of the difference in ELPD - which should not be
expected to equal the difference of the standard errors - we use a paired
estimate to take advantage of the fact that the same set of N data points was
used to fit both models.”
So, it seems, the standard error of the difference is so low (relative to the
individual standard errors), because the ELPD values are correlated. If we do
l1 <- loo(blim1)
l2 <- loo(slm1)
l3 <- loo(mnm1)
elpd <- cbind(l1 = l1$pointwise[, "elpd_loo"],
l2 = l2$pointwise[, "elpd_loo"],
l3 = l3$pointwise[, "elpd_loo"])
cor(elpd)
we get
# l1 l2 l3
# l1 1.000 0.942 0.995
# l2 0.942 1.000 0.937
# l3 0.995 0.937 1.000
I think you have a couple stray header columns (diag_diff, diag_elpd).
We added an explanation of what these columns mean and why they are empty
here.
When you say “Again, parameter recovery is not perfect.” I think this is
misleading. We don’t expect it to be perfect other than asymptotically as
we get more data. If we didn’t get variation within the posterior
intervals, we’d know the posterior intervals were too wide. I have the
same question about circles and disks in this plot.
We rephrased this part and also changed the phrase where it appeared earlier.
We reworked the graph, so it should be clearer now.
Any guesses as to how Lakshminarayan came up with the “knowledge graph”
in figure 9? Was it empirically from test results or just a guess based
on intuition? That is, if I have my own problem, where does this graph
come from? With no guessing, it seems like it’s imposing a strict
impossibility on some data responses.
We are not sure about the origin of the knowledge structure. In her thesis,
Lakshminarayan (1995, p. 118) just says that the “knowledge structure was
constructed using the Query algorithm” without giving any details as to who
was queried or how the procedure went. Doignon and Falmagne (1999, p. 201)
state that it was “postulated, based on the content of the items”. We adopted
the second explanation as it includes the first one (Query is based on the
item content).
If you want to build your own knowledge structure, you could use one of the
procedures that we now mention in the introduction. Item tree analysis is
implemented in the pks package. No guessing does not preclude certain response
patterns from occurring: One could always be in a higher knowledge state but
make many careless errors. Such responses, however, may have a low
probability.
I don’t think continuous responses are unguessable (the answer turns out
to be 0 or 1 or 1/2 or sqrt(2) or pi too much of the time for it to be
unguessable). I’d be inclined to fit the guessing model with a strong not
going to be guessable prior.
We agree this would be an alternative strategy. We mention this now in the
text.
I’m not sure what I should be comparing in figure 10. Are the shared
parameters intended to be identical between models? It looks like it.
Yes, the results for the careless-error parameters (beta) should be similar in
both models, which is the case here. We’ve made the comparison more explicit.
Are there frequentist packages to fit these models and if so, do you get
similar results? Are there advantages to the Bayesian approach you want
to include in the summary?
We added ML estimates to Fig. 10 for comparison, the results are quite
similar. We added a paragraph about what we see as the main advantages of the
Bayesian approach to the summary.
If whatever you’re using for layout can handle it, it’d be nice to hide
most of the plotting code behind a combo-box that lets you select to
show/not show the code. Same for the data munging code. I’d also hide the
raw LOO computations and just present the final comparison after the
first time (or even the first time as the details don’t matter so much
for the point of this case study).
We put all plotting code into code boxes. Also the data preparation code,
except for when it occurs for the first time as it may help user to understand
the output below. We hid all LOO raw output and only display the summary.
I’d also suggest just following the style guide in the index, which can
be applied automatically. As is, there’s a lot of stylized spacing here
that’s a bit distracting (like indenting priors and aligning things
vertically leaving lots of space inside parens).
We applied autoformating to all Stan code (with some small adjustments).
In writing math out, constants like function names are conventionally
roman (e.g., \textrm{beta} rather than \textit{beta}.
We use \textrm{} for the distribution names now.
Thanks for using reasonable amounts of precision in reporting.
Credit goes to the packages for setting reasonable defaults.
Rather than a separate code supplement, I’d suggest just putting whatever
you think people will need in an appendix. This isn’t to say you
shouldn’t link to a GitHub or other repo, just that it’s nice if these
things can be self contained. I hate coming to the end of a paper I’ve
downloaded and realize I didn’t get the supplementary material.
We agree, that document should be self-contained as it is. All Stan code is
displayed. The Quarto source code is available via the code button at the top
of the document. We don’t think an appendix is necessary.
Whenever there’s a fixed value, like here in the model block,
pi ~ dirichlet(rep_row_vector(1.0, nstate));It’s better to move the constant up to transformed data so it doesn’t
need to get reallocated each time (Stan doesn’t need the .0 like some
languages—the C++ compiler figures out what each scalar type should be
exactly).transformed data {
row_vector[nstate] ones_nstate = rep_row_vector(1, nstate);
}
…
pi ~ dirichlet(ones_nstate);Now having said that, even better is to just remove this sampling
statement for pi because it’s a uniform distribution and hence a no-op as
the simplex type automatically applies a uniform distribution over values
matching its constraints (like all Stan types).
We omitted the uniform prior prior, but stated in the text that a Dirichlet
distribution can be used if a more informative prior is desired.
If you need more speed, it’d be faster to recode this:
(1 - R) * diag_pre_multiply(log( beta), K' ) + R * diag_pre_multiply(log1m(beta), K' )as
diag_pre_multiply((1 - R) * log(beta) + R * log1m(beta), K’);
For the example where R is a (32 x 5)-matrix, beta a 5-element vector, and K a
(9 x 5)-matrix,
R * diag_pre_multiply(log1m(beta), K')
(32 x 5) * (5 x 5) * (5 x 9)
results in a (32 x 9)-matrix. But
(1 - R) * log(beta) + R * log1m(beta)
(32 x 5) * (5 x 1) + (32 x 5) * (5 x 1)
is a 32-element vector, and a (32 x 32) diagonal matrix is incompatible with
K’ (5 x 9). Please let us know if we misunderstood your comment.
But whenever you see yourself needing a transpose, the better thing to do
is just store K in transposed form so you don’t need the transpose as you
go. This is the only place K is used, so just declare the transpose of K
(you can call it K_t or K_tr if there would otherwise be confusion with
standard naming).
We define K_t = K' and K1m_t = (1 - K)' in the transformed data block.
I would suggest 0 as the lower bound for data size. Then you can set 0
size data to do posterior predictive checks.
We were pondering on this comment, but didn’t know what to do, we’re sorry. Do
you mean prior predictive checks (not posterior)? But even then, we
couldn’t figure out how to implement these checks without commenting out the
likelihood. Which component(s) would you set to zero? (Feel free to ignore if
it’s not so important.)
I think the check eta_filter[i, j] > 0 ? 0 : 1 is easier to understand as
!eta_filter[i, j] or eta_filter[i, j] == 0. You could get clever with
some floating point operations to vectorize, but I don’t think the
obfuscation is worth saving some vertical space.
We use eta_filter[i, j] == 0 now.