3.1 The polyserial correlation model

Suppose we observe a continuous real-valued random variable X with unknown population mean $\mu = \mu _X \in \mathbb {R}$ and unknown population variance $\sigma ^2 = \sigma ^2_X> 0$ . In addition, suppose that we also observe a polytomous ordinal random variable Y that takes values in some finite set $\mathcal {Y}$ of known cardinality r. Without loss of generality, we assume throughout this article that $\mathcal {Y} = \{1,2,\dots , r\}$ . Further assume that there exists a latent continuous random variable $\eta $ governing the ordinal variable Y through the unobserved discretization process

(3.1) $$ \begin{align} Y = \begin{cases} 1 & \mathrm{ if }\ \eta < \tau_1, \\2 & \mathrm{ if }\ \tau_1 \leq \eta < \tau_2,\\3 & \mathrm{ if }\ \tau_2 \leq \eta < \tau_3,\\\vdots & \\r & \mathrm{ if }\ \tau_{r-1} \leq \eta, \end{cases} \end{align} $$

where $-\infty < \tau _1 < \tau _2 < \cdots < \tau _{r-1} < +\infty $ are fixed but unknown threshold parameters. In practice, Y often denotes the responses to a Likert-type rating item with r response categories.

The primary object of interest is the fixed but unknown population correlation between the observed X and the latent $\eta $ , denoted by

$$\begin{align*}\rho = \mathbb{C}\mathrm{or} \left[X,\ \eta\right]. \end{align*}$$

To identify the correlation coefficient $\rho $ , one usually assumes that X and $\eta $ are jointly normally distributed according to

(3.2) $$ \begin{align} \begin{pmatrix} X \\ \eta \end{pmatrix} \sim \mathrm{N}_2 \left( \begin{pmatrix} \mu \\ 0 \end{pmatrix} , \begin{pmatrix} \sigma^2 & \rho\sigma \\ \rho\sigma & 1 \end{pmatrix} \right). \end{align} $$

The bivariate normality model (3.2) implies the marginal normality properties $X\sim \mathrm {N}(\mu , \sigma ^2)$ and $\eta \sim \mathrm {N}(0,1)$ . It also identifies the correlation coefficient $\rho \in (-1,1)$ through the familiar identity $\mathbb {C}\mathrm {or} \left [X,\ \eta \right ] = \mathbb {C}\mathrm {ov} \left [X,\ \eta \right ]\big /\sqrt {\mathbb {V}\mathrm {ar} \left [ X \right ]\mathbb {V}\mathrm {ar} \left [ \eta \right ]} = \rho $ . Note that since the variable $\eta $ is unobserved, its population mean and variance are not jointly identifiable, which is why they are fixed to 0 and 1, respectively.

Combining the discretization model (3.1) with the bivariate normality model (3.2) yields the polyserial correlation model (Pearson, Reference Pearson1913), or, in short, polyserial model. In this model, the correlation parameter $\rho = \mathbb {C}\mathrm {or} \left [X,\ \eta \right ]$ is referred to as the polyserial correlation coefficient. If Y is dichotomous, the polyserial model reduces to the biserial correlation model of Pearson (Reference Pearson1909). The theoretical properties of biserial correlation are studied by Tate (Reference Tate1955a, Reference Tate1955b) as well as Jaspen (Reference Jaspen1946), and those of polyserial correlation by Olsson et al. (Reference Olsson, Drasgow and Dorans1982).

The polyserial model is subject to $d = r + 2$ parameters, namely, the mean and variance parameters $(\mu , \sigma ^2)$ of the observed X, the polyserial correlation coefficient $\rho $ from the normality model (3.2), as well as the $r-1$ thresholds from the discretization process (3.1) of the latent $\eta $ . These parameters are jointly collected in a d-dimensional parameter vector

$$\begin{align*}\boldsymbol{\theta} = \left(\rho, \mu, \sigma^2, \boldsymbol{\tau}^\top \right)^\top, \end{align*}$$

where the vector $\boldsymbol {\tau } = (\tau _1,\dots , \tau _{r-1})^\top $ contains the $r-1$ thresholds.

Under the polyserial model evaluated at a parameter vector $\boldsymbol {\theta }\in \boldsymbol {\Theta }$ , the density of the observed-latent pair $(X,\eta )$ of continuous variables is given by

$$\begin{align*}{p}_{X\eta}^{}\left(x,v; \boldsymbol{\theta}\right) = \phi_2 \left( \begin{pmatrix} x \\ v \end{pmatrix}; \begin{pmatrix} \mu \\ 0 \end{pmatrix} , \begin{pmatrix} \sigma^2 & \rho\sigma \\ \rho\sigma & 1 \end{pmatrix} \right), \qquad x,v\in\mathbb{R}, \end{align*}$$

where $\phi _2(\cdot; \boldsymbol {m}, \boldsymbol {S})$ is the density of the bivariate normal distribution with population mean $\boldsymbol {m}$ and covariance matrix $\boldsymbol {S}$ . Further, denote by ${P}_{X\eta }^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )$ the distribution function corresponding to the normal density ${p}_{X\eta }^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )$ .

For the observed variables $(X,Y)$ , the joint density at a realization $x\in \mathbb {R}$ and a response $y\in \mathcal {Y} = \{1,\dots , r\}$ of the ordinal Y under the polyserial model at parameter $\boldsymbol {\theta }$ reads

(3.3) $$ \begin{align} {p}_{XY}^{}\left(x,y; \boldsymbol{\theta}\right) = \int_{\tau_{y-1}}^{\tau_y} {p}_{X\eta}^{}\left(x,v; \boldsymbol{\theta}\right) \mathrm{d} v, \end{align} $$

where we adopt the conventions $\tau _0 = -\infty $ and $\tau _r = +\infty $ .Footnote ² The joint distribution function of the observed $(X,Y)$ under the polyserial model can now be expressed as

(3.4) $$ \begin{align} {P}_{XY}^{}\left(x,y; \boldsymbol{\theta}\right) = \mathbb{P}_{\boldsymbol{\theta}} \left[ X \leq x, Y \leq y \right] = \int_{-\infty}^x \sum_{w\leq y} {p}_{XY}^{}\left(u,w; \boldsymbol{\theta}\right) \mathrm{d} u = \int_{-\infty}^{x}\int_{-\infty}^{\tau_y} {p}_{X\eta}^{}\left(u,v; \boldsymbol{\theta}\right) \mathrm{d} v\mathrm{d} u , \end{align} $$

where the third equality follows from (3.3) in conjunction with the linearity of the integral operator. As such, ${P}_{XY}^{}\left (x,y; \boldsymbol {\theta }\right )$ is the distribution function associated with density ${p}_{XY}^{}\left (x,y; \boldsymbol {\theta }\right )$ , so we refer to it as the polyserial model distribution.

Our expressions for the polyserial model distribution and density are different but equivalent to the more commonly used expressions in Olsson et al. (Reference Olsson, Drasgow and Dorans1982), which are provided in Section A of the Supplementary Material.

3.2 Point polyserial correlation

In addition to the correlation between the observed X and the latent $\eta $ , one might also be interested in the correlation between X and the observed ordinal Y. The correlation between the observed X and Y is known as point polyserial correlation. In order to identify the desired point polyserial correlation $\mathbb {C}\mathrm {or} \left [X,\ Y\right ]$ , one needs to assign a numerical interpretation to the r answer categories of Y, that is, introduce a scoring system. Given a scoring system, the point polyserial correlation coefficient $\widetilde {\rho } = \mathbb {C}\mathrm {or} \left [X,\ Y\right ]$ is identified by the polyserial model and can be estimated by using estimates of the polyserial model parameters $\boldsymbol {\theta }$ . We provide details in Section A of the Supplementary Material.

3.3 Maximum likelihood estimation

