Package 'ForceChoice'

ForceChoice: Forced-Choice Modeling Based on IRT and CDM

Description

The ForceChoice package provides a unified framework for fitting, simulating, and evaluating forced-choice and traditional item response models. It supports eight model families spanning dominance (MIRT), ideal-point/unfolding (MGGUM), polytomous (MGPCM), and forced-choice (FCMIRT, FCGGUM, TIRT, FCDCM, FCGDINA) paradigms, with full Bayesian estimation via Stan, fast improved stochastic EM (iStEM), and a deterministic EM backend for FCGDINA.

Model Families

Traditional (single-stimulus) models:

fit.MIRT

Multidimensional IRT (1PL–4PL). Binary responses with optional Q-matrix structure.

fit.MGPCM

Multidimensional Generalized Partial Credit Model. Polytomous ordered-category responses.

fit.MGGUM

Multidimensional Generalized Graded Unfolding Model. Ideal-point polytomous responses.

Forced-choice (comparative) models:

fit.FCMIRT

Forced-choice MIRT. Item-level dominance model with sequential ranking over item endorsement logits at the block level.

fit.FCGGUM

Forced-choice GGUM. Item-level ideal-point model with sequential ranking over binary GGUM endorsement logits.

fit.TIRT

Thurstonian IRT for forced-choice. Pairwise probit comparison of latent utility differences.

fit.FCDCM

Forced-choice DCM. Higher-order cognitive diagnostic model with exact attribute-profile marginalization.

fit.FCGDINA

Forced-choice GDINA. General CDM item-response structure with forced-choice block likelihood.

Estimation Methods

Model families support a common estimation interface:

Stan (method = "stan"):

Full Bayesian inference via Hamiltonian Monte Carlo (NUTS/HMC). Provides posterior means, standard deviations, and convergence diagnostics ($\hat{R}$).

iStEM (method = "iStEM"):

Improved Stochastic EM (iStEM) algorithm combining finite-grid block Gibbs person sampling with L-BFGS-B item-parameter optimization. Scales well to moderate data sets.

EM (method = "EM"):

Deterministic posterior-weight EM for FCGDINA. Useful as a reproducible baseline or fast diagnostic estimator.

Core Workflow

A typical analysis follows four steps:

1. Simulate (or load) data with sim.data.*().

2. Fit a model with fit.*().

3. Rotate the solution if needed with the rotate S3 generic.

4. Evaluate fit with get.fit.index() and summary.good.of.fit().

Goodness-of-Fit

The package uses the limited-information M2 framework (Maydeu-Olivares & Joe, 2005, 2006) implemented in good.of.fit. Fit indices include:

• Information criteria: AIC, AICc, BIC, CAIC, SABIC, HQIC

• Absolute fit: M2, RMSEA (with CI), SRMSR, McDonald NCI

• Comparative fit: CFI, TLI, IFI

• Pseudo-$R^2$: McFadden, Cox–Snell, Nagelkerke, etc.

• Local dependence: Yen's Q3 (raw and adjusted)

• Classification: posterior entropy, mean max posterior

Reproducibility

The package is designed so that simulation, estimation, extraction, and model checking can be scripted end to end. Simulation functions return the true parameters and the effective call arguments; fitted objects retain the original input arguments, estimates, convergence diagnostics, and likelihood-related quantities used by logLik and get.fit.index. For computationally expensive Stan or iStEM workflows, manuscript replication code should set seeds explicitly and use a small demonstration configuration in addition to any full-scale analysis.

Author(s)

Maintainer: Haijiang Qin [email protected] (ORCID) [copyright holder]

Authors:

• Haijiang Qin [email protected] (ORCID) [copyright holder]

• Lei Guo [email protected] (ORCID) [copyright holder]

Examples

set.seed(123)
sim <- sim.data.FCGDINA(N.person = 12, N.block = 2, I.block = 2,
                        D = 2, model = "DINA")
str(sim$response)
head(sim$Q.matrix)

---

S3 Methods: coef

Description

Extract estimated model coefficients (item, statement, or block parameters) from fitted ForceChoice model objects. This provides a standard R interface for retrieving parameter estimates, standard errors, and convergence diagnostics in a consistent format across all eight model families.

Usage

## S3 method for class 'MIRT'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'MGPCM'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'MGGUM'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'FCMIRT'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'FCDCM'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'FCGDINA'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'FCGGUM'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

## S3 method for class 'TIRT'
coef(object, type = c("par", "theta", "Corr", "all"), ...)

Arguments

object	An object of one of the following classes: "MIRT" — Multidimensional IRT model. "MGPCM" — Multidimensional Generalized Partial Credit Model. "MGGUM" — Multidimensional Generalized Graded Unfolding Model. "FCMIRT" — Forced-Choice Multidimensional IRT model. "FCDCM" — Forced-Choice Diagnostic Classification Model. "FCGDINA" — Forced-Choice GDINA model. "FCGGUM" — Forced-Choice Generalized Graded Unfolding Model. "TIRT" — Thurstonian IRT for Forced-Choice.
type	Character string specifying which parameters to extract: "par" (default) for item/statement parameters, "theta" for person trait estimates, "Corr" for factor correlations, or "all" for a list with all three.
...	Additional arguments (currently ignored).

Details

For FCDCM objects, type = "par" returns a list with components delta (higher-order parameters, $D \times 2$) and par (block $\eta$ parameters, $B \times 2$), reflecting the hierarchical structure of the model.

Value

type = "par"

A matrix of item/statement/block parameter estimates. Rows index items/statements/blocks and columns index parameter types. Column names indicate the parameter (e.g., "a1", "b", "d0", "tau1", "eta0").

type = "theta"

An $N \times D$ matrix of person parameter (latent trait) estimates.

type = "Corr"

A $D \times D$ inter-trait correlation matrix.

type = "all"

A named list with components par, theta, and Corr.

Methods (by class)

• coef(MIRT): Extract coefficients from MIRT objects. Returns $I \times (D+3)$ matrix with columns a1..aD, b, c, d.

• coef(MGPCM): Extract coefficients from MGPCM objects. Returns $I \times (D+K_{\max})$ matrix with columns a1..aD, d0..d(K-1).

• coef(MGGUM): Extract coefficients from MGGUM objects. Returns matrix with columns a1..aD, delta1..deltaD, tau0..tau(K-1).

• coef(FCMIRT): Extract coefficients from FCMIRT objects. Returns statement parameter matrix with columns a1..aD, b, c, d.

• coef(FCDCM): Extract coefficients from FCDCM objects. Returns a list with delta ($D \times 2$ higher-order parameters) and par ($B \times 2$ block $\eta$ parameters). type = "theta" returns the higher-order trait, and type = "Corr" returns a scalar 1.

• coef(FCGDINA): Extract coefficients from FCGDINA objects. type = "par" returns CDM delta parameters, type = "theta" returns posterior attribute mastery probabilities, type = "Corr" returns NA (no continuous latent trait).

• coef(FCGGUM): Extract coefficients from FCGGUM objects. Returns matrix with columns a1..aD, delta1..deltaD, tau0..tau(K-1).

