Throughout this section we will make use of matrices (i.e., two dimensional arrays of parameters) to simplify our explanations. We will denote such matrices by bold, italicized capital letters—for example, A for a matrix containing all slope parameters—and denote their individual elements by italicized lower case letters as we have done thus far, with row indices indicated by subscripts and column indices indicated by parenthesized superscripts. Bold lower case letters will indicate row or column vectors.

In the case of matrices with dimensions of size IC (corresponding to the set of all response categories on all items), we use the index pair (i,c) as a shorthand for the index corresponding to category c of item i. The definitions and shapes of several key matrices are provided in Table 1 for reference.

Specifying a procedure to uniquely identify the MNCM is required because the model has a number of symmetries that result in invariances—mathematical transformations to the model parameters that leave the predicted probabilities unchanged. Consequently, varying a set of free parameters to find the best fit to a data set does not uniquely determine those parameters, which are arbitrary within variations that honor the invariances. In order to compare results from different data sets—or even from different subsets of the same data—we must impose additional constraints on the parameters in order to identify (uniquely specify) the model. A familiar example from dichotomous 2PL IRT is that shifting item difficulties and student abilities by the same amount does not change the predicted probability for student s to answer item i correctly; typically, the corresponding identification constraint for this invariance is to set the mean ability to zero.

The greater mathematical complexity of the MNCM results in a greater number of symmetries, and hence a larger number of identification constraints are necessary to specify reproducible results. Additionally, for a multidimensional model in which the probabilities depend on a scalar product (like the MNCM; see Equation 4), rotations of the vector space leave the probabilities unchanged; we exploit several affordances of this when fitting data (both synthetic data and actual student data).

Our identification procedure below is a specification of certain constraints on model parameters (or on properties of sets of parameters such as means or standard deviations), and the parameter transformations we use to effect these constraints. Since by definition the identification procedures do not change the probabilities (and often do not even change the tendencies), they can be applied either while fitting the model to data as an integral part of that process, or afterward as a separate post-processing step on the estimated parameters.

Note that in order to leave the predicted probabilities unchanged, enforcing identification constraints on some parameters often requires making corresponding changes to other parameters linked by the model invariances.

As mentioned earlier, we take a Bayesian approach to fitting the MNCM. In the Bayesian paradigm, the probability of a student selecting a particular response category is still modeled by Equations 1, 4; however, every model parameter [each individual ad(i,c), b^(i,c), and θs(d)] is treated as having an entire probability distribution over possible values rather than a single “optimal” value. This provides a principled way of modeling the effects of parameter uncertainty. In addition, Bayesian methods allow us to incorporate reasonable prior beliefs about how parameters will be distributed. For example, it is typical to assume a priori that population values will be more-or-less normally distributed. These two qualities make Bayesian modeling well suited to the challenges of fitting the real data described earlier.

Our choice of prior probability distribution is based on the work of Natesan et al. (2016), who recommended the use of hierarchical models for Bayesian IRT:

In a hierarchical model, some of the prior probability distributions are not fully specified but are instead parameterized by additional “higher level” random variables [here, α_d and β, which serve as scale parameters for the priors of ad(i,c) and b^(i,c) respectively]. These variables are then learned from the data. Such an approach allows item parameters with high-confidence estimates to inform the scale of the priors, which in turn better stabilizes estimates for the parameters which have less supporting data.

Note that all prior means for both ability and item parameters are fixed at zero. Due to the invariances of the model, this does not result in any loss of generality. Instead, it serves as a “soft” identification constraint which stabilizes the location of the parameters during the fitting process. Similarly, the scale of Θ (and therefore of A) is identified by the fixed variance of the θs(d) prior, and the rotation and shearing of the model are partly identified by the use of independent priors for each dimension (which implicitly results in diagonal covariance matrices for both A and Θ, again without any loss of generality). The signs and ordering of dimensions are left unidentified until after the fitting process, though this does not seem to adversely affect convergence due to our choice of fitting method discussed below.

