fit.MGPCM: Fit the Multidimensional Generalized Partial Credit Model...

View source: R/fit.MGPCM.R

fit.MGPCMR Documentation

Fit the Multidimensional Generalized Partial Credit Model (MGPCM)

Description

Fits a multidimensional generalization of the partial credit model to polytomous (multi-category) response data. Two estimation backends are provided: full Bayesian inference via Hamiltonian Monte Carlo (Stan) and a fast stochastic-EM algorithm (iStEM).

Usage

fit.MGPCM(
  data,
  D = NULL,
  Q.matrix = NULL,
  length.poly = NULL,
  method = c("iStEM", "stan"),
  control.model = NULL,
  control.method = NULL
)

Arguments

data

An N \times I matrix of integer responses coded 0, 1, \dots, K_i - 1. Rows index persons, columns index items.

D

Integer; number of latent dimensions (D \ge 1). Default is 2.

Q.matrix

An optional I \times D binary matrix. Entry q_{id} = 1 frees a_{id}; q_{id} = 0 fixes it to 0. Default uses a triangular identification structure.

length.poly

Optional integer scalar or length-I vector giving the number of categories per item. If NULL, counts are inferred from observed maxima; supply this argument when a valid category is unobserved in the sample.

method

Estimation method: "iStEM" (default) or "stan".

control.model

A named list of model-level hyperparameters. Supported entries:

a.mu, a.sigma

Prior location and scale for \log a_{id} (log-normal). Defaults: 0.25, 0.25. Controls the prior mean and spread of the discrimination parameters across dimensions.

d.mu, d.sigma

Prior mean and SD for free category intercepts d_{ik} (k \ge 1; normal). Defaults: 0, 1. The first category intercept d_{i0}=0 is fixed for identification.

theta.mu

Prior mean vector for \boldsymbol{\theta}_j. Default: rep(0, D). A vector of length D specifying the prior mean for each latent dimension.

L

Theta grid size per dimension for marginal log-likelihood computation and iStEM block Gibbs sampling. Default adapts to D (e.g., 61 for D = 1, 31 for D = 2, 15 for D = 3). Larger grids increase numerical precision at the cost of exponential growth in computation (L^D total nodes).

theta.lower, theta.upper

Bounds for the theta grid used by marginal log-likelihood computation and iStEM block Gibbs sampling. Defaults: -6, 6.

control.method

A named list of method-specific tuning parameters. Common entries (used by both Stan and iStEM):

cores

Number of CPU cores for parallel chains (Stan) or ignored (iStEM). Default: the number of chains.

vis

Logical; if TRUE (default), prints progress information to the console.

seed

Random seed for reproducibility. Default: a random integer.

Stan-specific entries:

chains

Number of MCMC chains (default: 2).

iter

Total iterations per chain (default: 5000).

warmup

Warmup/burn-in iterations per chain (default: iter / 2).

thin

Thinning interval (default: 1).

init

Initial values: "random" (default) for uniform(-2, 2) initialization, or a list of initial values per chain.

algorithm

MCMC algorithm: "HMC" (default), "HMC", or "Fixed_param".

adapt_delta

Target average acceptance probability (NUTS; default: 0.95). Values closer to 1 reduce step size and improve sampling for difficult posteriors.

max_treedepth

Maximum tree depth (NUTS; default: 10). Increase if "max treedepth exceeded" warnings appear.

stepsize

Initial step size for the leapfrog integrator (auto-tuned by Stan if not set).

int_time

Total integration time for HMC trajectories (only when algorithm = "HMC").

metric

Mass matrix type: "unit_e", "diag_e" (default), or "dense_e".

adapt_engaged

Logical; if TRUE (default), warmup adaptation is enabled.

adapt_init_buffer, adapt_term_buffer, adapt_window

Warmup adaptation scheduling parameters (defaults: 25, 50, 25).

iStEM-specific entries:

M

Number of burn-in batches retained for Geweke convergence diagnosis (default: 10; must be \ge 2).

B

Batch size: MCMC iterations per batch (default: 20).

burnin.maxitr

Maximum burn-in batches (default: 100).

maxitr

Maximum total batches (default: 2000).

eps1

Geweke z-score convergence threshold (default: 1.5).

eps2

Monte Carlo error tolerance (default: 0.4).

frac1, frac2

Fractions for the Geweke diagnostic (defaults: 0.1, 0.5).