• coef(TIRT): Extract coefficients from TIRT objects. Returns $I \times (D+1)$ matrix with columns lambda1..lambdaD, psi2.

Examples

sim <- sim.data.MIRT(N = 20, I = 6, D = 2, model = "m2pl")
fit <- fit.MIRT(
  sim$response, model = "m2pl", 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)
)

coef(fit)                  # item parameter estimates
coef(fit, type = "theta")  # person trait estimates
coef(fit, type = "Corr")   # factor correlation matrix
coef(fit, type = "all")    # everything in a list

---

S3 Methods: confint

Description

Computes Wald-type confidence intervals for the model parameters of a fitted ForceChoice model object. Intervals are constructed as $\hat{\psi}_k \pm z_{1-\alpha/2}\,\mathrm{SE}_k$, where $\hat{\psi}_k$ is the point estimate and $\mathrm{SE}_k$ is the standard error of the $k$-th free parameter.

Usage

## S3 method for class 'MIRT'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'MGPCM'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'MGGUM'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'FCMIRT'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'FCDCM'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'FCGDINA'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'FCGGUM'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'TIRT'
confint(object, parm = NULL, level = 0.95, ...)

Arguments

object	A fitted model object.
parm	Optional numeric indices or character names selecting free parameters. By default, intervals are returned for all free parameters.
level	Numeric; confidence level (default: 0.95 for 95% intervals).
...	Additional arguments (currently ignored).

Value

A numeric matrix with two columns, 2.5% and 97.5% (column names adjust automatically for non-default level). Each row corresponds to one free parameter.

Methods (by class)

• confint(MIRT): Confidence intervals for MIRT parameters.

• confint(MGPCM): Confidence intervals for MGPCM parameters.

• confint(MGGUM): Confidence intervals for MGGUM parameters.

• confint(FCMIRT): Confidence intervals for FCMIRT parameters.

• confint(FCDCM): Confidence intervals for FCDCM parameters.

• confint(FCGDINA): Confidence intervals for FCGDINA parameters.

• confint(FCGGUM): Confidence intervals for FCGGUM parameters.

• confint(TIRT): Confidence intervals for TIRT parameters.

---

S3 Methods: deviance

Description

Returns the model deviance, defined as $-2 \ell$ where $\ell$ is the marginal log-likelihood evaluated at the parameter estimates.

Usage

## S3 method for class 'MIRT'
deviance(object, ...)

## S3 method for class 'MGPCM'
deviance(object, ...)

## S3 method for class 'MGGUM'
deviance(object, ...)

## S3 method for class 'FCMIRT'
deviance(object, ...)

## S3 method for class 'FCDCM'
deviance(object, ...)

## S3 method for class 'FCGDINA'
deviance(object, ...)

## S3 method for class 'FCGGUM'
deviance(object, ...)

## S3 method for class 'TIRT'
deviance(object, ...)

Arguments

object	A fitted model object.
...	Additional arguments (currently ignored).

Value

A numeric scalar: $-2 \log L$.

Methods (by class)

• deviance(MIRT): Model deviance ($-2 \ell$) for MIRT objects.

• deviance(MGPCM): Model deviance for MGPCM objects.

• deviance(MGGUM): Model deviance for MGGUM objects.

• deviance(FCMIRT): Model deviance for FCMIRT objects.

• deviance(FCDCM): Model deviance for FCDCM objects.

• deviance(FCGDINA): Model deviance for FCGDINA objects.

• deviance(FCGGUM): Model deviance for FCGGUM objects.

• deviance(TIRT): Model deviance for TIRT objects.

---

S3 Methods: extract

Description

A generic S3 extractor function designed to retrieve internal components from fitted model objects produced by the ForceChoice package. This function provides a consistent interface across all eight model classes, allowing users to access estimated parameters, fit statistics, latent trait estimates, convergence diagnostics, and model configuration data without directly manipulating the internal list structure.

Usage

## S3 method for class 'MIRT'
extract(object, what, ...)

## S3 method for class 'MGPCM'
extract(object, what, ...)

## S3 method for class 'MGGUM'
extract(object, what, ...)

## S3 method for class 'FCMIRT'
extract(object, what, ...)

## S3 method for class 'FCDCM'
extract(object, what, ...)

## S3 method for class 'FCGDINA'
extract(object, what, ...)

## S3 method for class 'FCGGUM'
extract(object, what, ...)

## S3 method for class 'TIRT'
extract(object, what, ...)

Arguments

object	An object of one of the following classes: "MIRT" — Multidimensional IRT model. "MGPCM" — Multidimensional Generalized Partial Credit Model. "MGGUM" — Multidimensional Generalized Graded Unfolding Model. "FCMIRT" — Forced-Choice Multidimensional IRT model. "FCDCM" — Forced-Choice Diagnostic Classification Model. "FCGDINA" — Forced-Choice GDINA model. "FCGGUM" — Forced-Choice Generalized Graded Unfolding Model. "TIRT" — Thurstonian IRT for Forced-Choice.
what	A character string specifying the name of the component to extract. Valid choices depend on the class of object. See Details section for full listings.
...	Additional arguments passed to methods (currently ignored).

Details

This function supports extraction from the eight ForceChoice model classes. Below are the available components for each:

MIRT, MGPCM, MGGUM

Traditional (single-stimulus) models. Available components:

par

List (est, se, Rhat, free) of item parameter arrays.

theta

List (est, se, Rhat) of $N \times D$ person parameter matrices.

Corr

List (est, se, Rhat) of $D \times D$ inter-trait correlation matrices.

logLik

Marginal log-likelihood (class "logLik").

npar

Number of free parameters.

method

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

Q.matrix

The $I \times D$ Q-matrix.

length.poly

Vector of category counts per item (MGPCM, MGGUM only).

stan.obj

The stanfit object (Stan only).

MCMC.obj

The MCMC extract list (Stan only).

iStEM

List of iStEM diagnostics (iStEM only).

call

The matched call.

arguments

List of arguments used for fitting.

FCMIRT, FCGGUM

Forced-choice models. In addition to the traditional-model components listed above, the following are available:

block.items

List of item indices in each block.

response

The $N \times B$ forced-choice response matrix.

patterns

List of permissible ranking patterns per block.

patterns.total

List of all possible ranking patterns per block.

fc.type

Character vector of forced-choice formats ("RANK", "MOLE", "PICK").

FCDCM

Forced-Choice DCM. Available components:

par

List (est, se, Rhat, free) of block $\eta$ parameter matrices ($B \times 2$).

delta

List (est, se, Rhat) of higher-order $\delta$ parameters ($D \times 2$).

theta

List (est, se, Rhat) of $N \times 1$ higher-order trait estimates.

alpha

Posterior mean attribute profile ($N \times D$).

alpha.patterns

Full enumeration of $2^D$ attribute mastery patterns.

zeta.patterns

Condensation outputs for each pattern.

dcm.type

DCM condensation rule ("DINA" or "DINO").

response