Suppose we observe a sample $\{(X_i, Y_i)\}_{i=1}^N$ of N independent copies of $(X,Y)$ generated by the polyserial model at some true parameter vector $\boldsymbol {\theta }_{\ast} = \left (\rho _{\ast}, \mu _{\ast}, \sigma ^2_{\ast}, \boldsymbol {\tau }_{\ast}^\top \right )^\top $ . The statistical problem is to estimate the true $\boldsymbol {\theta }_{\ast}$ from the observed sample, which is traditionally achieved by the ML estimator proposed by Cox (Reference Cox1974) and Olsson et al. (Reference Olsson, Drasgow and Dorans1982).

The ML estimator (MLE) of $\boldsymbol {\theta }_{\ast}$ is defined as the log-likelihood maximizer

(3.5) $$ \begin{align} \widehat{\boldsymbol{\theta}}_N^{\mathrm{\ MLE}} = \arg\max_{\boldsymbol{\theta}\in\boldsymbol{\Theta}} \left\{\sum_{i=1}^N \log\big( {p}_{XY}^{}\left(X_i, Y_i; \boldsymbol{\theta}\right) \big)\right\}, \end{align} $$

where the parameter space

(3.6) $$ \begin{align} \boldsymbol{\Theta} = \left\{ \left(\rho, \mu, \sigma^2, \boldsymbol{\tau}^\top\right)^\top\ \Big|\ \rho\in (-1,1),\ \mu\in\mathbb{R},\ \sigma> 0,\ -\infty < \tau_1 < \cdots < \tau_{r-1} < +\infty \right\} \end{align} $$

is the set of legal parameters $\boldsymbol {\theta }$ the MLE maximizes over. As such, $\boldsymbol {\Theta }$ rules out degenerate cases, such as $\rho = \pm 1$ , non-positive standard deviations (SDs), or non-monotonic thresholds. Assuming that the polyserial model is correctly specified for the data at hand, Cox (Reference Cox1974) and Olsson et al. (Reference Olsson, Drasgow and Dorans1982) show that the MLE is consistent for the true $\boldsymbol {\theta }_{\ast}$ , asymptotically normally distributed, and fully efficient.

As a computationally attractive alternative to ML estimation, Olsson et al. (Reference Olsson, Drasgow and Dorans1982) propose a two-step (TS) estimation procedure. In the first step, one computes as estimators of $\mu $ , $\sigma ^2$ , and $\tau _k,k=1,\dots ,r-1,$ the sample statistics

(3.7) $$ \begin{align} \widehat{\mu}_{\textrm{TS}} = \frac{1}{N}\sum_{i=1}^N X_i, \quad \widehat{\sigma}^2_{\textrm{TS}} = \frac{1}{N-1} \sum_{i=1}^N \left(X_i - \widehat{\mu}_{\textrm{TS}}\right)^2, \quad \widehat{\tau}_{k,\textrm{TS}} = \Phi^{-1}\left(\frac{1}{N}\sum_{j=1}^k N_{Y,j} \right), \end{align} $$

respectively, where denotes the empirical marginal frequency of the j-th response option of Y. In the second step, one substitutes for these sample statistics in the log-likelihood in (3.5) and maximizes the ensuing log-likelihood with respect to the remaining parameter, $\rho $ , via conditional ML (conditional on the sample statistics). The main advantage of the TS approach is reduced computing time because one needs to numerically solve the maximization problem (3.5) only with respect to one parameter, $\rho $ , rather than all d parameters in $\boldsymbol {\theta }$ . A drawback is diminished efficiency as compared to (joint) ML, where all d parameters are estimated simultaneously. By means of simulation experiments, Olsson et al. (Reference Olsson, Drasgow and Dorans1982) find that if the polyserial model is correctly specified, ML and the TS estimator produce similar results, with differences decreasing with an increasing number of response options for Y. The TS approach is the default estimator for polyserial correlation in the SEM software package lavaan (Rosseel, Reference Rosseel2012).

While it is in principle possible to robustify the TS approach against contamination by using the same estimation theory as in this article, doing so would result in substantial theoretical and computational drawbacks. We discuss this in detail in Section C.4 of the Supplementary Material.

Alternative estimators of polyserial correlation have been proposed by Bedrick and Breslin (Reference Bedrick and Breslin1996), Lord (Reference Lord1963), and Brogden (Reference Brogden1949), with Bedrick (Reference Bedrick1990, Reference Bedrick1992), Koopman (Reference Koopman1983) and Kraemer (Reference Kraemer1981) studying the theoretical properties of the latter two. We discuss these approaches in more detail in Section 4.2, after conceptualizing misspecification of the polyserial model.

4.1 Partial misspecification of the polyserial model

The polyserial model is misspecified if at least one observation for the continuous-ordinal variable pair $(X,Y)$ has not been generated by the bivariate partially-latent normality model in (3.2). Akin to Welz et al. (Reference Welz, Mair and Alfons2026b), we consider a partial misspecification framework where only a fraction $(1-\varepsilon )$ of observations in a given sample are generated by an underlying normal distribution $P_{X\eta }$ with true parameter $\boldsymbol {\theta }_{\ast}$ , whereas a fixed but unknown fraction $\varepsilon $ of observations are generated from some different but unspecified underlying distribution $H_{X\eta }$ .Footnote ³ Since $H_{X\eta }$ is unspecified, its correlation structure may differ from $P_{X\eta }$ so that observations $(X,Y)$ generated by the underlying $H_{X\eta }$ (after discretization) may be uninformative for the true polyserial correlation coefficient.

Formally, the polyserial model is said to be partially misspecified if the unknown sampling distribution of the observed-latent variable pair $(X,\eta )$ is given by

(4.1) $$ \begin{align} F_{\varepsilon,X\eta} \left( x,v\right) = (1-\varepsilon) {P}_{X\eta}^{}\left(x,v; \boldsymbol{\theta}_{\ast}\right) + \varepsilon H_{X\eta}\left(x,v\right) \end{align} $$

for $x,v\in \mathbb {R}$ . Correspondingly, the implied unknown sampling distribution of the observed variables $(X,Y)$ reads

(4.2) $$ \begin{align} F_{\varepsilon,XY} \left( x,y\right) = (1-\varepsilon) {P}_{XY}^{}\left(x,y; \boldsymbol{\theta}_{\ast}\right) + \varepsilon H_{XY}\left(x,y\right) \end{align} $$

for $x\in \mathbb {R}, y\in \mathcal {Y}$ , where

$$\begin{align*}H_{XY}\left(x,y\right) = \int_{-\infty}^{\tau_{\varepsilon,y}} H_{X\eta}\left(x,\mathrm{d} v\right) \end{align*}$$

is the distribution of the observations for which the polyserial model is (partially) misspecified. The unknown and unspecified discretization thresholds $-\infty = \tau _{\varepsilon ,0} < \tau _{\varepsilon ,1}<\dots < \tau _{\varepsilon ,r-1}< \tau _{\varepsilon ,r} = +\infty $ need not equal the true discretization thresholds $\tau _{{\ast},1} < \dots < \tau _{{\ast},r-1}$ of the polyserial model. Furthermore, denote by $f_{\varepsilon ,XY}$ the unknown density corresponding to the sampling distribution $F_{\varepsilon ,XY}$ .

Misspecification models of the type in (4.1) are standard in the robust statistics literature, where they are known as Huber contamination models, owing to pioneering work of Huber (Reference Huber1964). Following Welz et al. (Reference Welz, Mair and Alfons2026b), we therefore adopt terminology from robust statistics and call $\varepsilon $ the contamination fraction, the uninformative $H_{X\eta }$ the contamination distribution (or simply contamination), and $F_{\varepsilon ,X\eta }$ the contaminated distribution. In the case of a zero-valued contamination fraction ( $\varepsilon = 0$ ), there is no misspecification so that the polyserial model is correctly specified. Indeed, if $\varepsilon = 0$ , then $F_{0,X\eta }(\cdot ,\cdot ) = {P}_{X\eta }^{}\left (\cdot ,\cdot; \boldsymbol {\theta }_{\ast}\right )$ and $F_{0,XY}(\cdot ,\cdot ) = {P}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }_{\ast}\right )$ , so the partial misspecification framework in (4.1) nests the correctly specified polyserial model.