The process of fitting a Bayesian model to data is called inference, and its output is an updated joint probability distribution over all model parameters called the posterior. Given a set of observed responses R, a prior probability distribution p(Ƶ) over parameters Ƶ≡{Θ,A,b,α,β}, and a response model p(R∣Ƶ) describing the probabilities of the observations given the parameters (Equations 1, 3 for the MNCM), then the posterior p(Ƶ∣R) may be found by applying Bayes' rule:

This equation, however, is deceptive in its simplicity; except in special cases, it is not possible to find a closed-form analytical expression for the posterior, and it must be approximated numerically.

Our work uses a variational inference (VI) approach in order to approximately solve Equation 7. VI methods work by re-framing Bayesian inference as an optimization problem: given an approximate (but analytically tractable) parameterized probability distribution, a numerical optimizer searches for the parameters which bring this approximant closest to the true posterior. Surprisingly, it is possible to do this without ever computing the true posterior by instead maximizing a surrogate objective function known as the Evidence Lower Bound (ELBO); further details are widely available in the academic literature (e.g. Blei et al., 2017). In addition to being considerably faster than the more typical Markov Chain Monte Carlo methods used for Bayesian inference, this optimization-based approach allows VI to converge to just one of the many (arbitrary) permutation of dimensions and signs in an under-identified model.

While the form of the approximate posterior distribution in VI may be arbitrarily complex, a popular and quite effective simplification is to treat each random variable in the model as having its own independent univariate posterior distribution. This is known as the mean-field approximation. For all real-valued model variables (θs(d), ad(i,c), and b^(i,c)), we choose as the approximate posterior a normal distribution parameterized by a mean and a standard deviation, both of which may be freely varied by the optimizer.

The approximate posterior for the higher-level random variables α_d and β require additional care since they cannot take negative values. This constraint is handled by introducing surrogate real-valued variables—with posteriors modeled by unconstrained normal distributions as above—and mapping these to positive numbers using an appropriate bijection (e.g., the softplus function x ↦ log[1 + exp(x)]).

When selecting the dimensionality of a multidimensional IRT model, it is important to balance improvements in predictive ability from additional dimensions against the added degrees of freedom thereby introduced. Not doing so inevitably leads one to select a model which performs well only on the specific data used to estimate its parameters, but fails to explain new data generated from the same underlying statistic process or yields psychologically meaningless parameter estimates. This problem is known as overfactoring—a type of overfitting—and is often dealt with in a classical maximum-likelihood context by using likelihood-ratio tests or information criteria to compare models with different dimensionalities (van der Linden, 2016, chs. 17 & 18).

Bayesian methods are not entirely immune to overfactoring, though their principled inclusion of parameter uncertainty does provide some built-in protection against it. Several approaches to Bayesian dimensionality assessment were recently compared by Revuelta and Ximénez (2017) for non-hierarchical MNCM models with up to three dimensions estimated using Markov Chain Monte Carlo methods. The authors recommended the use of a standardized generalized dimensionality discrepancy measure (SGDDM) in this context, noting that alternatives based on discrepancy measures struggled to correctly identifying the dimensionality of synthetic data.

In our work, we find that using a hierarchical model in combination with variational inference performs simultaneous parameter estimation and dimensionality selection: even when a D-dimensional model is specified, some of those may be “turned off” during inference (by setting the corresponding α_d and a_d to ~0) when the observed data provide insufficient evidence to confidently estimate their slopes. This results in considerable robustness to overfactoring as illustrated by our simulation study in Section 5.

The method described above was implemented in the Python programming language (version 3.11). We leveraged the NumPyro probabilistic programming framework (Phan et al., 2019) to define the probabilistic model, automatically generate a corresponding mean-field normal approximating posterior (including surrogate variables and bijections for α_d and β), and perform approximate inference using NumPyro's built-in Stochastic Variational Inference (SVI) algorithm.

