Type:	Package
Title:	Robust Categorical Data Analysis
Version:	0.2
Date:	2026-02-10
Maintainer:	Max Welz <max.welz@uzh.ch>
Description:	Robust categorical data analysis based on the theory of C-estimation developed in Welz (2024) <doi:10.48550/arXiv.2403.11954>. For now, the package only implements robust estimation of polychoric correlation as proposed in Welz, Mair and Alfons (2026) <doi:10.1017/psy.2025.10066> and robust estimation of polyserial correlation (Welz, 2026 <doi:10.1017/psy.2026.10091>) with methods for printing and plotting. We will implement further models in future releases. In addition, the package is still experimental, so input arguments and class structure may change in future releases.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
Depends:	ggplot2
Imports:	Rcpp (≥ 1.0.10), stats, mvtnorm, stringr, parallel, Matrix, numDeriv, pracma, utils
Suggests:	testthat (≥ 3.0.0)
LinkingTo:	Rcpp
RoxygenNote:	7.3.2
NeedsCompilation:	yes
Packaged:	2026-02-10 10:45:25 UTC; mwelz
Author:	Max Welz [aut, cre], Andreas Alfons [aut], Patrick Mair [aut]
Repository:	CRAN
Date/Publication:	2026-02-10 11:50:02 UTC

Neutral initialization of starting values for polychoric correlation

Description

Initializes starting values for numerical optimization in a neutral way. The optimization problem itself is convex, so the initialization should not matter much.

Usage

initialize_param(x, y)

Arguments

Value

A vector of initial values for the polychoric correlation coefficient, the X-threshold parameters, and the Y-threshold parameters

Examples

## example data
set.seed(123)
x <- sample(c(1,2,3), size = 100, replace = TRUE)
y <- sample(c(1,2,3), size = 100, replace = TRUE)
initialize_param(x, y)

Plot method for classes "robpolycor" and "polycor".

Description

Plot method for classes "robpolycor" and "polycor".

Usage

## S3 method for class 'robpolycor'
plot(x, cutoff = 3, ...)

Arguments

Value

An object of class "ggplot".

Examples

## example data
set.seed(123)
x <- sample(c(1,2,3), size = 100, replace = TRUE)
y <- sample(c(1,2,3), size = 100, replace = TRUE)

fit <- polycor(x,y) 
plot(fit)

Robust estimation of polychoric correlation

Description

Implements to robust estimator of Welz, Mair and Alfons (2024, doi:10.48550/arXiv.2407.18835) for the polychoric correlation model, based on the general theory of C-estimation proposed by Welz (2024, doi:10.48550/arXiv.2403.11954).

Usage

polycor(
  x,
  y = NULL,
  c = 0.6,
  variance = TRUE,
  constrained = "ifneeded",
  method = NULL,
  maxcor = 0.999,
  tol_thresholds = 0.01,
  init = initialize_param(x, y)
)

Arguments

Value

An object of class "robpolycor", which is a list with the following components.

thetahat

A vector of estimates for the polychoric correlation coefficient (rho) as well as thresholds for x (named a1,a2,...,a_{Kx-1}) and y (named b1,b2,...,b_{Ky-1}).

stderr

A vector of standard errors for each estimate in thetahat.

vcov

Estimated asymptotic covariance matrix of thetahat. The matrix \Sigma in the paper (asymptotic covariance matrix of \sqrt{N} \hat{\theta}) can be obtained via vcov * N, where N is the sample size.

chisq,pval,df

Currently NULL, will in a future release be the test statistic, p-value, and degrees of freedom of a test for bivariate normality.

objective

Value of minimized loss function.

optim

Object of class optim.

Examples

## example data
set.seed(123)
x <- sample(c(1,2,3), size = 100, replace = TRUE)
y <- sample(c(1,2,3), size = 100, replace = TRUE)

polycor(x,y)     # robust
polycor_mle(x,y) # non-robust MLE

Maximum likelihood estimation of polychoric correlation coefficient

Description

Implements the maximum likelihood estimator of Olsson (1979, Psychometrika, doi:10.1007/BF02296207) for the polychoric correlation model.

Usage

polycor_mle(
  x,
  y = NULL,
  variance = TRUE,
  constrained = "ifneeded",
  twostep = FALSE,
  method = NULL,
  maxcor = 0.999,
  tol_thresholds = 0.01,
  init = initialize_param(x, y)
)