Neither the contamination fraction $\varepsilon $ nor the contamination distribution $H_{X\eta }$ in (4.1) is assumed to be known. Thus, both quantities are left completely unspecified in practice, which “means that we are not making any assumptions on the degree, magnitude, or type of contamination (which is possibly absent altogether)” (Welz et al., Reference Welz, Mair and Alfons2026b). Hence, the polyserial model may be misspecified due to an unlimited variety of reasons, for instance, but not limited to outliers or nonnormality in the continuous X, and/or careless responding or item misunderstanding in the ordinal Y. Because $\varepsilon $ and $H_{X\eta }$ are unspecified in practice, the polyserial model distribution ${P}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }_{\ast}\right )$ remains the distribution of interest: Our only aim is to estimate the parameter $\boldsymbol {\theta }_{\ast}$ of the polyserial model while reducing the impact of potential contamination in the observed data. Consequently, the contaminated distributions $F_{\varepsilon ,XY}$ and $F_{\varepsilon ,X\eta }$ are never estimated. They are purely theoretical objects to study the properties of estimators of the polyserial model when that model is partially misspecified due to data contamination.

While no assumption is made on the specific value of the contamination fraction $\varepsilon $ , we impose the restriction $\varepsilon \in [0,0.5)$ , which is common in robust statistics (e.g., Hampel et al., Reference Hampel, Ronchetti, Rousseeuw and Stahel1986, p. 67). This restriction ensures proper identification: If half or more of the data points were contaminated, it would no longer be possible to distinguish between contamination and observations generated by the polyserial model, at least not without further assumptions. Under such additional assumptions, also values of $\varepsilon \geq 0.5$ could be considered. We refer to Welz et al. (Reference Welz, Mair and Alfons2026b) for a more detailed discussion.

Figure 1 visualizes a simulated example data set generated by the contaminated distribution $F_{\varepsilon ,X\eta }$ in (4.1) with contamination fraction $\varepsilon = 0.15$ . In this example, the contamination distribution $H_{X\eta }$ , whose random draws are orange dots, is a shifted bivariate t-distribution with noncentrality parameter set to $(10,-2)^\top $ , scale matrix $\mathrm {diag}(0.25, 0.25)$ , and 10 degrees of freedom. The remaining realizations (gray dots) are generated by the bivariate normal distribution $P_{X\eta }$ with true parameters $\mu _{\ast} = 0$ , $\sigma _{\ast}^2 = 1$ , and $\rho _{\ast} = 0.5$ . Here, contamination manifests through somewhat larger values in the X-dimension and inflation of the first two response options in the Y-dimension (after discretization), resulting in a (partially) misspecified polyserial model. Applying the MLE in (3.5) to the observed data for $(X,Y)$ in Figure 1 yields a sign-flipped polyserial correlation estimate of $-0.522$ , representing a substantial bias with respect to the true value of $\rho _{\ast} = 0.5$ . In contrast, the robust estimator—which will be introduced in Section 5 and uses the exact same information as ML—estimates a correlation of $0.498$ , which is accurate for the true value of $\rho _{\ast} = 0.5$ .

Figure 1 Simulated data where the polyserial model is misspecified for a fraction $\varepsilon = 0.15$ of the $N=10,000$ points. The gray dots are draws of $(X,\eta )$ from the polyserial model with true parameters $\rho = 0.5, \mu = 0, \sigma ^2 = 1$ , while the orange dots are draws from a contamination distribution $H_{X\eta }$ , being a bivariate t-distribution here with noncentrality parameter $(10,-2)^\top $ , scale matrix $\mathrm {diag}(0.25, 0.25)$ , and 10 degrees of freedom. The horizontal lines mark the thresholds that discretize the latent $\eta $ to the observed ordinal Y with five response options. The numbers in parentheses indicate the population marginal probability of the respective response option under the true polyserial model.

We stress that there exist nonnormal joint distributions of $(X,\eta )$ that, after discretizing the latent $\eta $ variable with fixed thresholds, result in the same density for the observed $(X,Y)$ as if the partially-latent $(X,\eta )$ were jointly normal. This phenomenon is known as discretize equivalence of latent variables (Foldnes & Grønneberg, Reference Foldnes and Grønneberg2019). It follows that there exist contamination distributions $H_{X\eta }$ and contamination fractions $\varepsilon> 0$ in (4.1) under which the contaminated density $f_{\varepsilon ,XY}$ of the observed $(X,Y)$ is exactly equal to the polyserial model density in (3.3) at the true parameter, that is, $f_{\varepsilon ,XY} (x,y) = {p}_{XY}^{}\left (x,y; \boldsymbol {\theta }_{\ast}\right )$ for all $(x,y)\in \mathbb {R}\times \mathcal {Y}$ . Hence, in this scenario, the polyserial model is misspecified, but the misspecification does not have adverse effects because the true polyserial model density of $(X,Y)$ remains unaffected. To avoid cumbersome notation in our analysis, we follow Welz et al. (Reference Welz, Mair and Alfons2026b) and assume consequential misspecification throughout this article, that is, $f_{\varepsilon ,XY} (x,y) \neq {p}_{XY}^{}\left (x,y; \boldsymbol {\theta }_{\ast}\right )$ for at least one response $y\in \mathcal {Y}$ whenever $\varepsilon> 0$ , given $x\in \mathbb {R}$ . Nevertheless, “it is silently understood that misspecification need not be consequential” (Welz et al., Reference Welz, Mair and Alfons2026b), in which case there is no practical problem and both ML and robust estimator are consistent for the true $\boldsymbol {\theta }_{\ast}$ .

Before we define our robust estimator, we briefly juxtapose the partial misspecification framework to distributional misspecification.

4.2 Distributional misspecification

The polyserial model is said to be distributionally misspecified if it is misspecified for all observations in a sample. This contrasts the partial misspecification framework in (4.1) where the model is only misspecified for a (possibly zero-valued) fraction of the observations. Let $G = G_{X,\eta }$ denote the unknown joint distribution of X and the latent $\eta $ , which, under distributional misspecification, is nonnormal for all data points. In distributional misspecification, the object of interest is the correlation coefficient between X and $\eta $ under the distribution G, denoted $\rho _G = \mathbb {C}\mathrm {or}_{G} \left [X,\ \eta \right ]$ , instead of the correlation coefficient under bivariate normality (which would be the polyserial correlation coefficient). Although our robust estimator is designed for partial misspecification rather than distributional misspecification, it can in certain situations also offer a robustness gain under distributional misspecification. We discuss this in more detail in Section 7.3.

5.1 Maximum likelihood revisited

The observed sample can be uniquely characterized by a particular empirical density function

(5.2)

for $x\in \mathbb {R}, y\in \mathcal {Y}$ , where the indicator function takes value 1 if an event E is true, and 0 otherwise. As such, the empirical density $\widehat {f}_{N}$ only takes nonzero values when evaluated at points in the observed sample $\{(X_i,Y_i)\}_{i=1}^N$ . The MLE $\widehat {\boldsymbol {\theta }}_N^{\mathrm {\ MLE}}$ in (3.5) can be expressed as the minimizer of the KL divergence between the empirical density $\widehat {f}_{N}$ and model density ${p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )$ , which is given by

(5.3) $$ \begin{align} {\mathrm{KL}}\left(\widehat{f}_{N} \ \big|\big|\ {p}_{XY}^{}\left(\cdot,\cdot; \boldsymbol{\theta}\right)\right) = \frac{1}{N}\sum_{i=1}^N\Bigg( \log\left(\widehat{f}_{N}^{}\left(X_i,Y_i\right) \right) - \log\left( {p}_{XY}^{}\left(X_i,Y_i; \boldsymbol{\theta}\right)\right) \Bigg) \end{align} $$