As the name might suggest, SVI relies in part on random sampling to evaluate and optimize the variational objective function, and must therefore be paired with a noise-tolerant optimization routine. We use the Adam optimizer (Kingma and Ba, 2014) configured with a variable “learning rate” parameter which is programmed to decrease from 0.05 to zero on a predetermined 30,000-step schedule.³

While SVI will often converge even if the posterior means of all variables are randomly initialized, we attempt to provide a more reasonable starting point for optimization by computing initial guesses of these using a fast IRT approximation (Zhang et al., 2020),⁴ resorting to random initialization only when this approach does not yield a solution. The posterior standard deviations of all variables are initialized to a fixed value of 0.1, which is the default value for mean-field normal posteriors in NumPyro.

The raw, post-optimization outputs of the SVI algorithm consist of an estimated mean and standard deviation for each variable in the approximate mean-field posterior—that is, for each individual ad(i,c), b^(i,c), θs(d), α_d, and β. However, these raw outputs are only weakly identified by the Bayesian prior during inference, and require post-processing to exactly impose the identification constraints described in Section 4.2. We perform this step analytically, using the posterior standard deviations to account for uncertainty when finding the covariance matrices of A and Θ (which are needed to fully identify the scale and rotation of the model) and the variance of b (which is included for completeness). The standard deviations are not used further in this work; we examine only the means of the identified parameters, which correspond to an expected a posteriori (EAP) solution.

Lastly, any dimensions inactivated during inference (Section 4.5) are omitted from the final estimates. These dimensions may be easily detected by computing the magnitudes of the estimated slope vectors and applying a simple threshold criterion; we use ||a_d|| ≤ 0.005.

The final identified outputs of our method therefore include:

Our synthesized datasets each comprised a set of synthetic student and item parameters Θ, A, and b (all identified according to the constraints in Section 4.2) and a corresponding synthetic response matrix R. All datasets used I = 30 items and C = 5 categories per item, matching the structure of the real data we will analyze in Section 6. The standard deviation of the b vector was fixed at 1.5, which is also consistent with that found later for the real dataset; this yielded synthetic data with a range of observed responses fractions across the categories in each item.

For the multidimensional simulation study, datasets were generated with D = 9 synthesized dimensions. Six sample size levels S ∈ {50, 100, 200, 600, 2, 000, 10, 000} and two A-matrix covariance structures (described below) were explored according to a fully-crossed design. For each condition, 100 replications were generated with different pseudo-random parameter values for each replication, yielding a total of 1,200 datasets.

The invariances of the MNCM afford us substantial flexibility in choosing a covariance structure for the synthesized A matrices. Of note, all possible covariance matrices for Θ and A—even those with arbitrary structure and multiple highly-correlated dimensions—are expressible as diagonal matrices in some choice of reference frame due to these invariances,⁵ which allows us to synthesize A to have uncorrelated dimensions without any loss of generality. The standard deviations of these dimensions were selected to span a fairly wide gamut in order to explore the limits of parameter recovery using our method; we fix the standard deviation of the first dimension to 1.0 and specify that each subsequent dimension is smaller than the previous one by a factor γ ∈ {0.8, 0.512}:

The first condition, γ = 0.8, yields data in which the smallest dimension accounts for approximately 1% of the overall variance in the synthesized tendencies. The second condition yields data with a much more rapid decrease in variance (as might be expected in a real dataset with fewer significant factors), with the specific choice of γ = 0.8³ = 0.512 intended to facilitate comparison across the two conditions: the standard deviations of dimensions 1–3 in the γ = 0.512 data exactly match those of dimensions 1, 4, and 7 in the γ = 0.8 data.

For the unidimensional simulation study comparing the Bayesian method to an established software package implementing Marginal Maximum Likelihood Estimation (MMLE), datasets were generated with D = 1, sample sizes S ∈ {50, 100, 200, 600, 2, 000, 10, 000}, and a standard deviation of 1.0 for the sole slope vector a₁. Again, 100 replications were generated per condition, yielding an additional 600 datasets for this second study.