The $N \times B$ binary block response matrix.

block.items

List of two-statement blocks.

patterns

List of binary response patterns per block.

logLik

Marginal log-likelihood (class "logLik").

npar

Number of free parameters.

method

Estimation method.

Q.matrix

The $I \times D$ Q-matrix.

stan.obj, MCMC.obj, iStEM, EM

Estimation backend objects.

call

The matched call.

arguments

List of arguments used for fitting.

FCGDINA

Forced-Choice GDINA. Available components:

delta

List (est, se, Rhat) of item-level CDM delta parameters.

alpha

Posterior attribute summaries, including posterior class probabilities when available.

alpha.patterns

Full enumeration of $2^D$ attribute mastery patterns.

design.matrix.list

Per-item CDM design matrices.

response, block.items, patterns, patterns.total, fc.type

Forced-choice data structures.

TIRT

Thurstonian IRT. Available components:

par

List (est, se, Rhat, free) of statement-level parameter matrices ($I \times (D+1)$).

theta

List (est, se, Rhat) of $N \times D$ person parameter matrices.

gamma

List (est, se, Rhat, free) of pairwise $\gamma$ parameters.

gamma.matrix

List (est, se, Rhat) of $I \times I$ skew-symmetric $\gamma$ matrices.

Corr

List (est, se, Rhat) of $D \times D$ inter-trait correlation matrices.

pairs.matrix

Matrix of pairwise comparison indices.

pairs.value

Person-specific pair data (MOLE/PICK).

response

The $N \times I_{pairs}$ pairwise response matrix.

block.items

List of item indices in each block.

fc.type

Character vector of forced-choice formats.

Q.matrix

The statement-level Q-matrix.

logLik

Marginal log-likelihood (class "logLik").

npar

Number of free parameters.

method

Estimation method.

stan.obj, MCMC.obj, iStEM

Estimation backend objects.

call

The matched call.

arguments

List of arguments used for fitting.

Value

The requested component. Return type varies depending on what and the class of object. If an invalid what is provided, an informative error is thrown listing valid options.

Methods (by class)

• extract(MIRT): Extract components from MIRT objects.

• extract(MGPCM): Extract components from MGPCM objects.

• extract(MGGUM): Extract components from MGGUM objects.

• extract(FCMIRT): Extract components from FCMIRT objects.

• extract(FCDCM): Extract components from FCDCM objects.

• extract(FCGDINA): Extract components from FCGDINA objects.

• extract(FCGGUM): Extract components from FCGGUM objects.

• extract(TIRT): Extract components from TIRT objects.

Examples

sim <- sim.data.MIRT(N = 20, I = 6, D = 2, model = "m2pl")
fit <- fit.MIRT(
  sim$response, model = "m2pl", 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)
)

extract(fit, "par")       # item parameter estimates
extract(fit, "theta")     # person trait estimates
extract(fit, "Corr")      # factor correlation matrix
extract(fit, "npar")      # number of free parameters
extract(fit, "logLik")    # marginal log-likelihood
extract(fit, "iStEM")     # iStEM convergence diagnostics

---

Fit the Forced-Choice Diagnostic Classification Model (FCDCM)

Description

Fits a higher-order cognitive diagnostic model for forced-choice (paired- comparison) data. Each block presents two statements, and the respondent selects the one they agree with more. The discrete attribute mastery profile $\boldsymbol{\alpha}_j$ for each person is marginalized over all $2^D$ possible patterns, linking a continuous higher-order latent trait $\theta_j$ to the DCM condensation rule at the statement level.

Usage

fit.FCDCM(
  data,
  Q.matrix,
  block.items = NULL,
  dcm.type = "DINA",
  method = c("iStEM", "stan"),
  control.model = NULL,
  control.method = NULL
)

Arguments

data	An $N \times B$ forced-choice data matrix. Each cell must be a two-statement ranking string such as "1>3" (first is selected over second), or 0/1 coded responses where 1 means selecting the first statement in the block.
Q.matrix	An $I \times D$ binary Q-matrix for the statement- level DCM. $q_{id} = 1$ if attribute $d$ is required by statement $i$. Must have exactly $2B$ rows.
block.items	A list of length $B$; each element is an integer vector of the two global statement indices in that block. If NULL, parsed from data assuming "A>B" format where A and B are statement indices.
dcm.type	Character vector (length 1 or $I$): "DINA" (default, conjunctive) or "DINO" (disjunctive). Recycled if scalar.
method	Estimation method: "iStEM" (default) or "stan".
control.model	A named list of model-level hyperparameters for the higher-order DCM. Supported entries: delta1.mu, delta1.sigma Prior mean and SD for $\delta_{1d}$ (higher-order discrimination parameters; log-normal). Defaults: 0, 0.5. The log-normal prior constrains $\delta_{1d} > 0$, reflecting the assumption that higher $\theta_j$ increases the probability of mastering each attribute. Larger values of $\delta_{1d}$ indicate stronger relationships between the continuous trait $\theta_j$ and attribute mastery. delta0.mu, delta0.sigma Prior mean and SD for $\delta_{0d}$ (higher-order difficulty/threshold parameters; normal). Defaults: 0, 1. The probability of mastering attribute $d$ at $\theta_j = 0$ is $\{1+\exp(\delta_{1d}\delta_{0d})\}^{-1}$. Higher values of $\delta_{0d}$ indicate more difficult attributes. eta0.mu, eta0.sigma Prior mean and SD for $\eta_{0b}$ (truncated normal on (0, 0.5)). Defaults: 0.10, 0.10. The prior encodes the expectation that $\eta_{0b}$ is small (rarely choose the lower-mastery statement). etaAB.mu, etaAB.sigma Prior mean and SD for $\eta_{ABb}$ (truncated normal on (0, 0.5)). Defaults: 0.30, 0.10. The prior encodes the expectation that $\eta_{ABb}$ is larger (usually choose the higher-mastery statement). use.prior Logical; if TRUE (default), the normal priors on the block response parameters are applied during iStEM estimation. Set to FALSE to use flat likelihood updates. par Optional initial $B \times 2$ matrix for iStEM with columns eta0 and etaAB. L Theta grid size for numerical integration and the unidimensional iStEM theta Gibbs sampler. Default: 61. theta.lower, theta.upper Bounds for the unidimensional 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). The discrete attribute-profile marginalization can create a challenging posterior; consider increasing to 0.9–0.95 if divergences occur. 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 for burn-in (default: 1.5). eps2 Monte Carlo error tolerance for the final chain (default: 0.4). frac1, frac2 Fractions for the Geweke diagnostic (defaults: 0.1, 0.5). optim.maxit Maximum iterations for parameter optimization steps (default: 50). In FCDCM iStEM, the higher-order $\delta$ parameters and block $\eta$ parameters are updated via closed-form or numerical optimization. fix.corr Logical; in the FCDCM, the latent trait is unidimensional, so this parameter has no effect (default: FALSE). estimate.se Logical; if TRUE (default), standard errors are computed from the final Monte Carlo chain.

Value