and follows by definition of $\widehat {f}_{N}$ and the KL divergence in (5.1) as well as the convention $0\log (0)=0$ .Footnote ⁵ We refer to White (Reference White1982) for a detailed exposition of the connection between KL divergence and ML.

As alluded to earlier, ML estimation tends to be highly susceptible to contamination in the observed sample. It may therefore be desirable to consider alternatives that are designed to more robust against contamination. It turns out that the idea of minimizing a given divergence between the empirical density $\widehat {f}_{N}$ and the model density ${p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )$ can be exploited to construct contamination-robust estimators by choosing alternative divergences to the KL divergence. This is exactly what minimum power divergence estimation (Basu et al., Reference Basu, Harris, Hjort and Jones1998) does.

5.2 Minimum power divergence

Basu et al. (Reference Basu, Harris, Hjort and Jones1998) propose a class of divergences between two densities that generalizes the KL divergence, but is less affected by data contamination. Suppose $g_1$ and $g_2$ are bivariate densities supported on $\mathbb {R}^2$ or a common subset thereof. For a fixed tuning constant $\alpha> 0$ , the divergence of Basu et al. (Reference Basu, Harris, Hjort and Jones1998) between $g_1$ and $g_2$ is defined by

$$ \begin{align*} {\mathrm{D}_\alpha}\left(g_1 \ \big|\big|\ g_2\right) = \int\int \Bigg( g_2^{1+\alpha}(s,t) - \left(1+\frac{1}{\alpha}\right)g_1(s,t)g_2^\alpha(s,t) + \frac{1}{\alpha}g_1^{1+\alpha}(s,t) \Bigg)\mathrm{d} t\mathrm{d} s, \end{align*} $$

where each integral is taken over the densities’ domain. Note that the divergence $\mathrm {D}_\alpha $ is not defined at $\alpha = 0$ . To overcome this issue, Basu et al. (Reference Basu, Harris, Hjort and Jones1998) define $\mathrm {D}_0$ as the limit of $\mathrm {D}_\alpha $ as $\alpha \downarrow 0$ , that is,

(5.4) $$ \begin{align} {\mathrm{D}_0}\left(g_1 \ \big|\big|\ g_2\right) = \lim_{\alpha\downarrow 0}{\mathrm{D}_\alpha}\left(g_1 \ \big|\big|\ g_2\right) = \int\int g_1(s,t) \log\left( \frac{g_1(s,t)}{g_2(s,t)} \right) \mathrm{d} t \mathrm{d} s = {\mathrm{KL}}\left(g_1 \ \big|\big|\ g_2\right), \end{align} $$

where the second equality follows from the identity $z\mapsto \log z = \lim _{\alpha \downarrow 0} \alpha ^{-1}(z^\alpha - 1)$ and the third equality from the definition of the KL divergence in (5.1). Basu et al. (Reference Basu, Harris, Hjort and Jones1998) call the divergence $\mathrm {D}_\alpha , \alpha \geq 0$ , the DPD. It follows that the DPD generalizes the KL divergence. We explain in a moment the intuition behind the DPD.

As before, if the second dimension of the two densities $g_1, g_2$ corresponds to a discrete random variable, that is, the first variable in $g_j$ is continuous, but the second one is discrete, replace the inner integral in the DPD $\mathrm {D}_\alpha , \alpha \geq 0$ , by a summation over the second variable’s domain.

5.3 Proposed robust estimator

For $\alpha> 0$ , the DPD between the empirical density $\widehat {f}_{N}$ in (5.2) and the polyserial model’s density ${p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )$ reads

(5.5) $$ \begin{align} {\mathrm{D}_\alpha}\left(\widehat{f}_{N} \ \big|\big|\ {p}_{XY}^{}\left(\cdot,\cdot; \boldsymbol{\theta}\right)\right) = \int_{\mathbb{R}} \sum_{y\in\mathcal{Y}}{p}_{XY}^{1+\alpha}\left(x,y; \boldsymbol{\theta}\right)\mathrm{d} x - (1+\alpha^{-1}) \frac{1}{N}\sum_{i=1}^N {p}_{XY}^{\alpha}\left(X_i,Y_i; \boldsymbol{\theta}\right) + \alpha^{-1}, \end{align} $$

which follows from writing out $\widehat {f}_{N}$ . For the choice $\alpha = 0$ , the DPD ${\mathrm {D}_0}\left (\widehat {f}_{N} \ \big |\big |\ {p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )\right )$ reduces to the KL divergence in (5.3) by (5.4).

For a pre-specified tuning constant $\alpha \geq 0$ , our proposed estimator minimizes the divergence $\mathrm {D}_\alpha $ between $\widehat {f}_{N}$ and ${p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )$ with respect to the model parameter $\boldsymbol {\theta }\in \boldsymbol {\Theta }$ . Specifically, our proposed estimator is given by

(5.6) $$ \begin{align} \widehat{\boldsymbol{\theta}}_N = \arg\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\mathrm{D}_\alpha}\left(\widehat{f}_{N} \ \big|\big|\ {p}_{XY}^{}\left(\cdot,\cdot; \boldsymbol{\theta}\right)\right). \end{align} $$

In particular, if $\alpha = 0$ , the estimator $\widehat {\boldsymbol {\theta }}_N$ coincides with the MLE $\widehat {\boldsymbol {\theta }}_N^{\mathrm {\ MLE}}$ (see Section 5.1). Since it minimizes a DPD, estimator $\widehat {\boldsymbol {\theta }}_N$ constitutes a minimum DPD estimator (Basu et al., Reference Basu, Harris, Hjort and Jones1998). As such, a limit theory for $\widehat {\boldsymbol {\theta }}_N$ is readily available from Basu et al. (Reference Basu, Harris, Hjort and Jones1998), who derive the theoretical properties of such estimators for general models. Before we turn to the estimator’s limit theory, we first provide some intuition as to why minimum DPD estimators are more robust than ML whenever $\alpha> 0$ .

An estimator defined as the minimum of a (differentiable) loss function can be equivalently characterized as a root of the loss’ gradient; a characterization known as estimating equation. The estimating equation of the MLE $(\alpha = 0)$ can be written as

$$\begin{align*}\frac{1}{N} \sum_{i=1}^N \boldsymbol{u}_{\boldsymbol{\theta}}\left(X_i,Y_i\right) = \boldsymbol{0}, \end{align*}$$

which is satisfied for $\boldsymbol {\theta } = \widehat {\boldsymbol {\theta }}_N^{\mathrm {\ MLE}}$ , and where the d-vector

$$\begin{align*}\boldsymbol{u}_{\boldsymbol{\theta}}\left(x,y\right) = \frac{\partial }{\partial \boldsymbol{\theta}}\log{p}_{XY}^{}\left(x,y; \boldsymbol{\theta}\right), \qquad x\in\mathbb{R},y\in\mathcal{Y}, \end{align*}$$

denotes the log-likelihood score function. A closed-form expression of $\boldsymbol {u}_{\boldsymbol {\theta }}\left (x,y\right )$ is provided in Section A of the Supplementary Material. Analogously, the estimating equation of the proposed estimator $\widehat {\boldsymbol {\theta }}_N$ in (5.6) can be shown to read

(5.7) $$ \begin{align} \frac{1}{N}\sum_{i=1}^N {p}_{XY}^{\alpha}\left(X_i, Y_i; \boldsymbol{\theta}\right) \boldsymbol{u}_{\boldsymbol{\theta}}\left(X_i, Y_i\right) - \boldsymbol{c}_\alpha(\boldsymbol{\theta}) = \boldsymbol{0}, \end{align} $$

which is satisfied for $\boldsymbol {\theta } = \widehat {\boldsymbol {\theta }}_N$ and where $ \boldsymbol {c}_\alpha (\boldsymbol {\theta }) = \int _{\mathbb {R}}\sum _{y\in \mathcal {Y}}{p}_{XY}^{1+\alpha }\left (x,y; \boldsymbol {\theta }\right )\boldsymbol {u}_{\boldsymbol {\theta }}\left (x,y\right )\mathrm {d} x $ is a correction factor independent of observed data. Observe that for $\alpha = 0$ , the estimating equation (5.7) equals that of the MLE,Footnote ⁶ in which all observations in the sample $\{(X_i, Y_i)\}_{i=1}^N$ are weighted equally. Conversely, when $\alpha> 0$ , the term