Arguments

Value

An object of class "robpolycor". See polycor() for details.

Examples

## example data
set.seed(123)
x <- sample(c(1,2,3), size = 100, replace = TRUE)
y <- sample(c(1,2,3), size = 100, replace = TRUE)

polycor(x,y)     # robust
polycor_mle(x,y) # non-robust MLE

Robust estimation of polychoric correlation matrix

Description

A useful wrapper of polycor to robustly estimate a polychoric correlation matrix by calculating all unique pairwise polychoric correlation coefficients.

Usage

polycormat(
  data,
  c = 0.6,
  parallel = FALSE,
  num_cores = 1L,
  return_polycor = TRUE,
  variance = TRUE,
  constrained = "ifneeded",
  method = NULL,
  maxcor = 0.999,
  tol_thresholds = 0.01
)

Arguments

Value

If return_polycor = TRUE, returns a list with a polychoric correlation matrix and list of "polycor" objects. If return_polycor = FALSE, then only a correlation matrix is returned.

Examples

## example data
set.seed(123)
data <- matrix(sample(c(1,2,3), size = 3*100, replace = TRUE), nrow = 100)
polycormat(data)     # robust 
polycormat_mle(data) # non-robust MLE

Maximum likelihood estimation of polychoric correlation matrix

Description

A useful wrapper of polycor_mle to estimate a polychoric correlation matrix via maximum likelihood by calculating all unique pairwise polychoric correlation coefficients.

Usage

polycormat_mle(
  data,
  parallel = FALSE,
  num_cores = 1L,
  return_polycor = TRUE,
  variance = TRUE,
  constrained = "ifneeded",
  method = NULL,
  maxcor = 0.999,
  tol_thresholds = 0.01
)

Arguments

Value

If return_polycor = TRUE, returns a list with a polychoric correlation matrix and list of "polycor" objects. If return_polycor = FALSE, then only a correlation matrix is returned.

Examples

## example data
set.seed(123)
data <- matrix(sample(c(1,2,3), size = 3*100, replace = TRUE), nrow = 100)
polycormat(data)     # robust 
polycormat_mle(data) # non-robust MLE

Robust estimation of polyserial correlation

Description

Implements the robust estimator of Welz (2025, doi:10.48550/arXiv.2510.15632) for the polyserial correlation model.

Usage

polyserial(
  x,
  y,
  alpha = 0.5,
  num_y = max(y),
  constrained = "ifneeded",
  method = NULL,
  variance = TRUE,
  weights = TRUE,
  init = polyserial_initialize_param(x = x, num_y = num_y, robust = TRUE),
  maxcor = 0.999,
  tol_thresholds = 0.01,
  tol_sigma2 = 0.01
)

Arguments

Value

An object of class "polyserial", which is a list with the following components.

thetahat

A vector of estimates for the polyserial correlation coefficient (rho), population mean of X (mu), population variance of Y (sigma2), as well as thresholds for y (named tau1,tau2,...,tau_{r-1}).

stderr

A vector of standard errors for each estimate in thetahat.

vcov

Estimated asymptotic covariance matrix of thetahat. The matrix \Sigma in the paper (asymptotic covariance matrix of \sqrt{N} \hat{\theta}) can be obtained via vcov * N, where N is the sample size.

pointpolyserial

Estimated polyserial correlation coefficient, calculated with provided scoring of Y

weights

List of rescaled and raw outlyingness weights for each observation as well as maximum possible raw weight that was used for rescaling (sup).

objective

Value of minimized loss function.

optim

Object of class optim.

inputs

List of provided inputs.

Examples

## example data
set.seed(123)
x <- rnorm(n = 100)
y <- sample(c(1,2), size = 100, replace = TRUE)

polyserial(x,y)

Efficiency of minimum DPD estimators of polyserial model

Description

Calculate population asymptotic variance-covariance matrix associated with a parameter vector theta, assuming that the polyserial model is correctly specified and that theta is the true model parameter. May take a few moments to compute because a relatively large number of integrals need to be numerically solved.

Usage

polyserial_efficiency(theta, alpha)

Arguments

Value

A numeric matrix that is the population asymptotic variance-covariance matrix associated with a parameter vector theta and tuning constant alpha.

Examples

theta <- c(rho = 0, mu = 0, sigma2 = 1, tau1 = 0) # true parameter vector
polyserial_efficiency(theta, alpha = 0.5)