An object of class "FCDCM" containing:

npar

Number of free parameters ($2D + 2B$).

method

"stan" or "iStEM".

theta

List with est, se, Rhat ($N \times 1$); higher-order trait estimates.

delta

List with est, se, Rhat ($D \times 2$); higher-order $\delta_1$, $\delta_0$ parameters.

par

List with est, se, Rhat ($B \times 2$); block response parameters eta0 and etaAB.

alpha

List with est, se, Rhat, and prob; est is the binary classified attribute profile ($N \times D$) and prob contains posterior mastery probabilities.

class.post

Posterior probabilities for all $2^D$ attribute profiles.

alpha.patterns, zeta.patterns

Full enumeration of attribute mastery patterns and their condensation outputs.

logLik

Marginal log-likelihood (class "logLik").

Model Specification

Higher-order latent structure. A unidimensional continuous trait $\theta_j \sim N(0, 1)$ governs the probability of mastering each attribute $d = 1, \dots, D$:

$P(\alpha_{jd} = 1 \mid \theta_j) = \frac{1}{1 + \exp\bigl[-\bigl(\delta_{1d}\,\theta_j - \delta_{1d}\,\delta_{0d}\bigr)\bigr]},$

where $\delta_{1d} > 0$ is the discrimination and $\delta_{0d}$ is the difficulty (threshold) for attribute $d$. Attributes are conditionally independent given $\theta_j$.

Statement-level condensation. For each statement $i = 1, \dots, I$ with Q-vector $\mathbf{q}_i = (q_{i1}, \dots, q_{iD})'$, the latent mastery status $\zeta_{ij} \in \{0, 1\}$ is determined by one of two condensation rules:

DINA (conjunctive):

$\zeta_{ij} = \prod_{d: q_{id} = 1} \alpha_{jd}.$

All required attributes must be mastered.

DINO (disjunctive):

$\zeta_{ij} = 1 - \prod_{d: q_{id} = 1} (1 - \alpha_{jd}).$

At least one required attribute must be mastered.

Forced-choice block response. Block $b$ consists of two statements $(A_b, B_b)$. The probability that the respondent selects statement $A_b$ over $B_b$ follows the original FC-DCM specification and fixes equal-condensation probabilities to chance:

$P(\text{choose } A_b \mid \zeta_{A_b}, \zeta_{B_b}) = \begin{cases} \eta_{0b}, & \text{if } \zeta_{A_b} < \zeta_{B_b}, \\ 0.5, & \text{if } \zeta_{A_b} = \zeta_{B_b}, \\ 0.5+\eta_{ABb}, & \text{if } \zeta_{A_b} > \zeta_{B_b}. \end{cases}$

with $0 < \eta_{0b} < 0.5$ and $0 < \eta_{ABb} < 0.5$. Thus equal statement-level mastery produces chance choice, choosing the lower mastery statement has probability $\eta_{0b}$, and choosing the higher mastery statement has probability $0.5+\eta_{ABb}$. The probability of choosing $B_b$ is the complement.

The marginal likelihood integrates over all $2^D$ attribute patterns:

$P(\mathbf{Y}_j \mid \theta_j) = \sum_{\boldsymbol{\alpha} \in \{0,1\}^D} P(\boldsymbol{\alpha} \mid \theta_j) \times \prod_{b=1}^{B} P(Y_{jb} \mid \zeta_{A_b}, \zeta_{B_b}).$

Design and Identifiability

A block with $\zeta_{A_b}=\zeta_{B_b}$ contributes no direct information for distinguishing the two equal condensation states (both are fixed at 0.5). For single-attribute statements and $D=2$, if all blocks compare the same ordered attribute pair, the response distribution distinguishes 10 from 01 but leaves 00 and 11 to be separated only by the higher-order structural model. Good global fit can therefore coexist with weak profile classification.

A forced-choice Q/block design should provide repeated cross-attribute comparisons, avoid pairing statements with identical Q-vectors, and keep each attribute represented in both statement positions. A practical numerical target, consistent with the original FC-DCM simulation logic, is that each attribute appears multiple times as the first/predominant statement and multiple times as the second/inferior statement, with comparisons distributed across different attribute pairs. The automatic simulator uses this balanced assembly rule when block.items is not supplied.

Estimation Methods

Stan (method = "stan"):

Full Bayesian inference via HMC. The discrete attribute profile is exactly marginalized in the Stan likelihood.

iStEM (method = "iStEM"):

Improved Stochastic EM. Person sampling considers the joint space of $\theta_j$ and $\boldsymbol{\alpha}_j$. Parameter updates use closed-form or numerical optimization steps.

References

de la Torre, J., & Douglas, J. A. (2004). Higher-order latent trait models for cognitive diagnosis. Psychometrika, 69(3), 333–353. doi:10.1007/BF02295640

Junker, B. W., & Sijtsma, K. (2001). Cognitive assessment models with few assumptions, and connections with nonparametric item response theory. Applied Psychological Measurement, 25(3), 258–272.

See Also

sim.data.FCDCM, get.fit.index.FCDCM, logLik.FCDCM

Examples

sim <- sim.data.FCDCM(N.person = 20, N.block = 3, D = 2,
                      dcm.type = "DINA")
fit <- fit.FCDCM(sim$data, Q.matrix = sim$Q.matrix,
                 block.items = sim$block.items,
                 dcm.type = "DINA", method = "iStEM",
                 control.method = list(
                   seed = 123, vis = FALSE,
                   M = 2, B = 2, burnin.maxitr = 2,
                   maxitr = 3, eps1 = 10, eps2 = 10,
                   estimate.se = FALSE))
head(fit$alpha$est)                # binary attribute profiles (0/1)
head(fit$alpha$prob)               # marginal mastery probabilities
gof <- get.fit.index(fit)
summary(gof)

---

Fit the Forced-Choice GDINA Model

Description

Fits a forced-choice cognitive diagnostic model (FCGDINA) to comparative response data. The statement-level cognitive diagnosis component can be DINA, DINO, ACDM, or GDINA, and the block-level response model transforms statement endorsement probabilities into forced-choice probabilities for full rankings ("RANK"), most-least choices ("MOLE"), or best-only choices ("PICK").

Usage

fit.FCGDINA(
  data,
  Q.matrix,
  model = c("GDINA", "DINA", "DINO", "ACDM"),
  block.items = NULL,
  fc.type = NULL,
  method = c("EM", "stan", "iStEM"),
  control.model = NULL,
  control.method = NULL
)

Arguments