$$\begin{align*}\frac{1}{N}\sum_{i=1}^N {p}_{XY}^{\alpha}\left(X_i, Y_i; \boldsymbol{\theta}\right) \boldsymbol{u}_{\boldsymbol{\theta}}\left(X_i, Y_i\right) \end{align*}$$

in estimating equation (5.7) “provides a relative-to-the-model downweighting for [contaminated] Footnote ⁷ observations” (Basu et al., Reference Basu, Harris, Hjort and Jones1998, p. 551) with individual-specific weights

(5.8) $$ \begin{align} \widetilde{w}_{i,\alpha}(\boldsymbol{\theta}) = {p}_{XY}^{\alpha}\left(X_i, Y_i; \boldsymbol{\theta}\right). \end{align} $$

A weight $\widetilde {w}_{i,\alpha }(\boldsymbol {\theta })$ is bounded from below by 0 and takes values close to 0 for observations that the polyserial model cannot fit well, such as contaminated data points. As such, “observations that are widely discrepant to the [polyserial] model will get nearly zero weights” (Basu et al., Reference Basu, Harris, Hjort and Jones1998, pp. 551–552). Larger choices of $\alpha> 0$ lead to more stringent downweighting. For instance, a contaminated observation with an ordinal response $Y_i$ that strongly disagrees with the latent threshold structure implied by the model density $p_{XY}$ would receive a weight close to 0.

If contamination is absent $(\varepsilon = 0)$ , it can be shown that the robust estimator $\widehat {\boldsymbol {\theta }}_N$ and the MLE $\widehat {\boldsymbol {\theta }}_N^{\mathrm {\ MLE}}$ in (3.5) are asymptotically equal to one another for all choices $\alpha \geq 0$ . This property is implied by the property of Fisher consistency that minimum DPD estimators possess (implied by Theorem 1 in Basu et al., Reference Basu, Harris, Hjort and Jones1998, see page 551 therein).

Being expressible as a root of a summation over the sample where each summand is a certain function evaluated at a single observation (Eq. (5.7)), minimum DPD estimators are a special case of M-estimation (see Basu et al., Reference Basu, Harris, Hjort and Jones1998, Section 3.1). However, unlike classic M-estimation as proposed in Huber (Reference Huber1964), minimum DPD estimation does not require the data to admit a numerical interpretation due to operationalizing “outlyingness” through discrepant model fit rather than extreme values. Consequently, it can be directly applied to models for mixed data, such as the polyserial model.

5.4 Rescaling of weights

By construction, for a given tuning constant $\alpha> 0$ , the robust estimator’s raw weights $\widetilde {w}_{i,\alpha }(\boldsymbol {\theta })$ in (5.8) are bounded from below by 0 and from above by some data-independent finite value $M_{\alpha }(\boldsymbol {\theta })$ only depending on $\alpha $ and $\boldsymbol {\theta }$ , which is defined by $M_{\alpha }(\boldsymbol {\theta }) = \sup \{{p}_{XY}^{\alpha }\left (x, y; \boldsymbol {\theta }\right ) : x\in \mathbb {R}, y\in \mathcal {Y}\}$ . The fact that the upper bound $M_{\alpha }(\boldsymbol {\theta })$ varies for different values of $\boldsymbol {\theta }$ and $\alpha $ makes it difficult to compare raw weights across different estimates and specifications. Since the raw weights are intended as a tool for detecting contaminated observations, it would be beneficial to have such a comparability. In the following, we describe a simple novel procedure to rescale the weights to always be bounded from above by 1. Rescale the raw weights by their upper bound $M_{\alpha }(\boldsymbol {\theta })$ to obtain

(5.9) $$ \begin{align} w_{i,\alpha}(\boldsymbol{\theta}) = \widetilde{w}_{i,\alpha}(\boldsymbol{\theta}) \big/ M_{\alpha}(\boldsymbol{\theta}), \qquad i = 1,\dots, N. \end{align} $$

We henceforth refer to the rescaled weights $w_{i,\alpha }(\boldsymbol {\theta })$ simply as weights. The rescaled weights are contained in the interval $[0,1]$ and take values close to 0 for poorly fitting observations and value 1 for perfectly fitting observations. This rescaling is without loss of generality because the estimating equation (5.7) can equivalently be expressed as

$$ \begin{align*} \frac{1}{N}\sum_{i=1}^N w_{i,\alpha}(\boldsymbol{\theta}) \boldsymbol{u}_{\boldsymbol{\theta}}\left(X_i, Y_i\right) - \boldsymbol{c}_\alpha(\boldsymbol{\theta}) / M_{\alpha}(\boldsymbol{\theta}) = \boldsymbol{0}. \end{align*} $$

Section B of the Supplementary Material describes an algorithm for computing the upper bound $M_{\alpha }(\boldsymbol {\theta })$ . The usefulness of the rescaled weights as a tool for identifying contaminated data points will be demonstrated in an empirical application in Section 8, where they successfully detect observations discrepant from the model.

6.1 Estimand

It is instructive to discuss what the proposed minimum DPD estimator $\widehat {\boldsymbol {\theta }}_N$ in (5.6) actually estimates. Recall that $f_{\varepsilon ,XY}$ denotes the density of the sampling distribution $F_{\varepsilon ,XY}$ in (4.2) that the observed sample follows. Just like $F_{\varepsilon ,XY}$ , the population density $f_{\varepsilon ,XY}$ is completely unspecified and unknown in practice, including contamination fraction $\varepsilon $ . For a fixed tuning constant $\alpha \geq 0$ , the estimand $\boldsymbol {\theta }_0$ of estimator $\widehat {\boldsymbol {\theta }}_N$ is the parameter that minimizes the DPD between the population density and the polyserial model density, that is,

$$\begin{align*}\boldsymbol{\theta}_0 = \arg\min_{\boldsymbol{\theta}\in\boldsymbol{\Theta}}{\mathrm{D}_\alpha}\left(f_{\varepsilon,XY} \ \big|\big|\ {p}_{XY}^{}\left(\cdot,\cdot; \boldsymbol{\theta}\right)\right). \end{align*}$$

Observe that ${\mathrm {D}_\alpha }\left (f_{\varepsilon ,XY} \ \big |\big |\ {p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )\right )$ is simply the population analog of the empirical divergence ${\mathrm {D}_\alpha }\left (\widehat {f}_{N} \ \big |\big |\ {p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }\right )\right )$ that is minimized by the estimator $\widehat {\boldsymbol {\theta }}_N$ in (5.6).

In the absence of contamination ( $\varepsilon = 0$ ), one has that $f_{0,XY}(\cdot ,\cdot ) = {p}_{XY}^{}\left (\cdot ,\cdot; \boldsymbol {\theta }_{\ast}\right )$ by (4.2), and it follows from Theorem 1 in Basu et al. (Reference Basu, Harris, Hjort and Jones1998) that $\boldsymbol {\theta }_0 = \boldsymbol {\theta }_{\ast}$ , so that the population divergence is zero-valued. Subsequently, just like the MLE, $\widehat {\boldsymbol {\theta }}_N$ is unbiased for the true $\boldsymbol {\theta }_{\ast}$ (in population) if the model is correctly specified. This property is known as Fisher consistency. We stress that Fisher consistency holds true for all choices of $\alpha \geq 0$ . In other words, if contamination is absent, then the proposed estimator is unbiased (in population) no matter the choice of tuning constant $\alpha $ .

