Statistical Methods
Reference

Classical Test Theory

CTT characterizes items and tests using observed-score statistics from an \(N\times K\) response matrix (\(N\) respondents, \(K\) items). Items may be dichotomous (0/1) or polytomous (non-negative integers). All variance estimates use the sample \((n-1)\) denominator throughout. Missing data are handled pairwise (default) or listwise.

Item-Level Descriptive Statistics

Mean, SD, and Standard Error

Let \(\mathcal{V}_j\) be the set of respondents with non-missing responses to item \(j\) and \(n_j=|\mathcal{V}_j|\). Then

Relative Difficulty (p-value)

Items with \(p_j<0.15\) or \(p_j>0.85\) are flagged. For dichotomous items this equals the classical p-value.

Discrimination Indices

Adjusted Point-Biserial (Item-Rest Correlation)

Item \(j\) is correlated with the sum of the remaining \(K-1\) items, eliminating the inflation from item \(j\) appearing in its own total:

This is Pearson's \(r\), computed pairwise over non-missing observations. Items with \(r_{pb,j}<0.20\) are flagged.

Biserial and Polyserial Correlation

These indices estimate the latent linear association assuming the true variable is continuous and normal but observed only in coarsened categories. For an item with \(M\) categories, let \(\hat{F}_m\) be the cumulative proportion through category \(m\). The internal threshold is \(z_m=\Phi^{-1}(1-\hat{F}_m)\) and its standard normal ordinate is \(\varphi(z_m)\):

The last category is excluded because its threshold is \(+\infty\). For \(M=2\) this reduces to the classical biserial formula. \(\Phi^{-1}\) is computed via the Beasley-Springer-Moro rational approximation (accuracy \(\approx10^{-9}\)).

Reliability and Standard Error of Measurement

Cronbach's \(\alpha\)

where \(s_j^2\) is the sample variance of item \(j\) and \(s_T^2\) is the sample variance of total scores (Cronbach, 1951). Both use the \((n-1)\) denominator.

Alpha-if-Item-Deleted

For each item \(j\), \(\alpha\) is recomputed on the \((K-1)\)-item submatrix formed by removing column \(j\). An increase signals that item \(j\) is reducing internal consistency.

Standard Error of Measurement

Item Response Models

PACER supports binary and polytomous items in any mixture within a single test. The scaling constant \(D\) defaults to 1.7 (close to the normal-ogive metric; \(D=1\) gives the pure logistic metric). All probabilities are clamped to \([\varepsilon,\,1-\varepsilon]\) with \(\varepsilon\approx1.49\times10^{-8}\) to prevent \(\log 0\) in the likelihood.

Binary Item Response Models

Three-Parameter Logistic (3PL)

2PL and 1PL / Rasch

Setting \(c=0\) gives the 2PL. The 1PL additionally fixes \(a=1\) and estimates a free population standard deviation \(\sigma\) that scales the quadrature nodes: \(\theta_k^*=\sigma z_k\), where \(\{z_k\}\) are the standard GH nodes.

Graded Response Model (GRM)

The GRM (Samejima, 1969) models items with \(K+1\) ordered categories \((0,1,\ldots,K)\). The boundary response function gives the cumulative probability of responding in category \(k\) or higher:

with ordered thresholds \(d_1

Boundary probabilities use the sign convention \(p_{kq}=\sigma(Da(d_k-\theta_q))\), so larger \(d_k\) means a harder threshold. Category probabilities are formed by differencing consecutive rows of the \(K\times Q\) boundary matrix.

Generalized Partial Credit Model (GPCM) and PCM

The GPCM (Muraki, 1992) uses adjacent-category transitions. For an item with \(K+1\) categories \((0,\ldots,K)\), define the log-numerator of category \(k\):

with \(\eta_0\equiv0\) and a leading \(d_0=0\) prepended to the step-parameter vector. Category probabilities follow the multinomial logit (softmax):

Evaluated via log-sum-exp: \(P(X=k)=\exp(\eta_k-M)/\sum_j\exp(\eta_j-M)\) where \(M=\max_j\eta_j\), preventing overflow at extreme \(\theta\). The PCM (Masters, 1982) is the GPCM with \(a\equiv1\).

Rating Scale Model (RSM)

The RSM (Andrich, 1978) constrains the PCM so that all items share identical step-difficulty offsets \(\tau_1,\ldots,\tau_K\) (\(\tau_0=0\) prepended by convention). Each item has a unique location \(\delta_i\), and its step parameters are

Response probabilities are evaluated with the PCM formula using these compound parameters. The free parameter vector is \([\delta_1,\ldots,\delta_I,\,\tau_1,\ldots,\tau_K]\).

Item Calibration: Marginal Maximum Likelihood via EM

PACER estimates item parameters by marginal maximum likelihood (MML), integrating the latent trait \(\theta\) out of the likelihood. Because \(\theta\) is unknown, direct maximization of the item parameters is not feasible; instead PACER forms the marginal likelihood and maximizes it using the Bock-Aitkin (1981) EM algorithm. Unique response patterns are identified before the E-step to eliminate redundant likelihood evaluations.

The Marginal Likelihood

The goal is to estimate the full vector of item parameters \(\bm{\gamma}\). Since the person parameters \(\theta_i\) are also unknown, they are integrated out under an assumed population distribution \(g(\theta\mid\mu,\sigma)\):

The integral has no closed form and is evaluated by Gauss-Hermite quadrature over \(Q\) nodes:

In practice, logs of individual pattern likelihoods are accumulated rather than computing products, ensuring numerical stability. When the data contain both binary and polytomous items (\(J\) binary, \(M\) polytomous), the pattern likelihood factors as a product over both item types:

where \(P_j\) is the binary ICC (1PL/2PL/3PL) and \(P_l\) is the polytomous category probability (GPCM or GRM). Both item types are estimated within the same EM loop.

Model Identification

For the 2PL and 3PL, the population distribution is fixed at \(N(0,1)\) to identify the latent scale. For the 1PL (Rasch model), only the mean need be fixed; PACER therefore estimates \(N(0,\sigma)\), with \(\sigma\) as a free parameter interpretable as the common discrimination across all items.

Starting Values

Because the EM algorithm is iterative, starting values are required. Three strategies are available. Random draws each parameter from a uniform distribution over the admissible parameter space — useful when the likelihood may be multimodal; running the analysis multiple times with different random starts is a recommended diagnostic. Fixed sets \(a_j=1\), \(b_j=0\), \(c_j=0.25\) for all binary items and uses ordered random starting values for polytomous step parameters. Logistic uses logistic regression intercepts to initialize the difficulty parameters (and proportional-odds intercepts for polytomous step parameters), with \(a\) and \(c\) initialized as in the fixed option — this is the default and generally performs best.

E-Step: Posterior Weights and Sufficient Statistics

Given current item parameters, the likelihood of response pattern \(\bm{x}_i\) at node \(\theta_k\) is

The posterior weight is

The marginal log-likelihood used to monitor convergence is

Binary Sufficient Statistics

For each binary item \(j\) and node \(k\), summed over respondents \(\mathcal{C}_j\) who answered item \(j\):

\(r_{jk}\) is the expected number of correct responses at node \(k\); \(f_{jk}\) is the expected total count.

Polytomous Sufficient Statistics

For polytomous item \(j\), the sufficient statistic at node \(k\) for category \(s\) is

In the pattern-based implementation, the posterior weight for unique pattern \(p\) at node \(k\) is \(r_{p,k}=c_p\cdot w_k\cdot L_k/m_p\) where \(c_p\) is the pattern count and \(m_p=\sum_k L_k w_k\).

M-Step: Binary Items (2PL / 1PL / 3PL)

The M-step maximizes the expected complete-data log-likelihood for item \(j\):

PACER computes analytical gradients for the 2PL. The gradients with respect to \(a_j\) and \(b_j\) are:

In the multigroup path, the gradient accumulates across all groups \(g\) using group-scaled nodes \(\theta_{gk}=\mu_g+\sigma_g z_k\):

For the 3PL, \(c_j\) is parameterized on the logit scale internally; its gradient is computed by forward differences.

M-Step: GPCM and PCM

The complete-data log-likelihood for GPCM item \(j\) is the multinomial log-likelihood in expected category counts:

The analytical gradient with respect to \(a_j\) uses the conditional mean of the log-numerator:

For step parameter \(d_j\) (\(j=1,\ldots,K\)), let \(P_k^{\ge j}=\sum_{s\ge j}P(X=s\mid\theta_k)\) be the exceedance probability. Then:

This is a score-residual: model-expected exceedance minus the observed indicator. For the PCM, \(a_j\equiv1\) and the \(a\)-gradient is not computed.

M-Step: GRM

The complete-data log-likelihood for GRM item \(j\) is:

Let \(f_k^{(m)}=P^*(X\ge m\mid\theta_k)\cdot[1-P^*(X\ge m\mid\theta_k)]\) be the logistic slope of the \(m\)-th boundary. The gradient with respect to \(a_j\) depends on which category was observed:

The term \(\Delta_k(s)\) depends on the observed category \(s\):

The gradient with respect to threshold \(d_m\) is nonzero only for the two categories adjacent to it:

GRM item gradients are parallelized over items at the M-step.

M-Step: RSM

The RSM free parameter vector is \(\bm{\xi}=[\delta_1,\ldots,\delta_I,\,\tau_1,\ldots,\tau_K]\). Within the optimizer, item \(i\)'s step vector is \(\bm{d}_i=[0,\,\delta_i+\tau_1,\ldots,\delta_i+\tau_K]\) and probabilities are evaluated with the PCM formula (\(a=1\)). The gradients follow from the GPCM step gradient with \(a=1\):

The \(\tau_t\) gradient is the sum of per-item step gradients across all \(I\) items, since \(\partial d_{it}/\partial\tau_t=1\) for all \(i\). Item gradients are computed in parallel; \(\tau\) gradients are accumulated sequentially to avoid data races.

Parameter Priors

Priors are incorporated by augmenting the marginal log-likelihood with log-prior terms, one per item parameter:

Priors on \(a\) and \(b\) are optional. A prior on \(c\) is applied by default when estimating the 3PL and can be disabled for very large samples.

Beta Prior for the Guessing Parameter

The beta distribution is parameterized by shape parameters \(\alpha\) and \(\beta\), but it is more natural to specify a prior mean and standard deviation. PACER converts these to shape parameters via the method of moments. Starting from

solving for \(\alpha\) and \(\beta\) gives

The defaults \(\mu_c=0.20\) and \(\sigma_c=0.087\) correspond to \(\mathcal{B}(4,16)\), which places most mass below 0.35 and is appropriate for items with three to five response options. As an alternative, a logit-normal prior can be selected, which places the prior on \(\text{logit}(c)\) and may be preferable when the guessing parameter is expected to be very small.

Multigroup Calibration (Bock-Zimowski)

For data from multiple groups, PACER implements the Bock-Zimowski (1997) multigroup MML-EM. Each group \(g\) has distribution \(\theta\sim N(\mu_g,\sigma_g^2)\); the reference group (\(g=0\)) is fixed at \(N(0,1)\) for identification.

Group-Scaled Nodes

Weights \(\{w_k\}\) are invariant because the Jacobian and density ratio cancel exactly.

Multigroup E-Step

For respondent \(i\) in group \(g\), probabilities are evaluated at \(\theta_{gk}\):

Sufficient statistics \(r_g[j,k]\) and \(f_g[j,k]\) are accumulated per group on the same quadrature index \(k\).

Group Distribution Update (Closed Form)

Let \(W_{gk}=\sum_{i\in g}h_{ik}\). After the item M-step:

The outer loop iterates until

where \(\tau_\ell=\max(\text{Tol},\,10^{-5})\). The convergence condition on the group parameters requires that both the mean and standard deviation shift by less than \(10^{-4}\) for every non-reference group.

Optimizer: L-BFGS-B

Item parameters are estimated with L-BFGS-B, a limited-memory quasi-Newton method with box constraints. The inverse Hessian approximation is built from the most recent \(m=10\) curvature pairs \(\{(\bm{s}_t,\bm{y}_t)\}\) where \(\bm{s}_t=\bm{x}_{t+1}-\bm{x}_t\) and \(\bm{y}_t=\nabla f_{t+1}-\nabla f_t\). The two-loop recursion computes the search direction \(\bm{p}_t=-H_t^{-1}\nabla f_t\) without explicitly forming \(H_t^{-1}\). Box constraints are enforced by projecting onto the active constraint set. Convergence: \(|\Delta f|/(0.1+|f|)<\tau\) (default \(\tau=10^{-8}\)). Convergence is evaluated relative to the current function value to ensure scale-invariant termination.

Ability Estimation

Given calibrated item parameters, PACER estimates each respondent's latent trait \(\hat\theta\) from their observed response vector. Four methods are available. All operate on sparse response vectors so that respondents with missing items are scored without imputation.

Maximum Likelihood Estimation (MLE)

A dense grid from \(-8\) to \(+8\) in steps of 0.05 locates the global basin; gradient descent with backtracking then refines from the grid minimum. Extreme scores: perfect and null scores produce unbounded estimates. PACER subtracts (or adds) 0.5 from the response of the item with the smallest discrimination \(a_j\). Standard error from observed information via numerical second differentiation:

Maximum A Posteriori (MAP)

The normal prior \(N(\mu,\sigma^2)\) (defaults: \(\mu=0\), \(\sigma=1\)) is added to the log-likelihood and the same minimizer is used. MAP is always finite, including for extreme scores.

Expected A Posteriori (EAP)

EAP is the posterior mean under a normal prior, approximated on a uniform quadrature grid spanning \(\mu\pm4\sigma\) in \(Q\) steps with weights proportional to \(\phi((\theta_k-\mu)/\sigma)\) (normalized to sum to one):

The posterior standard deviation serves as the standard error:

TCC Scoring

TCC scoring (Lord, 1980) inverts the Test Characteristic Curve \(\tau(\theta)=\sum_j E[X_j\mid\theta]\) to find

solved on \([-8,8]\) by Brent's root-finding method. The SE is \(1/\sqrt{I(\hat\theta_{\text{TCC}})}\) from the test information function.

Boundary Condition

A unique TCC estimate exists only when the raw score exceeds the sum of the guessing parameters across all binary items:

When this condition is not met — that is, when the observed score is no higher than what could be achieved by guessing alone — TCC assigns the same lower-bound \(\hat\theta\) to all such examinees. Similarly, a perfect raw score produces the upper-bound value defined by the search interval. The default interval \([-8,8]\) covers the range used in practice; if scores fall outside this range the interval endpoints can be widened.

Post-Calibration EAP Thetas

After calibration, PACER computes a \(\hat\theta\) for every respondent using 11-point GH quadrature with a diffuse prior \(N(0,10)\) (\(\sigma=\sqrt{10}\approx3.162\)). These thetas are used for DIF matching and item-fit analysis, not score reporting.

Test Characteristic Curve, Information, and SEM

Test Characteristic Curve (TCC)

When proportional TCC is enabled, output is divided by \(X_{\max}=\sum_j\max_k k\) to give a 0–1 scale.

Binary Item Information

For 1PL and 2PL:

For 3PL:

GRM Item Information

Samejima (1969):

where \(f_k(\theta)=P^*(X\ge k\mid\theta)[1-P^*(X\ge k\mid\theta)]\) is the logistic slope of the \(k\)-th boundary, with sentinels \(f_0\equiv0\) and \(f_{K+1}\equiv0\). Concretely, consecutive logistic-slope values are differenced and each squared difference is divided by the corresponding category probability, then summed and scaled by \(D^2a_j^2\).

GPCM / PCM / RSM Item Information

Information equals \(D^2a^2\) times the conditional variance of the item score (Muraki, 1992):

The test-level SEM at each \(\theta\) is

IRT Equating and Linking

IRT linking places parameter estimates from a new test form onto the scale of an existing item bank using common anchor items. PACER supports six methods in three families. Mixed binary and polytomous items are supported; GRM uses the cumulative-logistic expected score; GPCM/PCM/RSM use the adjacent-category expected score.

The Linear Scale Transformation

All methods estimate \((A,B)\) mapping the new-form theta scale to the bank scale: \(\theta_{\text{bank}}=A\theta_{\text{form}}+B\). Item parameters transform as

The reverse transformation is \(a_{\text{rev}}=Aa\), \(b_{\text{rev}}=(b-B)/A\), \(d_{\text{rev},k}=(d_k-B)/A\).

Stocking-Lord Methods

Stocking and Lord (1983) minimize the weighted integral of squared TCC differences over \(Q\) GH nodes. Let \(\tau_1(\theta)\) be the bank TCC and \(\tau_2^*(A,B;\theta)\) the linked new-form TCC.

Haebara Methods

Haebara (1980) replaces TCC comparison with item-by-item ICC comparison. For binary items \(P_{jk}\) is the 3PL ICC; for polytomous items it is \(E[X_j\mid\theta_k]\).

All four characteristic-curve objectives are minimized by two-dimensional Nelder-Mead simplex (starting simplex side 0.5, up to 2,000 iterations, tolerance \(10^{-10}\)).

Mean/Sigma and Mean/Mean

Mean/Sigma (MS)

The location pool: \(b\)-parameters for binary items; all \(d\)-parameters pooled as a flat vector for GPCM/PCM (or averaged per item when the single-location option is selected); GRM boundary parameters; item location \(\delta_i\) only for RSM items.

Mean/Mean (MM)

dropD² Iterative Anchor Purification

After the initial equating, a per-item discrepancy statistic is computed:

For polytomous items \(P_j\) is the expected score \(E[X_j\mid\theta_k]\). Values are standardized to

The item with the largest \(|\tilde{D}^2_j|\) exceeding the cutoff (default 2.5) is removed; equating is re-run on the reduced anchor set. The process repeats until no item exceeds the cutoff or fewer than two anchors remain.

Cognitive Diagnostic Models

CDMs classify examinees into mastery profiles over \(K\) binary latent attributes \(\bm{\alpha}_i=(\alpha_{i1},\ldots,\alpha_{iK})\in\{0,1\}^K\). The Q-matrix \(\bm{Q}\in\{0,1\}^{J\times K}\) specifies which attributes are required by each item.

DINA Model

The Deterministic Inputs, Noisy And-gate (DINA) model sets the ideal response indicator as \(\xi_{ij}=\prod_{k=1}^{K}\alpha_{ik}^{q_{jk}}\). A respondent achieves \(\xi_{ij}=1\) only if they have mastered all required attributes. The response probability accommodates slipping (\(s_j\)) and guessing (\(g_j\)):

DINO Model

The Deterministic Inputs, Noisy Or-gate (DINO) model requires mastery of at least one required attribute:

EM Estimation of Item Parameters

Item parameters (\(s_j\), \(g_j\)) and class probabilities (\(\pi_c\)) are estimated jointly by marginal maximum likelihood over the \(L=2^K\) latent classes \(\bm{\alpha}^{(c)}\) with prior class probabilities \(\pi_c\ge0\), \(\sum_c\pi_c=1\). The algorithm is initialized with equal class probabilities (\(\pi_c=1/L\)) and uniform starting values for \(g_j\) and \(s_j\). The population distribution is re-estimated at each step, so the posteriors are empirical Bayes estimates.

E-Step

Let \(\Pr(\bm{\alpha}_l\mid\bm{x}_i)\) be the posterior probability that examinee \(i\) has skill profile \(\bm{\alpha}_l\). The expected number of examinees with profile \(\bm{\alpha}_l\) is

and the expected number of those examinees who answered item \(j\) correctly is

Let \(I^0_{jl}\) be the expected count of examinees who lack at least one attribute required for item \(j\) (i.e., \(\xi_{ijl}=0\)) and \(R^0_{jl}\) be the subset of those who answered correctly. Analogously, \(I^1_{jl}\) and \(R^1_{jl}\) correspond to examinees who possess all required attributes (\(\xi_{ijl}=1\)).

M-Step

With these expected counts, the guessing and slip parameters are updated in closed form:

Class probabilities are updated as \(\hat\pi_c=N^{-1}\sum_i\tau_{ic}\). The EM loop alternates between E-step and M-step until the change in parameters across iterations falls below the convergence tolerance.

Skill Profile Estimation

Once item parameters are calibrated, each examinee's skill profile is estimated from the likelihood of their response pattern \(\bm{x}_i\) given each of the \(L=2^K\) possible profiles:

Three estimators are available.

Maximum Likelihood (MLE)

The MLE selects the profile that maximizes the likelihood across all \(L\) possible classifications:

Because skill profiles are discrete, this is computed by exhaustive evaluation over all \(L\) profiles. The MLE may be multimodal — two or more profiles can share the same maximum likelihood. PACER flags such cases and reports all profiles achieving the maximum.

Maximum A Posteriori (MAP)

The MAP incorporates a prior distribution \(\Pr(\bm{\alpha}_l)\) over skill profiles (supplied by the user), maximizing the posterior:

where \(\Pr(\bm{\alpha}_l)\) is a discrete prior such that \(\sum_l\Pr(\bm{\alpha}_l)=1\). Like the MLE, the MAP can be multimodal.

Expected A Posteriori (EAP)

The EAP uses the full posterior distribution over profiles. The posterior for profile \(\bm{\alpha}_l\) given \(\bm{x}_i\) is

Let \(\bm{\tau}_i = [\Pr(\bm{\alpha}_1\mid\bm{x}_i),\ldots,\Pr(\bm{\alpha}_L\mid\bm{x}_i)]\) be the \(1\times L\) posterior vector and \(\bm{P}\) the \(L\times K\) binary matrix of all possible skill profiles. The EAP estimate is

This is a \(1\times K\) vector where each element is the marginal posterior probability of mastery for the corresponding attribute — equivalent to \(\sum_{l:\,\alpha_k^{(l)}=1}\Pr(\bm{\alpha}_l\mid\bm{x}_i)\). Unlike the MLE and MAP, the EAP yields continuous mastery probabilities between 0 and 1 and is always unique.

Examinee Skill Profile Estimation

Once item parameters are calibrated, each examinee's mastery profile is characterized by their full posterior distribution over the \(2^K\) latent classes. Three quantities are reported.

Modal Profile (MAP Classification)

The modal-class estimate assigns the examinee to the latent class with the highest posterior probability:

The estimated attribute profile is \(\hat{\bm{\alpha}}_i = \bm{\alpha}^{(\hat{c}_i)}\), a binary vector of length \(K\) indicating mastery (1) or non-mastery (0) of each attribute. This is the primary classification output.

Marginal Attribute Mastery Probabilities

The marginal probability of mastering attribute \(k\) is obtained by summing the posterior over all latent classes in which attribute \(k\) is mastered:

These \(K\) probabilities provide a continuous measure of mastery confidence. A value near 1 indicates strong evidence of mastery; near 0 indicates strong evidence of non-mastery; intermediate values indicate classification uncertainty for that attribute.

Classification Certainty

The posterior entropy of the class distribution for examinee \(i\) quantifies overall classification certainty:

Low entropy indicates that most posterior weight is concentrated on a single class (high certainty). High entropy indicates spread across multiple classes, flagging examinees whose profiles should be interpreted cautiously. The minimum possible entropy is 0 (perfect certainty); the maximum is \(\ln(2^K)=K\ln 2\) (uniform uncertainty).

Differential Item Functioning

DIF occurs when respondents from different groups who have the same ability have systematically different probabilities of answering an item correctly. PACER implements the Mantel-Haenszel procedure (Holland & Thayer, 1988) for binary items and the Liu-Agresti (1996) generalization for polytomous items. All procedures first stratify respondents into \(K\) equal-frequency bins (default \(K=10\)) using a matching variable — either total score or IRT theta. Binning uses linear interpolation on the empirical quantile function so that equal proportions of respondents fall in each stratum; the lowest observed score is included in the first bin.

Mantel-Haenszel Procedure (Binary Items)

Within each stratum \(s\), a \(2\times2\) contingency table is formed:

The MH log-odds ratio pools across strata:

The ETS D-DIF statistic:

Negative D-DIF means the reference group has higher odds of answering correctly at the same ability — the item favors the reference group. Note: some published implementations invert this sign; PACER applies the convention stated above.

Chi-square with Yates continuity correction (1 df):

The p-value is computed from the regularized incomplete gamma function via series and continued-fraction expansions (accurate to \(\approx10^{-10}\)).

Generalized MH Procedure (Polytomous Items)

For polytomous items with score values \(v_1<\cdots

The test statistic \(T^2/\mathrm{Var}(T)\) is chi-squared with one degree of freedom.

Standardized Mean Difference (SMD)

Positive SMD means the focal group scores higher (item favors focal); negative means the item favors the reference group. Some published implementations use \(|\text{SMD}|\) for the direction determination; PACER uses the raw signed value as defined above. The standardized SMD is

ETS DIF Classification

Gauss-Hermite Quadrature

All integrals over the latent-trait distribution are approximated by Gauss-Hermite quadrature with \(Q\) nodes (default \(Q=21\)):

where \(\phi\) is the standard normal density. Nodes \(\{z_k\}\) are roots of the probabilists' Hermite polynomial \(\text{He}_Q(x)\); weights satisfy \(w_k=v_{k,1}^2\sqrt{2\pi}\) where \(v_{k,1}\) is the first eigenvector component. For common values of \(Q\in\{11,15,21,31,41,51\}\), PACER uses a pre-computed node-and-weight lookup table for speed. For arbitrary \(Q\), nodes and weights are computed on demand via the symmetric QL algorithm applied to the tridiagonal Jacobi matrix with off-diagonal elements \(\beta_i=\sqrt{i}\). For multigroup calibration, nodes are shifted and scaled: \(\theta_{gk}=\mu_g+\sigma_g z_k\). Weights are invariant because the Jacobian and density ratio cancel exactly.

Optimization Routines

L-BFGS-B (Item Calibration)

Item parameters are estimated by L-BFGS-B. The two-loop recursion computes the quasi-Newton direction \(\bm{p}_t=-H_t^{-1}\nabla f_t\) from \(m=10\) curvature pairs without forming \(H_t^{-1}\) explicitly. Box constraints are enforced by projecting onto the active constraint set. Step length satisfies the strong Wolfe conditions via backtracking. Convergence: \(|\Delta f|/(0.1+|f|)<\tau\) (default \(\tau=10^{-8}\)).

Gradient Descent with Backtracking (Theta Estimation)

MLE and MAP theta estimation use gradient descent with backtracking. At each step: \(\theta_{\text{new}}=\theta-\alpha[-\ell'(\theta)]\) where \(\ell'(\theta)\approx[\ell(\theta+h)-\ell(\theta-h)]/(2h)\), \(h=10^{-6}\). Step size \(\alpha\) is increased by 1.2 on progress and halved on rejection; convergence when \(\alpha<10^{-10}\).

Brent's Method (TCC Inversion)

TCC scoring inverts \(\tau(\theta)=X_{\text{obs}}\) on \([-8,8]\). Brent's method combines bisection (guaranteed progress), secant (linear interpolation), and inverse quadratic interpolation (parabolic fit through three points). The method falls back to bisection whenever the other two strategies would not make sufficient progress, giving guaranteed \(O(\log\varepsilon^{-1})\) worst-case convergence. Tolerance: \(10^{-10}\).

Nelder-Mead Simplex (Equating Objectives)

The two-parameter SL and HB equating objectives are minimized over \((A,B)\) by Nelder-Mead simplex. The simplex of three vertices is updated by reflection, expansion, contraction, or shrinkage at each iteration. Initial simplex side: 0.5; maximum iterations: 2,000; convergence tolerance on objective range: \(10^{-10}\).

References

Andrich, D. (1978). A rating formulation for ordered response categories. Psychometrika, 43, 561–573.

Bock, R. D., & Aitkin, M. (1981). Marginal maximum likelihood estimation of item parameters: Application of an EM algorithm. Psychometrika, 46, 443–459.

Bock, R. D., & Zimowski, M. F. (1997). Multiple group IRT. In W. J. van der Linden & R. K. Hambleton (Eds.), Handbook of modern item response theory. Springer.

Cronbach, L. J. (1951). Coefficient alpha and the internal structure of tests. Psychometrika, 16, 297–334.

Doran, H. (2023). A collection of numerical recipes useful for building scalable psychometric applications. Journal of Educational and Behavioral Statistics. https://doi.org/10.3102/10769986221116905

Haebara, T. (1980). Equating logistic ability scales by a weighted least squares method. Japanese Psychological Research, 22, 144–149.

Holland, P. W., & Thayer, D. T. (1988). Differential item performance and the Mantel-Haenszel procedure. In H. Wainer & H. I. Braun (Eds.), Test validity. Lawrence Erlbaum.

Liu, I. M., & Agresti, A. (1996). Mantel-Haenszel-type inference for cumulative odds ratios with a stratified ordinal response. Biometrics, 52, 1223–1234.

Lord, F. M. (1980). Applications of item response theory to practical testing problems. Erlbaum.

Masters, G. N. (1982). A Rasch model for partial credit scoring. Psychometrika, 47, 149–174.

Muraki, E. (1992). A generalized partial credit model: Application of an EM algorithm. Applied Psychological Measurement, 16, 159–176.

Samejima, F. (1969). Estimation of latent ability using a response pattern of graded scores. Psychometrika Monograph Supplement, No. 17.

Stocking, M. L., & Lord, F. M. (1983). Developing a common metric in item response theory. Applied Psychological Measurement, 7, 201–210.