data	An $N \times B$ forced-choice data matrix.
Q.matrix	An $I \times D$ binary Q-matrix. Rows correspond to statements, columns correspond to attributes, and every row must contain at least one 1.
model	CDM item model: "GDINA", "DINA", "DINO", or "ACDM".
block.items	A list of length $B$. Each element gives the global statement indices in the corresponding forced-choice block. If NULL, the block structure is parsed from full-ranking character data when possible.
fc.type	Forced-choice response type: "RANK", "MOLE", or "PICK". A scalar is recycled to all blocks. If NULL, the type is inferred from character data when possible and otherwise defaults to "RANK".
method	Estimation method. "EM" (default) runs a deterministic posterior-weight EM algorithm, "stan" runs full Bayesian MCMC, and "iStEM" runs improved stochastic EM.
control.model	Optional list of model-level controls. For Stan, delta.mu and delta.sigma set the normal prior $\delta \sim N(\mu, \sigma)$ for the free Stan FCGDINA logit-scale delta effects (defaults: 0 and 1). pi.prior.alpha sets the Dirichlet concentration for the structural parameters $\pi$ (default 1, uniform prior over the simplex). For EM and iStEM, Stan priors are ignored; use delta.lower and delta.upper in control.method instead.
control.method	Optional list of method-specific controls.

Value

An object of class "FCGDINA", a list containing:

npar

Number of free CDM delta parameters.

method

Estimation method actually used: "stan", "iStEM", or "EM".

alpha

List with est, se, Rhat, and prob. est is an $N \times D$ matrix of posterior attribute mastery probabilities; prob contains posterior probabilities for all $2^D$ attribute profiles.

delta

List with est, se, and Rhat. For Stan fits, each element is the probability-scale transform of the identified minimum-norm logit coefficient vector. Absolute statement endorsement probabilities are not identified by forced-choice data alone.

delta.identification

Stan-only rank, nullity, blockwise rank diagnostics, and the identifying convention.

delta.store

Stored posterior, iStEM, or bootstrap delta draws when available.

Q.matrix, block.items, patterns, patterns.total, design.matrix.list, alpha.patterns, fc.type, response

Processed model and data structures used by prediction, fit indices, and S3 methods.

stan.obj, MCMC.obj

Stan fit and extracted posterior draws for method = "stan"; otherwise NULL.

iStEM

iStEM convergence diagnostics and chains for method = "iStEM"; otherwise NULL.

EM

EM convergence diagnostics, log-likelihood trace, and bootstrap summary for method = "EM"; otherwise NULL.

logLik

Marginal log-likelihood with class "logLik".

call, arguments

Matched call and effective fitting arguments.

Model Specification

Let $\boldsymbol{\alpha}_n \in \{0, 1\}^D$ denote the latent attribute profile for person $n$. For statement $i$, the binary Q-vector $\mathbf{q}_i$ selects the required attributes and the reduced attribute pattern is mapped to a CDM design row $\mathbf{x}_{ic}$. For Stan estimation, $\boldsymbol{\delta}_i$ is represented directly as the item-specific CDM effect vector on the logit scale: intercept, attribute main effects, and interaction effects. The item-level endorsement probability for attribute class $c$ is

$p_{ic} = \mathrm{logit}^{-1}(\mathbf{x}_{ic}' \boldsymbol{\delta}_i),$

so main effects and interactions may be negative while probabilities remain in $(0, 1)$. Forced-choice observations identify only within-block utility contrasts, not the absolute endorsement logits of individual statements. The Stan implementation therefore removes the exact null space of those contrasts and reports the unique minimum-norm coefficient representative. The design matrix is generated according to model: "DINA", "DINO", "ACDM", or "GDINA". EM and iStEM retain the existing probability-scale CDM parameterisation.

For each forced-choice block, the statement endorsement probabilities are converted to ranking-pattern probabilities by the same sequential Luce–Plackett transformation used by fit.FCMIRT and fit.FCGGUM. The likelihood marginalizes over all $2^D$ attribute profiles using a uniform latent-class prior:

$P(\mathbf{Y}_n) = \sum_{c=1}^{2^D} \pi_c \prod_{b=1}^{B} P(Y_{nb} \mid \boldsymbol{\alpha}_c, \boldsymbol{\delta}, \mathrm{block}_b).$

Data Coding

data must be an $N \times B$ matrix, with one column per forced-choice block. Entries may be character rankings such as "2>1>3" or 1-based integer pattern indices. If integer pattern indices are supplied, block.items must be provided because item identities cannot be recovered from indices alone. For "MOLE" and "PICK" character data, supplying block.items is recommended and required whenever partial rankings do not identify all items in a block.

Estimation Controls

Common control.method entries are:

seed

Integer random seed.

cores

Number of cores used by Stan chains; ignored by the deterministic EM updates.

vis

Logical progress flag.

delta.lower, delta.upper

Lower and upper bounds for delta optimization in EM and iStEM. Defaults are -4 and 4.

Stan-specific entries include chains, iter, warmup, thin, init, algorithm, and the usual Stan control entries such as adapt_delta, max_treedepth, stepsize, and metric. Defaults are inherited from the common ForceChoice Stan helper: 4 chains, 2000 iterations, and half the iterations used as warmup. The default init = 0 is valid for Stan FCGDINA because it implies $p_{ic} = 0.5$ under the logit link.

iStEM-specific entries include M, B, burnin.maxitr, maxitr, eps1, eps2, frac1, frac2, optim.maxit, and estimate.se. Each iStEM update samples discrete attribute classes from the current posterior and optimises $\delta$ parameters given the sampled class memberships.

EM-specific entries are:

maxitr

Maximum EM iterations; default 500.

minitr

Minimum EM iterations before convergence checks; default 2.

tol

Absolute log-likelihood change tolerance; default 1e-6.

par.tol

Maximum absolute delta-parameter change tolerance; default 1e-4.

optim.maxit

Maximum L-BFGS-B iterations for blockwise delta updates; default 200.

estimate.se

Logical; whether to estimate standard errors by nonparametric bootstrap. Default FALSE.

bootstrap

Number of bootstrap samples when estimate.se = TRUE. Default 100; must be at least 2.

References

de la Torre, J. (2011). The generalized DINA model framework. Psychometrika, 76(2), 179–199. doi:10.1007/s11336-011-9207-7

Templin, J., & Henson, R. A. (2006). Measurement of psychological disorders using cognitive diagnosis models. Psychological Methods, 11(3), 287–305. doi:10.1037/1082-989X.11.3.287

See Also

sim.data.FCGDINA, model.FCGDINA, get.fit.index.FCGDINA, fit.FCDCM, fit.FCMIRT, fit.FCGGUM

Examples

sim <- sim.data.FCGDINA(N.person = 20, N.block = 2, I.block = 2,
                        D = 2, model = "DINA", fc.type = "RANK")

fit <- fit.FCGDINA(
  sim$data,
  Q.matrix = sim$Q.matrix,
  block.items = sim$block.items,
  model = "DINA",
  fc.type = sim$fc.type,
  method = "EM",
  control.method = list(seed = 123, vis = FALSE,
                        maxitr = 2, estimate.se = FALSE)
)

coef(fit, type = "par")
head(fit$alpha$est)
logLik(fit)


# stan code, long time
fit.stan <- fit.FCGDINA(
  sim$data,
  Q.matrix = sim$Q.matrix,
  block.items = sim$block.items,
  model = "GDINA",
  method = "stan",
  control.method = list(seed = 123, cores = 2,
                        chains = 1, iter = 200, warmup = 100)
)

---

Fit the Forced-Choice Generalized Graded Unfolding Model (FCGGUM)

Description

Fits a forced-choice variant of the multidimensional generalized graded unfolding model (GGUM) to comparative ranking data. Unlike dominance-based FC models, the FCGGUM uses an ideal-point (unfolding) item response process: endorsement probability peaks when the person and item locations coincide, and the block-level ranking follows the Luce–Plackett model.

Usage

fit.FCGGUM(
  data,
  Q.matrix = NULL,
  block.items = NULL,
  D = NULL,
  fc.type = NULL,
  method = c("iStEM", "stan"),
  control.model = NULL,
  control.method = NULL
)

Arguments

data	An $N \times B$ forced-choice data matrix with ranking strings or 1-based pattern indices.
Q.matrix	An optional $I \times D$ matrix with entries $-1, 0, 1$. See fit.MGGUM for the sign convention.
block.items	A list of length $B$ giving global item indices in each block. Parsed from data if NULL and data contains ranking strings.
D	Integer; number of latent dimensions. If NULL, inferred from Q.matrix when available, otherwise from the block structure.
fc.type	Forced-choice response type: "RANK" (default), "MOLE", or "PICK".
method	Estimation method: "iStEM" (default) or "stan".
control.model	A named list of model-level hyperparameters for the GGUM unfolding model in a forced-choice framework. Supported entries: a.mu, a.sigma Prior location and scale for $\log a_{id}$ (log-normal). Defaults: 0, 0.5. The log-normal prior ensures positivity of discrimination parameters in the unfolding model. The zero mean-log encourages slopes near 1 unless the data strongly indicate otherwise. delta.pos.mu, delta.pos.sigma Prior mean and SD for $\delta_{id}$ when $q_{id} = 1$ (positive-side item locations; normal). Defaults: 1, 0.5. In the FC context, these priors regularize the location of statements on the positive side of each latent dimension. delta.neg.mu, delta.neg.sigma Prior mean and SD for $\delta_{id}$ when $q_{id} = -1$ (negative-side item locations; normal). Defaults: -1, 0.5. These priors regularize statement locations on the negative side of each dimension. tau.mu, tau.sigma Location and scale for the log-normal prior on the positive threshold magnitude $-\tau_{i1}$. Defaults: 0, 0.5. In the binary FCGGUM, each statement has a single free threshold $\tau_{i1}<0$. theta.mu Prior mean vector for $\boldsymbol{\theta}_j$. Default: rep(0, D). L Theta grid size per dimension for marginal log-likelihood computation and iStEM block Gibbs sampling. Default adapts to $D$. 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). The GGUM likelihood involves summation terms that can create challenging posterior geometries; consider increasing to 0.9–0.95 if divergences occur. 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 for burn-in (default: 1.5). eps2 Monte Carlo error tolerance for the final chain (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). GGUM item parameters ($a$, $\delta$, $\tau$) are optimized jointly under monotonicity and sign constraints. Increase if optimization warnings appear. fix.corr Logical; if TRUE, the inter-trait correlation matrix is fixed to the identity (default: FALSE). estimate.se Logical; if TRUE (default), standard errors are computed from the final Monte Carlo chain. delta.lower Lower bound on $|\delta_{id}|$ (default: 1e-4). Ensures item locations are bounded away from zero for active dimensions. tau.gap.lower, tau.gap.upper Lower and upper bounds on the positive gap used to construct $\tau_{i1}<0$. Defaults are 1e-4 and 6.