Neutral initialization of starting values for polyserial correlation

Description

Initializes starting values for numerical optimization in a neutral way.

Usage

polyserial_initialize_param(x, num_y, robust = FALSE)

Arguments

Value

A vector of initial values for the polyserial correlation coefficient, mu, sigma2, and Y-threshold parameters

Examples

## example data
set.seed(123)
x <- rnorm(100)
polyserial_initialize_param(x = x, num_y = 3)

Maximum likelihood estimation of polyserial correlation

Description

Implements maximum likelihood estimation of the polyserial correlation model.

Usage

polyserial_mle(
  x,
  y,
  num_y = max(y),
  constrained = "ifneeded",
  method = NULL,
  variance = TRUE,
  init = polyserial_initialize_param(x = x, num_y = num_y, robust = FALSE),
  maxcor = 0.999,
  tol_thresholds = 0.01,
  tol_sigma2 = 0.01
)

Arguments

Value

An object of class "polyserial", which is a list with the following components.

thetahat

A vector of estimates for the polyserial correlation coefficient (rho), population mean of X (mu), population variance of Y (sigma2), as well as thresholds for y (named tau1,tau2,...,tau_{r-1}).

stderr

A vector of standard errors for each estimate in thetahat.

vcov

Estimated asymptotic covariance matrix of thetahat. The matrix \Sigma in the paper (asymptotic covariance matrix of \sqrt{N} \hat{\theta}) can be obtained via vcov * N, where N is the sample size.

pointpolyserial

Estimated polyserial correlation coefficient, calculated with provided scoring of Y

objective

Value of minimized loss function.

optim

Object of class optim.

inputs

List of provided inputs.

Examples

## example data
set.seed(123)
x <- rnorm(n = 100)
y <- sample(c(1,2), size = 100, replace = TRUE)

polyserial_mle(x,y)

Print method for classes "robpolycor" and "polycor".

Description

Print method for classes "robpolycor" and "polycor".

Usage

## S3 method for class 'robpolycor'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

A print to the console.

Examples

set.seed(123)
x <- sample(c(1,2,3), size = 100, replace = TRUE)
y <- sample(c(1,2,3), size = 100, replace = TRUE)
fit <- polycor(x,y) 

print(fit)
fit # equivalent

Print method for classes "robpolyserial" and "polyserial".

Description

Print method for classes "robpolyserial" and "polyserial".

Usage

## S3 method for class 'robpolyserial'
print(x, digits = max(3L, getOption("digits") - 3L), ...)

Arguments

Value

A print to the console.

Examples

set.seed(123)
x <- rnorm(100)
y <- sample(c(1,2), size = 100, replace = TRUE)
fit <- polyserial(x,y) 

print(fit)
fit # equivalent

Obtain estimated asymptotic variance-covariance matrix

Description

Estimate asymptotic variance-covariance matrix of polychoric model.

Usage

## S3 method for class 'robpolycor'
vcov(object, ...)

Arguments

Details

Method for classes "robpolycor" and "polycor". Returns the estimated asymptotic variance-covariance matrix of a point estimate thetahat. The matrix \Sigma in the paper (asymptotic variance-covariance matrix of \sqrt{N} \hat{\theta}) can be obtained via multiplying the returned matrix by the sample size.

Value

A numeric matrix, being the estimated asymptotic covariance matrix for the model parameters

Examples

set.seed(123)
x <- sample(c(1,2,3), size = 100, replace = TRUE)
y <- sample(c(1,2,3), size = 100, replace = TRUE)
fit <- polycor(x,y) 

vcov(fit)

Obtain estimated asymptotic variance-covariance matrix

Description

Estimate asymptotic variance-covariance matrix of polyserial model.

Usage

## S3 method for class 'robpolyserial'
vcov(object, ...)

Arguments

Details

Method for classes "robpolyserial" and "polyserial". Returns the estimated asymptotic variance-covariance matrix of a point estimate thetahat. The matrix \Sigma in the paper (asymptotic variance-covariance matrix of \sqrt{N} \hat{\theta}) can be obtained via multiplying the returned matrix by the sample size.

Value

A numeric matrix, being the estimated asymptotic covariance matrix for the model parameters

Examples

## example data
set.seed(123)
x <- rnorm(n = 100)
y <- sample(c(1,2), size = 100, replace = TRUE)

fit <- polyserial(x,y)
vcov(fit)