Our interest in finding misconceptions means that, when studying real data, we will need to rotate the coordinate system of our results to associate each dimension with an interpretable concept (common practice in exploratory analysis methods; see Section 6.1). We are therefore interested in evaluating how well the method recovers the basis-independent information present in our synthesized parameters and especially in the A matrix, rather than the particular coordinate system in which it initially extracts this information (which is arbitrary and does not affect the response probabilities due to the rotational invariance property of the MNCM). We achieving this by aligning the coordinate systems of the synthesized and recovered parameters prior to computing any evaluation metrics.

Note that using the same identification criteria for both sets of parameters does somewhat succeed in aligning their coordinate systems. Even so, the identified coordinate directions may themselves be sensitive to small errors in the parameters; this procedure may therefore lead to inflated errors which are due more to (arbitrary) differences in rotation than to differences which actually affect the response probabilities.

In this work, we use an orthogonal Procrustes procedure (Schönemann, 1966) to align the recovered and synthesized parameters while respecting the invariances of the model. We first find an orthogonal matrix Q which best maps Θ to Θ^ in the least-squares sense, i.e., which minimizes

where ||·||F2 is the squared Frobenius norm (equal to the sum of squares of all matrix elements). This matrix is then used to rotate the synthesized parameters to allow direct comparisons with the estimates:

Note that such rotations do not alter the predicted probabilities in Equations 1, 4.

The choice to align parameters to the reference frame of the estimates—as opposed to that of the synthesized parameters—is a deliberate one. Because of the self-limiting dimensionality of the MNCM-Bayes method, the number of dimensions extracted may be smaller than the number of dimensions synthesized. Remaining in the reference frame of the estimates allows us to limit our analysis to only these D^≤D extracted dimensions, and also allows us to examine parameter recovery metrics for each extracted dimension individually. In the reference frame of the synthesized parameters however, there is no longer a one-to-one correspondence between these extracted dimensions and individual coordinate directions, and limiting the analysis to the D^-dimensional subspace of extracted dimensions becomes difficult or impossible.

For this paper, we focus exclusively on the recovery of the slope parameters ad(i,c) as these provide the most information about the structure and content of an assessment instrument (which is our current research focus). We leave exploration of Θ and b for future work.

We evaluate parameter recovery separately for each dimension. Our metric of choice is the squared Pearson correlation coefficient—also known as the coefficient of determination—computed between the synthesized, Procrustes-aligned ad⋆(i,c) values and their estimated counterparts a^d(i,c):

We intentionally use a correlation-based metric instead of the more common root mean squared error (RMSE) in order to forgive differences in overall scale, which do not affect the subjective interpretation of the recovered parameters.

The r² metric ranges from 0 to 1 and may be understood as the fraction of variation in the slope estimates attributable to variation in the (rotated) ground-truth values. This definition implies that values of r² are analogous to reliability coefficients, except that these apply to recovered slope rather than test scores. We therefore suggest similar norms be used when evaluating r²: values greater than 0.9 suggest sufficient accuracy for interpreting individual slopes, while those as low as 0.7 may still provide value when interpreting multiple slopes in aggregate.

A key challenge that arises in applying multi-dimensional IRT models like the MNCM to real data is that of ascribing meaning to the dimensions thereby discovered. As with principal component analysis, the dimensions extracted by MNCM-Bayes are initially identified in a principal basis. Such solutions are generally not easy to interpret, as each dimension loads on many categories across many items.

The standard approach for increasing the interpretability of such results (both in multidimensional IRT and in classical factor analysis techniques) is to find a transformation which increases some measure of “simplicity” while not changing the modeled probabilities. Many such measures exist, as do standard methods for transforming solutions to maximize them given a matrix of slopes or factor loadings. For the sake of brevity, we will forgo a principled evaluation of these methods for now and present only one as an example: a bi-factor rotation method proposed by Jennrich and Bentler (2011) and implemented in the GPArotation R package (Bernaards and Jennrich, 2005).