Value

An object of class "FCGGUM" containing:

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 (2D + 2)$).

Corr

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

block.items, patterns, patterns.total, response

Block structure.

logLik

Marginal log-likelihood (class "logLik").

Model Specification

Item-level unfolding process. For statement $i$ with Q-vector $\mathbf{q}_i \in \{-1, 0, 1\}^D$, the category response function is the GGUM (see fit.MGGUM for the full polytomous specification). For the forced-choice context, each item is binary ($K_i=2$). Define

$r_{ij} = \left[ \sum_{d=1}^{D} a_{id}^2 \left(\theta_{jd} - \delta_{id}\right)^2 \right]^{1/2}, \qquad S_i = \sum_{d=1}^{D} a_{id}, \qquad \psi_{i1} = \tau_{i1}S_i .$

The binary GGUM endorsement probability used by the code is:

$P_i(\boldsymbol{\theta}_j) = \frac{ \exp\{r_{ij}-\psi_{i1}\} + \exp\{2r_{ij}-\psi_{i1}\} }{ 1 + \exp\{3r_{ij}\} + \exp\{r_{ij}-\psi_{i1}\} + \exp\{2r_{ij}-\psi_{i1}\} },$

where $\tau_{i1} < 0$ is the first (and only free) threshold.

Block-level forced-choice ranking. Within each block $b$, the probability of the observed ranking is computed by applying the sequential Luce/Plackett rule to $u_i(\boldsymbol{\theta}_j)=\mathrm{logit} \{P_i(\boldsymbol{\theta}_j)\}$:

$P(\text{ranking} \mid \boldsymbol{\theta}_j) = \prod_{m=1}^{K_b-1} \frac{\exp\{u_{i_{(m)}}(\boldsymbol{\theta}_j)\}} {\sum_{r=m}^{K_b}\exp\{u_{i_{(r)}}(\boldsymbol{\theta}_j)\}}.$

MOLE and PICK probabilities are sums over compatible full rankings, normalized over the observed block patterns by the C++ probability helper.

Q-matrix for unfolding. The Q-matrix serves a dual role: entries of $\pm 1$ activate both the slope $a_{id}$ and the location $\delta_{id}$; the sign of $q_{id}$ determines the sign of $\delta_{id}$ (positive vs. negative side of dimension $d$).

Estimation Methods

Stan (method = "stan"):

Full Bayesian inference via HMC. The GGUM likelihood involves summation terms that are handled in the Stan model block.

iStEM (method = "iStEM"):

Improved Stochastic EM. Person sampling uses finite-grid block Gibbs; GGUM item parameters updated via constrained L-BFGS-B with monotonicity constraints on $\tau$ and sign constraints on $\delta$.

References

Roberts, J. S., Donoghue, J. R., & Laughlin, J. E. (2000). A general item response theory model for unfolding unidimensional polytomous responses. Applied Psychological Measurement, 24(1), 3–32.

See Also

sim.data.FCGGUM, get.fit.index.FCGGUM, logLik.FCGGUM, fit.MGGUM

Examples

sim <- sim.data.FCGGUM(N.person = 20, N.block = 3, I.block = 2,
                       D = 2, fc.type = "RANK")
fit <- fit.FCGGUM(sim$data, Q.matrix = sim$Q.matrix,
                  block.items = sim$block.items,
                  D = 2, fc.type = "RANK", 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)
gof <- get.fit.index(fit)
summary(gof)

---

Fit the Forced-Choice Multidimensional IRT (FCMIRT) Model

Description