corr.optim.maxit

Maximum L-BFGS-B iterations for the constrained unit-diagonal correlation update (default: 50).

optim.maxit

Maximum L-BFGS-B iterations per item (default: 50).

fix.corr

Logical; fix correlations to identity (default: FALSE).

estimate.se

Logical; compute standard errors from final MC chain (default: TRUE).

a.lower, a.upper

Bounds on a_{id} (defaults: 1e-4, 6).

d.lower, d.upper

Bounds on free category intercepts d_{ik} for k \ge 1. Defaults are -8 and 8 in the iStEM backend.

Value

An object of class "MGPCM" with components:

npar

Number of free parameters.

method

"stan" or "iStEM".

theta

List with est, se, Rhat (N \times D).

par

List with est, se, Rhat, free (I \times (D + K_{max})). Columns are a1..aD, d0, d1, ..., d_{K_{max}-1}.

Corr

List with est, se, Rhat (D \times D).

length.poly

Integer vector of per-item category counts.

logLik

Marginal log-likelihood (class "logLik").

call, arguments

Call and argument records.

Model Specification

Let Y_{ij} \in \{0, 1, \dots, K_i - 1\} denote the categorical response of person j = 1, \dots, N to item i = 1, \dots, I, where K_i \ge 2 is the number of response categories for item i. Let \boldsymbol{\theta}_j be the D-dimensional latent trait vector.

The category response probability implemented in both the C++ and Stan backends is a softmax over category scores:

P(Y_{ij} = k \mid \boldsymbol{\theta}_j) = \frac{ \exp\{k\,\eta_{ij} + d_{ik}\} }{ \sum_{r=0}^{K_i-1} \exp\{r\,\eta_{ij} + d_{ir}\} }, \qquad k = 0, 1, \dots, K_i - 1,

where

\eta_{ij} = \sum_{d=1}^{D} a_{id}\,\theta_{jd}.

The discrimination parameters satisfy a_{id} > 0 when q_{id} = 1 and a_{id} = 0 when q_{id} = 0. The d_{ik} values are category intercepts, not cumulative step difficulties; d_{i0} = 0 is fixed for identification and d_{i1}, \dots, d_{i,K_i-1} are free.

This formulation nests the standard (unidimensional) generalized partial credit model (Muraki, 1992) when D = 1 and all q_{i1} = 1.

Prior distributions:

a_{id} (free)

Log-normal: \log a_{id} \sim N(\mu_a, \sigma_a^2)

d_{ik} (free)

Normal: d_{ik} \sim N(\mu_d, \sigma_d^2) for k \ge 1; d_{i0}=0 is fixed. The same normal penalty is used by Stan and by the iStEM item update when use.prior = TRUE.

\boldsymbol{\theta}_j

Multivariate normal: \boldsymbol{\theta}_j \sim N_D(\mathbf{0}, \boldsymbol{\Sigma})

Estimation Methods

Stan (method = "stan"):

Full Bayesian inference via HMC. The joint posterior is sampled with multiple chains, providing posterior means, standard deviations, and \hat{R} diagnostics.

iStEM (method = "iStEM"):

Improved Stochastic EM alternating finite-grid block Gibbs person-sampling and item-parameter optimization (L-BFGS-B). Convergence monitored via Geweke diagnostics and batch-means MC error.

References

Muraki, E. (1992). A generalized partial credit model: Application of an EM algorithm. Applied Psychological Measurement, 16(2), 159–176. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1177/014662169201600206")}

Yao, L., & Schwarz, R. D. (2006). A multidimensional partial credit model for polytomous data. Applied Psychological Measurement, 30(4), 295–318.

See Also

sim.data.MGPCM, get.fit.index.MGPCM, logLik.MGPCM, rotate

Examples

sim <- sim.data.MGPCM(N = 20, I = 6, D = 2, length.poly = 4)
fit <- fit.MGPCM(sim$response, D = 2, method = "iStEM",
                 control.method = list(
                   vis = FALSE, seed = 123,
                   M = 2, B = 2, burnin.maxitr = 2,
                   maxitr = 3, eps1 = 10, eps2 = 10,
                   estimate.se = FALSE))
head(fit$theta$est)
fit$par$est[1:5, ]
gof <- get.fit.index(fit)
summary(gof)


ForceChoice documentation built on Sept. 13, 2026, 1:06 a.m.