On the other hand, if contamination is present $(\varepsilon> 0)$ , then the estimand $\boldsymbol {\theta }_0$ will generally differ from the true $\boldsymbol {\theta }_{\ast}$ . How much they differ depends on the contamination fraction $\varepsilon $ and contamination type $H_{X\eta }$ in the contaminated distribution (4.1), as well as the choice of tuning constant $\alpha \geq 0$ . Roughly speaking, holding the unknown contaminated distribution $F_{\varepsilon ,X\eta }$ constant, the larger $\alpha $ , the closer $\boldsymbol {\theta }_0$ tends to be to the true $\boldsymbol {\theta }_{\ast}$ .Footnote ⁸ In other words, if the polyserial model is misspecified, larger values of $\alpha $ asymptotically tend to lead to less bias for the true $\boldsymbol {\theta }_{\ast}$ .

6.2 Asymptotic analysis

We are now ready to study the asymptotic properties of the proposed estimator $\widehat {\boldsymbol {\theta }}_N$ in (5.6). Under mild standard regularity conditions, it can be shown that $\widehat {\boldsymbol {\theta }}_N$ is asymptotically consistent for estimand $\boldsymbol {\theta }_0$ , that is,

$$\begin{align*}\widehat{\boldsymbol{\theta}}_N\stackrel{\mathbb{P}}{\longrightarrow}\boldsymbol{\theta}_0, \end{align*}$$

as $N\to \infty $ , and is furthermore asymptotically Gaussian,

$$\begin{align*}\sqrt{N} \left( \widehat{\boldsymbol{\theta}}_N - \boldsymbol{\theta}_0 \right) \stackrel{\mathrm{d}}{\longrightarrow} \mathrm{N}_d\Big(\boldsymbol{0}, \boldsymbol{\Sigma}\left( \boldsymbol{\theta}_0 \right) \Big), \end{align*}$$

as $N\to \infty $ , where “ $\stackrel {\mathbb {P}}{\longrightarrow }$ ” and “ $\stackrel {\mathrm {d}}{\longrightarrow }$ ” denote convergence in probability and distribution, respectively. These two stochastic convergence results are rigorously established in Theorem A.1 in the Supplementary Material. The theorem follows immediately from a more general result in Basu et al. (Reference Basu, Harris, Hjort and Jones1998). The asymptotic covariance matrix $\boldsymbol {\Sigma }\left ( \boldsymbol {\theta }_0 \right )$ has a closed-form sandwich-type construction

(6.1) $$ \begin{align} \boldsymbol{\Sigma}\left( \boldsymbol{\theta}_0 \right) = \boldsymbol{J}\left( \boldsymbol{\theta}_0 \right)^{-1} \boldsymbol{K}\left( \boldsymbol{\theta}_0 \right)\boldsymbol{J}\left( \boldsymbol{\theta}_0 \right)^{-1}, \end{align} $$

where the $d\times d$ full-rank matrices $\boldsymbol {\theta }\mapsto \boldsymbol {J}\left ( \boldsymbol {\theta } \right ), \boldsymbol {K}\left ( \boldsymbol {\theta } \right )$ are defined in Section A of the Supplementary Material. Neither $\boldsymbol {J}\left ( \boldsymbol {\theta } \right )$ nor $\boldsymbol {K}\left ( \boldsymbol {\theta } \right )$ are observed in practice, so neither is $\boldsymbol {\Sigma }\left ( \boldsymbol {\theta }_0 \right )$ , but each of these matrices can be consistently estimated, as we will explain momentarily. Owing to this asymptotic normality result, one can construct standard errors and confidence intervals for the estimand $\boldsymbol {\theta }_0$ . We stress that the estimator remains asymptotically normal even when the polyserial model is misspecified. Similar results for inference in general misspecified models with sandwich-type covariance constructions have been derived by, for example, White (Reference White1982) and Huber (Reference Huber, Cam and Neyman1967) for the MLE, Basu et al. (Reference Basu, Harris, Hjort and Jones1998) for robust minimum DPD estimators (of which our proposed estimator is a special case), and Welz (Reference Welz2024) for a class of robust estimators for general categorical data.

Being a function of the unobserved estimand $\boldsymbol {\theta }_0$ as well as additional unobserved population quantities, the population covariance matrix $\boldsymbol {\Sigma }\left ( \boldsymbol {\theta }_0 \right )$ in (6.1) is also unobserved in practice. Nevertheless, we can consistently estimate it by constructing certain estimators $\widehat {\boldsymbol {J}}_N\left ( \boldsymbol {\theta } \right )$ and $\widehat {\boldsymbol {K}}_N\left ( \boldsymbol {\theta } \right )$ , which are pointwise consistent for $\boldsymbol {J}\left ( \boldsymbol {\theta } \right )$ and $\boldsymbol {K}\left ( \boldsymbol {\theta } \right )$ , respectively, for a given parameter vector $\boldsymbol {\theta }\in \boldsymbol {\Theta }$ . Since both estimators are continuous in $\boldsymbol {\theta }$ , it follows from the continuous mapping theorem that the sample analog of the sandwich-type construction in (6.1)

$$\begin{align*}\widehat{\boldsymbol{\Sigma}}_N\left( \boldsymbol{\theta} \right) = \widehat{\boldsymbol{J}}_N\left( \boldsymbol{\theta} \right)^{-1} \widehat{\boldsymbol{K}}_N\left( \boldsymbol{\theta} \right)\widehat{\boldsymbol{J}}_N\left( \boldsymbol{\theta} \right)^{-1} \end{align*}$$

is pointwise consistent for $\boldsymbol {\Sigma }\left ( \boldsymbol {\theta } \right )$ . Combining this result with the convergence result $\widehat {\boldsymbol {\theta }}_N\stackrel {\mathbb {P}}{\longrightarrow }\boldsymbol {\theta }_0$ , the plug-in estimator $\widehat {\boldsymbol {\Sigma }}_N\left ( \widehat {\boldsymbol {\theta }}_N \right )$ is consistent for the asymptotic covariance matrix $\boldsymbol {\Sigma }\left ( \boldsymbol {\theta }_0 \right )$ . We refer to Section A of the Supplementary Material for the definition of $\,\widehat {\boldsymbol {J}}_N\left ( \boldsymbol {\theta } \right )$ and $\widehat {\boldsymbol {K}}_N\left ( \boldsymbol {\theta } \right )$ as well as details.

6.3 Efficiency

For the tuning constant $\alpha =0$ , which corresponds to the MLE $\widehat {\boldsymbol {\theta }}_N^{\mathrm {\ MLE}}$ , it can be shown (see Section A.5 of the Supplementary Material) that if the polyserial model is correctly specified ( $\varepsilon = 0$ ), the asymptotic covariance matrix in (6.1) reduces to the inverted Fisher information of the polychoric model, $\boldsymbol {I}_{\boldsymbol {\theta }_0}^{-1}$ , where

$$\begin{align*}\boldsymbol{I}_{\boldsymbol{\theta}} = \int_{\mathbb{R}} \sum_{y\in\mathcal{Y}} {p}_{XY}^{}\left(x,y; \boldsymbol{\theta}\right)\boldsymbol{u}_{\boldsymbol{\theta}}\left(x,y\right)\boldsymbol{u}_{\boldsymbol{\theta}}\left(x,y\right)^\top\mathrm{d} x \end{align*}$$

denotes the Fisher information of the polyserial model at $\boldsymbol {\theta }\in \boldsymbol {\Theta }$ . Hence, (6.1) nests the well-known result that the MLE’s asymptotic covariance matrix is equal to the inverted Fisher information matrix. It is furthermore well-known that ML estimation is fully efficient, meaning that no unbiased estimator has a smaller variance (e.g., Theorem 3.10 in Lehmann & Casella, Reference Lehmann and Casella1998).

In contrast, for strictly positive tuning constants $\alpha> 0$ , the asymptotic covariance matrix of the corresponding estimator $\widehat {\boldsymbol {\theta }}_N$ in (6.1) is generally not equal to the inverse Fisher information matrix at the polyserial model. Thus, recalling the full efficiency property of the MLE, our robust estimator is less efficient than the MLE. The efficiency loss is due to downweighting of observations with low probability under the polyserial model, which happens even when the model is correctly specified. Consequently, while the robust estimator and MLE are both consistent and unbiased for the true parameter value as long as contamination is absent, the former has a larger estimation variance.