Fits a forced-choice variant of the multidimensional item response theory model to comparative ranking data. In forced-choice (FC) formats, respondents compare items within blocks rather than rating each item in isolation. The FCMIRT model couples an item-level MIRT endorsement model with a block-level ranking mechanism to recover latent trait estimates from comparative responses.

Usage

fit.FCMIRT(
  data,
  model = "2PL",
  Q.matrix = NULL,
  block.items = NULL,
  D = NULL,
  fc.type = NULL,
  method = c("iStEM", "stan"),
  control.model = NULL,
  control.method = NULL
)

Arguments

data	An $N \times B$ forced-choice data object. Each entry is either a ranking string (e.g., "2>1>3") or a 1-based pattern index into the block's ranking patterns. If pattern indices are used, block.items must be supplied.
model	MIRT model type: "m1pl", "m2pl" (default), "m3pl", or "m4pl".
Q.matrix	An optional $I \times D$ binary (0/1) matrix indicating active item loadings.
block.items	A list of length $B$ where each element is an integer vector of global item indices belonging to that block. If NULL and data contains ranking strings, automatically parsed from the data.
D	Integer; number of latent dimensions. If NULL, inferred from Q.matrix; if both are NULL, an error is raised.
fc.type	Character vector (length 1 or $B$): "RANK" (full ranking, default), "MOLE" (most-least), or "PICK" (best only). Recycled to all blocks if scalar.
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, 1.0. Important: The Stan default for a.sigma is deliberately set to 1.0 (wider than the MIRT default of 0.25) because forced-choice blocks identify relative utilities; tight common slope priors can spuriously shrink within-block slopes toward equality, degrading trait recovery. For iStEM, the prior defaults match the standard MIRT values (a.mu = 0.25, a.sigma = 0.25) unless overridden. b.mu, b.sigma Prior mean and SD for $b_i$ (normal). Defaults: 0, 1. Under the block zero-sum constraint $\sum_{i \in \text{block } b} b_i = 0$, these priors regularize the free difficulty parameters. c.mu, c.sigma Uniform support bounds $[c_{\min}, c_{\max}]$ for lower-asymptote parameters $c_i$. Defaults: 0, 0.35. Only used for 3PL/4PL models. d.mu, d.sigma Uniform support bounds $[d_{\min}, d_{\max}]$ for upper-asymptote parameters $d_i$. Defaults: 0.65, 1. Only used for 4PL models. theta.mu Prior mean vector for $\boldsymbol{\theta}_j$. Default: rep(0, D). 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). theta.lower, theta.upper Bounds for the theta grid used by marginal log-likelihood computation and iStEM block Gibbs sampling. Defaults: -6, 6. use.prior Logical (iStEM only). If FALSE, the item-prior penalty term is omitted from the item update, recovering an approximate maximum-likelihood item step. Default is TRUE.
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). Forced-choice models often have challenging posterior geometries; consider increasing to 0.9 or higher if divergent transitions appear. 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 for burn-in (default: 1.5). eps2 Monte Carlo error tolerance for the final chain (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}$ during optimization (defaults: 1e-4, 6). b.lower, b.upper Bounds on $b_i$ during optimization (defaults: -6, 6). c.lower, c.upper Bounds on $c_i$ for 3PL/4PL models. Defaults from prior support. d.lower, d.upper Bounds on $d_i$ for 4PL models. Defaults from prior support.

Value

An object of class "FCMIRT" containing:

npar

Number of free parameters.

method

"stan" or "iStEM".

theta

List with est, se, Rhat ($N \times D$); person trait estimates.

par

List with est, se, Rhat, free ($I \times (D+3)$); item parameters a1..aD, b, c, d.

Corr

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

block.items, patterns, patterns.total, response

Block structure and response data.

logLik

Marginal log-likelihood (class "logLik").

Model Specification

Item-level endorsement. Each statement $i = 1, \dots, I$ is governed by a MIRT model (1PL through 4PL). For a person with trait $\boldsymbol{\theta}_j$, the probability of endorsing statement $i$ in isolation is:

$P_i(\boldsymbol{\theta}_j) = c_i + (d_i - c_i) \times \frac{1}{1 + \exp\bigl[-\bigl(\sum_{d=1}^{D} a_{id}\,\theta_{jd} - b_i\bigr)\bigr]}.$

Block-level ranking. Let block $b = 1, \dots, B$ consist of items $\{i_1, i_2, \dots, i_{K_b}\}$ and define the implemented statement utility as the logit of the item-level endorsement probability,

$u_i(\boldsymbol{\theta}_j) = \log\left[ \frac{P_i(\boldsymbol{\theta}_j)} {1-P_i(\boldsymbol{\theta}_j)} \right].$

The probability that an examinee produces the ranking $i_{(1)} \succ i_{(2)} \succ \dots \succ i_{(K_b)}$ (read as "$i_{(1)}$ is preferred over $i_{(2)}$ over ...") is:

$P\bigl(i_{(1)} \succ \dots \succ i_{(K_b)} \mid \boldsymbol{\theta}_j\bigr) = \prod_{m=1}^{K_b-1} \frac{ \exp\{u_{i_{(m)}}(\boldsymbol{\theta}_j)\} }{ \sum_{r=m}^{K_b} \exp\{u_{i_{(r)}}(\boldsymbol{\theta}_j)\} }.$

The Stan code evaluates the same choice kernel as $\log P_i + \sum_{h \ne i}\log(1-P_h)$ within each remaining set, which is algebraically equivalent to a softmax over $\mathrm{logit}(P_i)$.

Partial rankings (MOLE / PICK). For MOLE (most-least) data, only the best and worst items in each block are identified; for PICK data, only the best item. The probability of a partial ranking is obtained by marginalizing (summing) the full-ranking probabilities over all completions consistent with the partial constraint; the C++ probability helper returns probabilities normalized over the observed RANK/MOLE/PICK patterns for each block.

Identification constraint (2PLM-RANK). Block intercepts $b_i$ are constrained to sum to zero within each block:

$\sum_{i \in \text{block } b} b_i = 0, \qquad b = 1, \dots, B,$

so that only $K_b - 1$ difficulties are freely estimated per block.

Estimation Methods

Stan (method = "stan"):

Full Bayesian inference via HMC. The forced-choice likelihood is evaluated over the ranking pattern probabilities. By default, a.sigma is set to 1.0 (wider prior) to avoid spurious shrinkage of within-block slopes toward equality.

iStEM (method = "iStEM"):

Improved Stochastic EM with finite-grid block Gibbs person-sampling and item optimization (L-BFGS-B). Block-level identification constraints are enforced during each item-parameter update.

References

Brown, A., & Maydeu-Olivares, A. (2011). Item response modeling of forced-choice questionnaires. Educational and Psychological Measurement, 71(3), 460–502. doi:10.1177/0013164410375112

Luce, R. D. (1959). Individual choice behavior: A theoretical analysis. Wiley.

Plackett, R. L. (1975). The analysis of permutations. Journal of the Royal Statistical Society: Series C, 24(2), 193–202.

See Also

sim.data.FCMIRT, get.fit.index.FCMIRT, logLik.FCMIRT, fit.MIRT