Bi-factor rotation methods such as the one above transform the coordinate system of Θ and A such that each distractor has non-negligible slopes along only a small number of dimensions (ideally just two). The first of these dimensions is always a ‘general' factor on which all distractors load—we have found the ability in this dimension to be strongly associated with students' overall test score (Spearman rank correlation coefficient ρ_s = 0.98). The remainder are “group” factors associated with only a small set of distractors. We call the slope vectors for each of these sparse distractor vectors, and define their positive directions as those which yield a positive Spearman correlation between the corresponding rows of Θ and the test scores. (This sign determination method is sometimes marginal as some dimensions do not correlate strongly with score; for those that do, however, this causes the slope components which most characterize the non-Newtonian nature of each distractor to have negative signs).

To demonstrate the results of this transformation, we plot in Figure 3 the first two principal and sparse distractor vectors found for the large post-instruction data set described above. We display the components of each vector as a pattern of dots on an I × C grid, each dot having a size and intensity proportional to the magnitude of that component and being colored red when negative (blue if positive).

Each sparse distractor vector loads heavily on just a few distractors. To determine whether a particular vector represents a misconception, we examine the most heavily weighted distractors (typically the largest ~half-dozen assuming these stand out prominently from the background when plotted as in Figure 3) and see if selecting them would indicate consistent application of some alternate hypothesis to Newtonian mechanics. This process is admittedly quite subjective and could likely be improved upon in future work (e.g., by using outlier analysis to differentiate between prominent and background values), but nevertheless identifies several clear examples of misconceptions in our present results.

The limited results we just discussed are clear-cut examples showing that the combination of the MNCM-Bayes method followed by rotations to find sparse distractor vectors can discover known misconceptions—indeed, the two just discussed are among the top three found by Wheatley et al. (2022) using modified module analysis (Wells et al., 2019). Our results also illustrate potential improvements to our understanding of existing misconceptions, for example by showing that the impetus concept applies more strongly for motion in a curved path than in a straight path, or that the belief that only the last force applied determines an object's motion is also associated with a peculiar view of what happens when such forces stop acting.

The remainder of our results in this example application are included in the supplemental materials for this paper (both as figures and as tabulated slope coefficients). While several additional dimensions in these results seem to have precedents in prior misconceptions research, we postpone further discussion to a later application-focused paper with a more thorough analysis of our full FCI dataset (~34,000 exams, including pre- and post-instruction data from eight colleges and universities).

Many opportunities exist for improving the MNCM-Bayes method or extending its capabilities. One that we are already investigating is the choice of factor rotation method, which must ideally balance ease of interpretation with consistency of the discovered misconceptions. By comparing results obtained using different rotations across data from several colleges, we hope to better inform this choice in future version of our method. Another modification would be to allow partially-specified patterns of loadings. This would allow us to manually associate some dimensions with particular known misconceptions (by constraining any irrelevant slopes to zero) and could aid in steering any unconstrained dimensions toward novel misconceptions instead of existing ones.

On the applications side, a compelling first question for future study is whether the various misconceptions that we found for the FCI are consistent across different colleges, and whether similar misconceptions are as robustly found in pre-instruction data rather than post-instruction data; we plan to address both topics in a forthcoming paper (in preparation). One could also study how student ability scores for each misconception depends on factors such as gender, preparation, or background. And, by comparing results from pre- and post-instruction data, one might determine which instructional approaches are most effective in reducing the persistence of various misconceptions.

As our methods are applied to other research-designed instruments, it seems likely that new misconceptions—or re-contextualized versions of existing ones—will be discovered, especially in subject areas where few studies of student misconceptions exist. Finally, we need not limit ourselves to data from traditional concept tests: modified versions of this method could be applied to entire online courses where frequently-given wrong answers have already been mined from student data (as is done by e.g., MasteringPhysics.com and ExpertTA.com) and could be treated as distractors. Further applications await.