To quantify the efficiency loss, we calculate at the population level the relative efficiency of our robust estimator compared to the MLE when the polyserial model is correctly specified ( $\varepsilon = 0$ ). Recall that in this zero-contamination case, the estimand $\boldsymbol {\theta }_0$ is equal to the true parameter value $\boldsymbol {\theta }_{\ast}$ for all $\alpha \geq 0$ due to the property of Fisher consistency. For a given true value $\boldsymbol {\theta }_{\ast}$ , the relative efficiency (e.g., Van der Vaart, Reference Van der Vaart1998, Chapter 8.2) of our estimator for the true polyserial correlation $\rho _{\ast}$ with respect to the fully efficient MLE is given by

$$\begin{align*}\frac{\mathbb{V}\mathrm{ar} \left[ \widehat{\rho}_N^{\mathrm{\ MLE}} \right]}{\mathbb{V}\mathrm{ar} \left[ \widehat{\rho}_N \right]} = \left(\boldsymbol{I}_{\boldsymbol{\theta}_{\ast}}^{-1}\right)_{1,1} \Big/ \left(\boldsymbol{\Sigma}\left( \boldsymbol{\theta}_{\ast} \right)\right)_{1,1}, \end{align*}$$

where the operator $(\cdot )_{1,1}$ picks out the top left element of a matrix. The efficiency loss of minimum DPD estimators can also be expressed as a function of $\alpha $ for a given model (see Basu et al., Reference Basu, Harris, Hjort and Jones1998, Section 4.2).

Table 1 lists the relative efficiencies of various choices of the tuning constant $\alpha $ for an ordinal Y with $r=5$ response options and true parameter vector $\boldsymbol {\theta }_{\ast} = \left (\rho _{\ast}, \mu _{\ast}, \sigma _{\ast}^2, \tau _{{\ast},1}, \tau _{{\ast},2}, \tau _{{\ast},3}, \tau _{{\ast},4} \right )^\top = \left (0.5, 0, 1, -1.5, -0.5, 0.5, 1.5 \right )^\top $ . Up to about $\alpha = 0.25$ , the relative efficiency of the robust estimator stays above 90%, and then decreases roughly linearly to slightly below 50% for $\alpha = 1$ . Choices of $\alpha> 1$ , while yielding even more robustness, would result in a comparatively poor relative efficiency of considerably less than 50% and are therefore not considered in this article. The efficiency results for other choices of the true parameter vector $\boldsymbol {\theta }_{\ast}$ are very similar and are therefore deferred to Section C of the Supplementary Material. Our software implementation, discussed in more detail in Section 6.4, provides functionality to compute the relative efficiency at arbitrary user-specified parameter vectors and tuning constants $\alpha $ .

Table 1 Relative efficiency for estimating the polyserial correlation coefficient $\rho _{\ast}$ when Y has $r=5$ response options, for various choices of the tuning constant $\alpha $ at a true parameter vector $\boldsymbol {\theta }_{\ast} = \left (\rho _{\ast}, \mu _{\ast}, \sigma _{\ast}^2, \tau _{{\ast},1}, \tau _{{\ast},2}, \tau _{{\ast},3}, \tau _{{\ast},4} \right )^\top = \left (0.5, 0, 1, -1.5, -0.5, 0.5, 1.5 \right )^\top $

A loss of efficiency is a common property of many robust estimators involving continuous random variables, and, as demonstrated in this section, our proposed estimator is no exception.Footnote ⁹ As such, our estimator is based on a compromise between efficiency and robustness. Borrowing an insurance metaphor from Anscombe (Reference Anscombe1960), we “sacrifice some efficiency at the model, in order to insure against accidents caused by deviations from the model” (Huber & Ronchetti, Reference Huber and Ronchetti2009, p. 5). However, it turns out that the efficiency loss is relatively minor for reasonably small positive tuning constants like $\alpha = 0.1$ (Table 1), which, as the simulation experiments in Section 7 will reveal, suffices to gain substantial robustness over ML estimation.

6.4 Implementation

We provide an open-source implementation of our robust estimator as part of the package robcat (for “ROBust CATegorical data analysis”; Welz et al., Reference Welz, Alfons and Mair2026a) for the statistical programming environment R (R Core Team, 2024). The package is freely available from CRAN (the Comprehensive R Archive Network) at https://CRAN.R-project.org/package=robcat. We used this package for obtaining all numerical results in this article.

In principle, every suitable method for numerical optimization can be used to minimize the robust estimator’s minimization problem in (5.6). In our experience, unconstrained optimization via the BFGS algorithm (e.g., Nocedal & Wright, Reference Nocedal and Wright2006, Section 6.1) works well. Yet, additional stability could be gained from making explicit the constraints in the parameter space $\boldsymbol {\Theta }$ in (3.6), for which standard algorithms for constrained optimization can be used, such as the simplex method of Nelder and Mead (Reference Nelder and Mead1965). In our implementation, we adopt a similar default behavior as the implementation of robust polychoric correlation estimation (Welz et al., Reference Welz, Mair and Alfons2026b) in package robcat. Specifically, the implementation first tries unconstrained optimization with the BFGS algorithm. If instability or numerical nonconvergence are encountered or a constraint is violated, it instead uses the algorithm of Nelder and Mead (Reference Nelder and Mead1965) for constrained optimization.Footnote ¹⁰ Nevertheless, numerous other optimization algorithms are supported, and users can freely choose their preferred option.

As for the choice of tuning constant $\alpha \geq 0$ in the DPD function (5.5), higher values lead to greater robustness but less efficiency, whereas values closer to 0 are less robust but more efficient. We shall see in the simulation experiments in the next section that the choice $\alpha = 0.5$ constitutes a good compromise between robustness against contamination up to about $\varepsilon = 0.3$ while being reasonably efficient with about 76% of the MLE’s efficiency (Table 1). Therefore, $\alpha = 0.5$ is the default choice in our implementation. While this is the default choice, we acknowledge that there might be situations in which choices that lead to more (or less) robustness are more appropriate, so we advise users to decide on $\alpha $ on a case-by-case basis.

Moreover, we do not recommend the TS estimation procedure in (3.7) for robust estimation. Recall from (3.7) that the TS estimators for $\mu _{\ast}$ and $\sigma _{\ast}^2$ are given by the sample mean and sample variance, respectively, of the $X_i, i=1,\dots ,N,$ while the estimators for the thresholds are computed from the empirical cumulative frequencies of the $Y_i,i=1,\dots ,N$ . However, if there is contamination in the sample, such as outliers in the data for X and careless responses in the data for Y, all of these three sample statistics might become heavily biased. This bias is then possibly inherited by the estimate of the correlation coefficient in the second stage. To avoid this problem, our robust estimator estimates all model parameters simultaneously.

7.1 Design

We consider a polyserial model under which the random vector $(X, \eta )$ is jointly normally distributed according to (3.2) with true correlation parameter $\rho _{\ast} = 0.5$ and the true first two moments of the observed X are $\mu _{\ast} = \mathbb {E} \left [ X \right ] = 0$ as well as $\sigma ^2_{\ast} = \mathbb {V}\mathrm {ar} \left [ X \right ] = 1$ . The latent variable $\eta $ also has zero mean and unit variance and governs the observed ordinal Y through the discretization process in (3.1) with true threshold parameters $\tau _{{\ast},1} = -1.5, \tau _{{\ast},2} = -0.5, \tau _{{\ast},3} = 0.5, \mathrm { and } \tau _{{\ast},4} = 1.5,$ so that Y has $r=5$ response categories. Using integer scoring $\mathcal {Y} = \{1,2,\dots ,5\}$ , the corresponding true point polyserial correlation coefficient in this setting amounts to $\widetilde {\rho }_{\ast} = 0.477$ (see Section 3.2).