Examples

# Simulate forced-choice data
sim <- sim.data.FCMIRT(N.person = 20, N.block = 3, I.block = 2,
                       D = 2, model = "m2pl", fc.type = "RANK")

# Fit via iStEM
fit <- fit.FCMIRT(sim$data, model = "m2pl",
                  Q.matrix = sim$Q.matrix,
                  block.items = sim$block.items,
                  D = 2, fc.type = "RANK",
                  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))

# Examine trait recovery
cor(fit$theta$est, sim$theta)

# Fit indices (nominal binary expansion)
gof <- get.fit.index(fit)
summary(gof)

---

Fit the Multidimensional Generalized Graded Unfolding Model (MGGUM)

Description

Fits a multidimensional generalization of the generalized graded unfolding model (GGUM) to polytomous response data. Unlike dominance (cumulative) IRT models, the MGGUM belongs to the ideal-point (unfolding) model family where the probability of endorsement is highest when the person and item locations match, decreasing as the person moves away in either direction. Two estimation backends are provided: Bayesian MCMC (Stan) and iStEM.

Usage

fit.MGGUM(
  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$.
D	Integer; number of latent dimensions ($D \ge 1$). Default is 2.
Q.matrix	An optional $I \times D$ matrix with entries $-1, 0, 1$. See section Q-matrix sign convention. Default generates a random signed matrix.
length.poly	Integer vector (length $I$) or scalar giving the number of categories per item. If NULL, inferred from the observed response maxima.
method	Estimation method: "iStEM" (default) or "stan".
control.model	A named list of model-level hyperparameters for the GGUM unfolding model. Supported entries: a.mu, a.sigma Prior location and scale for $\log a_{id}$ (log-normal). Defaults: 0, 0.5. The log-normal prior ensures positivity of discrimination parameters. The zero mean-log encourages slopes near 1 unless the data strongly indicate otherwise. delta.pos.mu, delta.pos.sigma Prior mean and SD for $\delta_{id}$ when $q_{id} = 1$ (positive-side item locations; normal). Defaults: 1, 0.5. These priors regularize item locations on the positive side of each dimension. delta.neg.mu, delta.neg.sigma Prior mean and SD for $\delta_{id}$ when $q_{id} = -1$ (negative-side item locations; normal). Defaults: -1, 0.5. These priors regularize item locations on the negative side of each dimension. tau.mu, tau.sigma Location and scale for the log-normal prior on positive threshold magnitudes $-\tau_{ik}$ for $k \ge 1$. Defaults: 0, 0.5. The thresholds satisfy $\tau_{i1} < \tau_{i2} < \dots < \tau_{i,K_i-1} < 0$ with $\tau_{i0} = 0$ fixed for identification. theta.mu Prior mean vector for $\boldsymbol{\theta}_j$. Default: rep(0, D). L Theta grid size per dimension for marginal log-likelihood computation and iStEM block Gibbs sampling. Default adapts to $D$. 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). The burn-in phase ends when $\sum z^2 / K < \epsilon_1$ or the MC error criterion is satisfied. eps2 Monte Carlo error tolerance (default: 0.4). The algorithm continues until $max(d_k \cdot N) < \epsilon_2$. 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 during the item-parameter update (default: 50). Increase if item optimization does not converge within this limit. fix.corr Logical; if TRUE, the inter-trait correlation matrix is fixed to the identity (default: FALSE). estimate.se Logical; if TRUE (default), standard errors are computed from the final Monte Carlo chain. delta.lower Lower bound on $|\delta_{id}|$ (default: 1e-4). Ensures item locations are bounded away from zero for active dimensions. tau.gap.lower, tau.gap.upper Lower and upper bounds on the positive gaps used to construct ordered thresholds. Defaults are 1e-4 and 6.

Value

An object of class "MGGUM" containing:

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 (2D + K_{max})$). Columns: a1..aD, delta1..deltaD, tau0..tau_{K_{max}-1}.

Corr

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

length.poly

Per-item category counts.

logLik

Marginal log-likelihood (class "logLik").

Model Specification

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

The category response probability implemented by model.MGGUM, iStEM, and Stan is the distance-based MGGUM form. Define

$r_{ij} = \left[ \sum_{d=1}^{D} a_{id}^2 \left(\theta_{jd} - \delta_{id}\right)^2 \right]^{1/2}, \qquad S_i = \sum_{d=1}^{D} a_{id}.$

Let $\tau_{i0}=0$ and $\psi_{ik}=S_i\sum_{v=0}^{k}\tau_{iv}$. With $M_i=2K_i-1$, the probability of category $k = 0,\dots,K_i-1$ is

$P(Y_{ij} = k \mid \boldsymbol{\theta}_j) = \frac{ \exp\{k r_{ij}-\psi_{ik}\} + \exp\{(M_i-k)r_{ij}-\psi_{ik}\} }{ \sum_{r=0}^{K_i-1} \bigl[ \exp\{r r_{ij}-\psi_{ir}\} + \exp\{(M_i-r)r_{ij}-\psi_{ir}\} \bigr] }.$

Here $a_{id} \ge 0$ are discrimination parameters, $\delta_{id}$ are item-location parameters (signed positive when $q_{id} = 1$, negative when $q_{id} = -1$), and threshold parameters $\tau_{ik}$ satisfying $\tau_{i0} = 0$ and $\tau_{i1} < \tau_{i2} < \dots < \tau_{i,K_i-1} < 0$.

Q-matrix sign convention:

$q_{id} = 1$

Active dimension; $\delta_{id} > 0$ (item located on the positive side of dimension $d$).

$q_{id} = -1$

Active dimension; $\delta_{id} < 0$ (item located on the negative side of dimension $d$).

$q_{id} = 0$

Inactive dimension; $a_{id} = \delta_{id} = 0$.

The sign of $q_{id}$ controls the item-location side, not the sign of the slope. Both $a_{id}$ and $|\delta_{id}|$ reflect the item's relevance to dimension $d$.

Estimation Methods

Stan (method = "stan"):

Full Bayesian inference via HMC. All parameters are jointly sampled from the posterior distribution.

iStEM (method = "iStEM"):

Improved Stochastic EM with finite-grid block Gibbs person-sampling and item optimization via L-BFGS-B. The threshold parameters $\tau_{ik}$ are constrained to maintain monotonicity during optimization.

References

Roberts, J. S., Donoghue, J. R., & Laughlin, J. E. (2000). A general item response theory model for unfolding unidimensional polytomous responses. Applied Psychological Measurement, 24(1), 3–32. doi:10.1177/01466216000241001

Usami, S. (2011). Generalized graded unfolding model with structural equation for subject parameters. Japanese Psychological Research, 53(3), 221–232.

See Also

sim.data.MGGUM, get.fit.index.MGGUM, logLik.MGGUM

Examples

sim <- sim.data.MGGUM(N = 20, I = 6, D = 2, length.poly = 4)
fit <- fit.MGGUM(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)
gof <- get.fit.index(fit)
summary(gof)

---

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. 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)

---