To simulate contamination, we replace a fraction $\varepsilon $ of the data for $(X, \eta )$ with draws from a particular contamination distribution $H_{X\eta }$ , which is not modeled by our robust estimator nor ever assumed to be known. Specifically, the contamination distribution $H_{X\eta }$ in this simulation is set to a shifted bivariate t-distribution with noncentrality parameter $(10,-2)^\top $ , scale matrix $\mathrm {diag}(0.25, 0.25)$ , and 10 degrees of freedom. To obtain contaminated data for the ordinal Y, we discretize this contamination distribution’s realizations in the $\eta $ -dimension according to the same thresholds $\tau _{{\ast},1}, \dots , \tau _{{\ast},4}$ as the uncontaminated realizations from the polyserial model. The ensuing contaminated observations are mean-shifted to the right in the continuous X-dimension and primarily inflate the first and second response options in the ordinal Y-dimension. For instance, the contaminated data in Figure 1 were generated by this process (with contamination fraction $\varepsilon = 0.15$ ) and correspondingly exhibit these features. In the robust statistics literature, such contaminating data are an example of negative leverage points because they drag positive correlational estimates toward zero or even negative values, thereby creating negative leverage.

We sample $N=500$ observations $(X_i, Y_i), i=1,\dots , N$ , from this process with contamination fraction $\varepsilon \in \{0, 0.002, 0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.49\}$ . Note that $\varepsilon = 0.002$ corresponds to only one single contaminated data point for the sample size $N=500$ . For each simulated data set, we estimate the true parameter $\boldsymbol {\theta }_{\ast}$ by means of a minimum DPD estimator $\widehat {\boldsymbol {\theta }}_N$ with tuning constants $\alpha \in \{0,0.1,0.25,0.5, 0.75,1\}$ . This procedure is repeated 5,000 times.

Recall that the tuning constant choice $\alpha = 0$ corresponds to the MLE, for which we use the usual Fisher-information-based asymptotic covariance matrix to compute standard errors instead of the sandwich-type construction in (6.1). Further, recall that larger values of $\alpha $ lead to a more robust estimator, at the cost of a loss in efficiency. In fact, the relative efficiencies in Table 1 were computed for the same parameter value that serves as true value $\boldsymbol {\theta }_{\ast}$ in this simulation. We do not consider $\alpha $ -values beyond 1 due to their poor efficiency properties of having substantially less than 50% of the MLE’s efficiency (see Table 1). In addition to an estimate $\widehat {\rho }_N$ of the true polyserial correlation coefficient $\rho _{\ast}$ , we also construct an estimate $\widehat {\widetilde {\rho }}_N$ of the true point polyserial correlation coefficient $\widetilde {\rho }_{\ast}$ from the individual parameter estimates in $\widehat {\boldsymbol {\theta }}_N$ (see Section 3.2 for details).

Let $\widehat {\rho }_N$ be a polyserial correlation estimate computed from a given simulated data set, and $\widehat {\mathrm {SE}}\left (\widehat {\rho }_N\right )$ be its associated standard error estimate constructed from the limit theory in Theorem A.1 in the Supplementary Material. We evaluate performance by means of the following metrics:

• • Bias of the polyserial correlation estimate, $\widehat {\rho }_N - \rho _{\ast}$ , and the point polyserial correlation estimate, $\widehat {\widetilde {\rho }}_N - \widetilde {\rho }_{\ast}$ , averaged over the 5,000 repetitions.

• • Average approximate bias of the standard error estimate, defined as the difference between $\widehat {\mathrm {SE}}\left (\widehat {\rho }_N\right )$ and the sample SD of the individual correlation estimates $\widehat {\rho }_N$ across the 5,000 repetitions, where the latter is a finite-sample approximation of the true population standard error. A rigorous definition is provided in Section C.2 of the Supplementary Material.

• • Coverage at the significance level $\gamma = 0.05$ , defined as the proportion (across repetitions) of confidence intervals $\left [ \widehat {\rho }_N \mp q_{1-\gamma /2} \cdot \widehat {\mathrm {SE}}\left (\widehat {\rho }_N\right ) \right ]$ that contain the true $\rho _{\ast}$ , where $q_{1-\gamma /2}$ denotes the $(1-\gamma /2)$ quantile of the standard normal distribution.

• • Average confidence interval length at significance level $\gamma = 0.05$ , where the length of an individual confidence interval is given by $2\cdot q_{1-\gamma /2} \cdot \widehat {\mathrm {SE}}\left (\widehat {\rho }_N\right )$ .

7.3 Additional simulations

We perform three additional simulation studies, all described in detail in Section D of the Supplementary Material. These simulations are intended to explore the benefits and limitations of our proposed methodology.

The first additional simulation, described in Section D.1 of the Supplementary Material, has the same setup as the design in Section 7.1, but the contamination manifests through extreme outliers in both the X- and $\eta $ -dimension. This simulation is primarily intended to study the behavior of the MLE and robust estimator under gross errors in the data. In brief, one single gross error observation in a sample suffices for the MLE to produce sign-flipped estimates that are extremely biased. In contrast, the robust estimator remains nearly unaffected by small to moderate contamination fractions with such gross errors, and only starts exhibiting considerable biases in high-contamination settings.

The second additional simulation, described in Section D.2 of the Supplementary Material, is motivated by the fact that the simulation design in Section 7.1 is primarily concerned with mean-shifted contamination. Mean-shifted contamination is characterized by the population mean of the contamination distribution $H_{X\eta }$ being markedly different than that of the true normal distribution $P_{X\eta }$ . However, contamination may also manifest through changes in correlation rather than means. The second additional simulation therefore focuses on a correlation-shifted contamination distribution $H_{X\eta }$ , which is equal to the true normal distribution $P_{X\eta }$ , except for a sign-flipped correlation coefficient $-\rho _{\ast}$ . In brief, for moderate polyserial correlation $\rho _{\ast}$ (where $H_{X\eta }$ and $P_{X\eta }$ substantially overlap), our robust estimator does not yield an improvement over ML. Conversely, it does provide a notable gain in robustness for larger true correlations $\rho _{\ast}$ (where $H_{X\eta }$ and $P_{X\eta }$ barely overlap).

Combined with the simulations on mean-shifted contamination in Section 7.1 and Section D.1 of the Supplementary Material, the results of the simulation on correlation-shifted contamination suggest that the potential for robustness gain depends on the overlap between the contamination distribution $H_{X\eta }$ and true normal distribution $P_{X\eta }$ . If there is much overlap, the robust estimator cannot distinguish well between contamination and regular observations since it does not make any assumptions on the former, so it may not improve upon ML. In contrast, if $H_{X\eta }$ and $P_{X\eta }$ do not overlap much—like with mean-shifted contamination or correlation-shifted contamination with strong true correlation—the robust estimator can identify contaminated observations and subsequently downweigh their influence to achieve considerable robustness gains.

The third additional simulation, described in Section D.3 of the Supplementary Material, is concerned with distributional misspecification, where the polyserial model is misspecified for all observations in a sample (cf. Section 4.2). The simulation shows that the potential of robustness gain with our robust estimator depends on the case-specific characteristics of the nonnormal sampling distribution of $(X,\eta )$ for which the polyserial model is distributionally misspecified. In brief, if the sampling distribution can be reasonably well approximated by a mixture of a normal distribution and some other unknown distribution to emulate the contamination model in (4.1), then the robust estimator can improve upon the MLE. If the sampling distribution does not admit such an approximation, the robust estimator may not yield an improvement but produce similar estimates as ML, which can be quite poor under distributional misspecification (e.g., Bedrick, Reference Bedrick1995, and our simulations in Section D.3 of the Supplementary Material). Thus, we advice applied researchers to always test for partially-latent normality when ML and robust estimator yield similar estimates. While such similarity may be due to normality indeed holding true, it may also be due to distributional misspecification from a specific nonnormal distribution for which the robust estimator cannot improve over ML. We provide details and more explicit guidance in Section D.3 of the Supplementary Material.