# Robust Estimation and Inference in Categorical Data††thanks: I thank Alberto
Abadie, Andreas Alfons, Isaiah Andrews, Marco Avella Medina, Patrick Groenen,
Nick Koning, Patrick Mair, Anna Mikusheva, Whitney Newey, Alberto Quaini,
Elvezio Ronchetti, Stefanie Stantcheva, and Ming Yuan, as well as the
participants of the MIT Econometrics lunch seminar and the Econometric
Institute internal seminar for helpful comments and suggestions. Most of this
work was carried out when I was visiting the economics department at MIT; I
thank them for their hospitality. This work was supported by a grant from the
Dutch Research Council (NWO), research program Vidi (Project No.
VI.Vidi.195.141).

Max Welz
<EMAIL_ADDRESS>
(Econometric Institute
Erasmus School of Economics
)

###### Abstract

In empirical science, many variables of interest are categorical. Like any
model, models for categorical responses can be misspecified, leading to
possibly large biases in estimation. One particularly troublesome source of
misspecification is inattentive responding in questionnaires, which is well-
known to jeopardize the validity of structural equation models (SEMs) and
other survey-based analyses. I propose a general estimator that is designed to
be robust to misspecification of models for categorical responses. Unlike
hitherto approaches, the estimator makes no assumption whatsoever on the
degree, magnitude, or type of misspecification. The proposed estimator
generalizes maximum likelihood estimation, is strongly consistent,
asymptotically Gaussian, has the same time complexity as maximum likelihood,
and can be applied to any model for categorical responses. In addition, I
develop a novel test that tests whether a given response can be fitted well by
the assumed model, which allows one to trace back possible sources of
misspecification. I verify the attractive theoretical properties of the
proposed methodology in Monte Carlo experiments, and demonstrate its practical
usefulness in an empirical application on a SEM of personality traits, where I
find compelling evidence for the presence of inattentive responding whose
adverse effects the proposed estimator can withstand, unlike maximum
likelihood.

Keywords: Model misspecification, robust estimation, categorical variables,
inattentive responding

##  1 Introduction

In the social, economic, psychological, and behavioral sciences, many
quantities of interest are measured by means of categorical variables, often
through surveys or questionnaires. After collecting such data, researchers
typically model the quantities of interest by employing models specifically
designed for the categorical measurements thereof, for example structural
equation models for latent character traits or item response models for item
difficulty. However, such models can be misspecified, and it is well-known
that model misspecification leads to inconstant estimation of model
parameters, which may ultimately result in incorrect research findings. In
survey-collected data, a particularly relevant source of model
misspecification is inattentive or careless responding by survey participants
(e.g. Stantcheva,, 2022; Meade & Craig,, 2012; Huang et al., 2015b, ). While
there are many estimation methods that are designed to be robust to model
misspecification (see Maronna et al.,, 2018, for a recent overview), these
methods almost exclusively designed for continuous data, and may fail to be
effective or even computable when the data are categorical.

Motivated by the empirical relevance of this gap in the literature due to the
increasing popularity of surveys in economics (Stantcheva,, 2022) and their
long-standing importance in behavioral research (e.g. Rust et al.,, 2020), I
propose a novel estimation framework for models of categorical data that is
designed to be robust to misspecification of that model. The estimator is very
general in the sense that it can be applied to nearly any model of categorical
data, and, crucially, it makes no assumption whatsoever on the type,
magnitude, or location of potential misspecification. Hence, the proposed
estimator is robust to an unlimited and unspecified variety of possible
sources of model misspecification. I show that the estimator generalizes
commonly employed maximum likelihood estimation and possesses attractive
theoretical and computational properties. For instance, the estimator is
strongly consistent and asymptotically normally distributed, while it is of
the same time complexity as maximum likelihood, meaning that it comes at no
additional computational cost. Furthermore, I develop a novel statistical test
that test if a given categorical data point can be fitted well by the presumed
model. The test rejects this null hypothesis if the data point in question
cannot be fitted sufficiently well, and thereby helps pinpoint potential
sources of model misspecification.

As such, the proposed methodology is particularly attractive to fit models for
survey responses such as structural equation models or item response models
because it is likely that a certain proportion of survey participants responds
inattentively or carelessly (Ward & Meade,, 2023). The estimator can not only
withstand the harmful effects of inattentive respondents, but also identify
them through the diagnostic test.

I verify the attractive theoretical and robustness properties of the proposed
estimator by means of extensive simulation studies and demonstrate its
practical usefulness in an empirical application on a structural equation
model on a measurement of the Big-5 personality traits (Goldberg,, 1992). I
find compelling evidence for the presence of inattentive respondents that, if
unaccounted for, have a sizable effect on parameter estimates. For instance,
the correlation coefficient between two mutually contradic·tory items in a
neuroticism scale is estimated as $-0.62$ by hitherto estimation methods,
whereas my estimator yields a substantially stronger correlational estimate of
$-0.93$, which is, unlike the former estimate, in line with literature on this
scale. Likewise, a structural equation model fitted by my robust estimator
yields stronger and internally more consistent factor loadings that can
explain substantially more variation than a fit based on commonly employed
methods. In addition, I demonstrate how the estimator can be used to robustly
estimate scale reliability coefficients like Cronbach’s $\alpha$. As such, the
entire process of building, testing, and verifying scales can be robustified
against inattentive responding without making any assumptions on the nature of
inattentive responding, a departure from previous literature.

With its focus on models for questionnaire responses, this paper ties into a
growing literature concerned with the validity of research findings when
employed models are misspecified due to responses that do not follow the
assumed model and subsequently cannot be fitted well by that model. An ensuing
poor model fit may lead a researcher to doubt or even reject the theoretical
model (Lai & Green,, 2016). Poor model fit can occur for two reasons. First,
the theory behind the model may simply be wrong, in case of which the model is
also wrong and the theory is correctly rejected by the data. Second, theory
and model may in fact be correct, yet a misfit occurs due to the presence of a
limited number of observations from a different population, such as
inattentive respondents or heterogeneous subgroups, which may lead one to
incorrectly reject a theory (Arias et al.,, 2020). Indeed, already a small
proportion of inattentive responses can substantially deteriorate model fit
(Arias et al.,, 2020; Huang et al., 2015a, ; Woods,, 2006), and ultimately
lead a researcher to reject a correct hypothesis or sustain an incorrect
hypothesis (Arias et al.,, 2020; Huang et al., 2015b, ; Maniaci & Rogge,,
2014; McGrath et al.,, 2010; Woods,, 2006; Schmitt & Stults,, 1985).
Inattentive responding itself is widely prevalent (Ward & Meade,, 2023;
Bowling et al.,, 2016; Meade & Craig,, 2012), with most estimates on its
prevalence ranging from 10–15% of study participants (Curran,, 2016; Huang et
al., 2015b, ; Huang et al.,, 2012; Meade & Craig,, 2012), while already a
prevalence 5–10% can jeopardize the validity of research findings (Arias et
al.,, 2020; Credé,, 2010; Woods,, 2006). In fact, Ward & Meade, (2023)
conjecture that inattentive responding is likely present in all survey data.
Due to the damaging effects of inattentive responses, a large number of
methods for their detection has emerged, ranging from consistency indicators
such as psychometric antonyms/synonyms (Meade & Craig,, 2012) over response
times (e.g. Bowling et al.,, 2023) to model-based techniques, such as person-
fit statistics (e.g. Drasgow et al.,, 1985), structural equation models (e.g.
Kim et al.,, 2018), mixture models (e.g. Arias et al.,, 2020), or attention
check items (e.g., Section 3 in Stantcheva,, 2022). More recently, machine
learning techniques have been proposed (Welz & Alfons,, 2023; Schroeders et
al.,, 2022). I refer to Alfons & Welz, (2024) for a recent review of
inattention in survey data. I deviate from hitherto methods by not making any
assumptions on the nature of potential misspecification and deriving
statistical guarantees on the performance of my approach.

In econometrics, this paper ties into a recent stream of methodological
literature on misspecification-robust estimation (e.g. Andrews et al.,, 2017;
Armstrong & Kolesár,, 2021; Bugni & Ura,, 2019; Bonhomme & Weidner,, 2022;
Chernozhukov et al.,, 2018, 2022; Kitamura et al.,, 2013), and develops a
categorical analogue to robust $M$-estimation (e.g. Van de Geer,, 2000), which
is primarily intended for models of continuous data.

This paper is organizes as follows. Section 2 describes the modeling framework
and Section 3 the proposed methodology, while Sections 4 and 5 contain Monte
Carlo experiments and an empirical application, respectively. Section 6
concludes.

I employ the following notation throughout. A vector $\bm{a}\in\mathbb{R}^{d}$
is understood as a $(d\times 1)$ column vector, and $\bm{a}^{\top}$ denotes
its $(1\times d)$ transpose. Further, denote by
$\|\bm{a}\|=\sqrt{\bm{a}^{\top}\bm{a}}$ its Euclidean norm, and, for a square
matrix $\bm{A}\in\mathbb{R}^{d\times d}$, denote by
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{A}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}=\sqrt{\lambda_{\max}(\bm{A}^{\top}\bm{A})}$
its operator norm, where $\lambda_{\max}(\cdot)$ is the largest eigenvalue of
a square matrix. Assuming that $\bm{A}$ is symmetric positive-definite, let
$\bm{A}^{1/2}$ denote a lower triangular $(d\times d)$ matrix that satisfies
the Cholesky factorization
$\bm{A}=\bm{A}^{1/2}\left(\bm{A}^{1/2}\right)^{\top}$. Furthermore, let
$\Phi_{d}$ denote the distribution function of the $d$-variate standard normal
distribution, and let $\mathbb{P}$ and $\mathbb{E}$ be a generic probability
measure and a generic expectation operator, respectively. Throughout this
paper, I make the implicit assumption that all stochastic objects are well-
defined on an appropriately chosen measure space. Furthermore, let the symbols
“$\stackrel{{\scriptstyle\mathrm{d}}}{{\longrightarrow}}$” and
“$\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}$” denote
convergence in distribution and almost sure convergence, respectively.

##  2 Setup and modeling framework

This section introduces the modeling framework that will be employed
throughout the paper. The considered framework is very general and many well-
known models emerge as special cases.

###  2.1 Notation and examples

Suppose one observes a $k$-dimensional categorical random variable
$\bm{Z}=(Z_{1},Z_{2},\dots,Z_{k})^{\top}$ that takes values in a finite sample
space $\bm{\mathcal{Z}}=\\{\bm{z}_{1},\bm{z}_{2},\dots,\bm{z}_{m}\\}$ with $m$
possible vector-valued outcomes of dimension $k$. Hence, realizations of this
random variable can be represented in a $k$-way contingency table that cross-
tabulates empirical frequencies for each outcome. The individual categories in
$\bm{\mathcal{Z}}$ can either be ordered (so that $\bm{Z}$ is
ordinal/quantitative, like educational attainment) or unordered (so that
$\bm{Z}$ is nominal/qualitative, like marital status).

Models for categorical outcomes typically parametrize the probability mass
function (henceforth density) of each of the outcomes, and the primary
statistical problem is to estimate the model’s parameters from observed data.
Formally, a model
$\big{\\{}\bm{p}(\bm{\theta})=(p_{\bm{z}}\left(\bm{\theta}\right))_{\bm{z}\in\bm{\mathcal{Z}}}:\bm{\theta}\in\bm{\Theta}\big{\\}}$
subject to $d$ parameters and parameter space
$\bm{\Theta}\subset\mathbb{R}^{d}$ assigns to each categorical outcome
$\bm{z}=(z_{1},z_{2},\dots,z_{k})^{\top}\in\bm{\mathcal{Z}}$ a nonnegative
probability

$p_{\bm{z}}\left(\bm{\theta}\right)=\mathbb{P}_{\bm{\theta}}\left[\bm{Z}=\bm{z}\right]=\mathbb{P}_{\bm{\theta}}\left[Z_{1}=z_{1},\dots,Z_{k}=z_{k}\right]\quad\text{such
that}\quad\sum_{\bm{z}\in\bm{\mathcal{Z}}}p_{\bm{z}}\left(\bm{\theta}\right)=1,$
(1)

which depends on a constant but unobserved parameter
$\bm{\theta}\in\bm{\Theta}$ that one wishes to estimate by using realizations
of categorical $\bm{Z}$. Such models can be used for predicting test results
(like the Rasch model in item response theory), modeling the association
between latent personality traits, categorical regression, or predicting how
often a certain event occurs, just to name a few. I present a selection of
such models in Appendix A. Textbook treatments can be found in Agresti, (2010)
for ordinal variables and Agresti, (2012) for nominal ones. To fix ideas, the
following section focuses in more detail on a particularly relevant type of
categorical model, namely structural equation modeling of questionnaire
responses, the harmful effects that misspecification of such models can
entail, and how issues stemming from misspecification can be alleviated with
the methodology proposed in this paper.

###  2.2 Lead example: Structural Equation Modeling

Structural equation models (SEMs) are a fundamental tool to analyze
questionnaire responses in the social, behavioral, and business sciences. SEMs
are typically based on factor analyses. Factor models are designed to explain
the correlation between a set of modeled variables in terms of a small number
of latent variables called factors (Mardia et al.,, 1979, Ch. 9). Such models
are attractive for questionnaire data because questionnaires are typically
designed to measure latent constructs such as attitudes or personality traits
by means of multiple rating items for each construct (repeated elicitation).
Multiple items measuring the same construct are called a scale. Formally, a
factor model for $r$ factors (like constructs) modeling $q$ variables (like
rating responses) is in its most basic form given by

$\bm{\Sigma}=\bm{\Lambda}\bm{\Lambda}^{\top}+\bm{\Psi},$ (2)

where $\bm{\Sigma}$ is the population correlation matrix of the $q$ modeled
variables,111One can also use covariance matrices for factor models, but it is
generally recommended to use correlation matrices to account for scale
differences in the measurement of the modeled variables (e.g. Mardia et al.,,
1979, Ch. 9). $\bm{\Lambda}$ is a $q\times r$ matrix of coefficients (factor
loadings) that summarize the variation explained by the factors, and
$\bm{\Psi}$ is a diagonal covariance matrix capturing unexplained variation.
Factor models like (2) are primarily used to develop, test, and validate
behavioral theory (e.g. Rust et al.,, 2020), are frequently utilized in
economics (e.g. Almlund et al.,, 2011; Borghans et al.,, 2008; Heckman et
al.,, 2006; Osborne-Groves,, 2005), and offer many econometric extensions (see
Ch. 3 in Schennach,, 2022, for an overview).

In practice, one first estimates population correlation matrix $\bm{\Sigma}$
from some sample of the $q$ modeled variables and then fits factor model (2)
to the estimated correlation matrix by estimating $\bm{\Lambda}$ and
$\bm{\Psi}$, subject to some identification constraints. However, hitherto
estimation methods like maximum likelihood are highly susceptible to
misspecification of factor model (2), resulting in inconsistent estimates of
$\bm{\Lambda}$ and $\bm{\Psi}$, which in turn can lead to incorrect research
findings. In the context of factor models for questionnaire responses, a
prominent source of misspecification are inattentive or careless responding,
which have been shown to be a serious threat for the validity of SEM analyses
due to a variety of psychometric issues, such as, for example but not limited
to, reduced scale reliability (Arias et al.,, 2020), reduced construct
validity (Kam & Meyer,, 2015), improper factor structure (Arias et al.,, 2020;
Huang et al., 2015b, ; Woods,, 2006), as well as errors in hypothesis testing
(Maniaci & Rogge,, 2014; McGrath et al.,, 2010; Woods,, 2006), and there is
evidence that already a small proportion of inattentive respondents of 5–10%
are problematic (Arias et al.,, 2020; Credé,, 2010; Schmitt & Stults,, 1985;
Woods,, 2006).

Since factor analyses are fundamentally a model of the correlation structure
of the modeled variables (see eq. 2), one approach to robustify factor
analyses against inattentive responding (or other sources of misspecification)
is to robustify the estimation of the population correlation matrix
$\bm{\Sigma}$ that is being modeled.222For factor models of continuous
observed variables, Pison et al., (2003) propose to robustly estimate the
correlation matrix by means of the outlier-robust Minimum Covariance
Determinant estimator (Rousseeuw,, 1985). However, this estimator crucially
relies on continuity and Gaussianity of the modeled variables, and is
therefore not applicable to categorical data (and may not even be computable
for such). The commonly employed Pearson sample correlation matrix of rating
items (where each response option is encoded by an integer) is strongly
affected by inattentive responses (Credé,, 2010; Raymaekers & Rousseeuw,,
2021), so a different correlation estimator is needed. An often recommended
alternative to sample correlation between rating items is polychoric
correlation (e.g. Foldnes & Grønneberg,, 2022; Garrido et al.,, 2013; Holgado-
Tello et al.,, 2010), but polychoric correlation estimated by maximum
likelihood is also sensitive to inattentive responding (cf. Foldnes &
Grønneberg,, 2022) and generally yields similar estimates as sample
correlation when there are five or more rating options (Rhemtulla et al.,,
2012). However, it turns out that estimating the correlation matrix
$\bm{\Sigma}$ by means of polychoric correlation coefficients is robust to
inattentive responding if the latter are estimated with the robust estimator
proposed in this paper. It follows that one can make factor analyses on
questionnaire responses robust to inattentive responding by fitting the factor
model to a polychoric correlation matrix estimated by my robust estimator.
Consequently, subsequent analyses that rely on estimated correlation or factor
structure such as scale reliability (e.g. coefficients $\alpha$ (Cronbach,,
1951) or $\omega$ (McDonald,, 1970)) can also be robustified against
inattention with this approach.

In the pairwise polychoric correlation model (Pearson & Pearson,, 1922;
Pearson,, 1901), one observes a bivariate variable $\bm{Z}=(X,Y)^{\top}$ of
ratings that takes values in the set
$\bm{\mathcal{Z}}=\mathcal{X}\times\mathcal{Y}=\\{1,2,\dots,J_{x}\\}\times\\{1,2,\dots,J_{y}\\}$,
that is, the items associated with $X$ and $Y$ have $J_{x}$ and $J_{y}$
response categories, respectively. The polychoric model assumes that the
probability of observing responses
$\bm{z}=(x,y)\in\mathcal{X}\times\mathcal{Y}$ is given by

$p_{\bm{z}}\left(\bm{\theta}\right)=\mathbb{P}_{\bm{\theta}}\left[X=x,Y=y\right]=\int_{a_{x-1}}^{a_{x}}\int_{b_{y-1}}^{b_{y}}\phi_{2}(t,s;\varrho)\
\textnormal{d}s\ \textnormal{d}t,$ (3)

where $\phi_{2}(\cdot,\cdot;\varrho)$ is the bivariate standard normal density
with correlation parameter $\varrho\in(-1,1)$, and real-valued threshold
parameters $a_{1}<\dots<a_{J_{x}-1}$, $b_{1}<\dots<b_{J_{y}-1}$ with the
convention $a_{0}=b_{0}=-\infty,a_{J_{x}}=b_{J_{y}}=+\infty$. This model is
parametrized by a vector
$\bm{\theta}=(\varrho,a_{1},\dots,a_{J_{x}-1},b_{a},\dots,b_{J_{y}-1})^{\top}$
of dimension $d=J_{x}+J_{y}-1$, where the object of primary interest is the
polychoric correlation coefficient $\varrho$, which is the correlation of the
two latent variables that $X$ and $Y$ measure. For parameter estimation in the
polychoric model, maximum likelihood is typically employed (Olsson,, 1979). If
a given dataset comprises rating responses to $q$ items, the associated
$q\times q$ polychoric correlation matrix of all items holds the $q(q-1)/2$
unique pairwise polychoric correlation coefficients $\varrho_{ij}$ for
distinct items $i,j=1,\dots,q$. Guidelines for polychoric correlation matrices
can be found in Gadermann et al., (2012).

###  2.3 Conceptualization of model misspecification

In the lead example of SEMs, the primarily discussed source of model
misspecification are inattentive respondents in questionnaire studies. Yet,
misspecifiation of models for categorical outcomes (not necessarily restricted
to questionnaire responses) may manifest in numerous ways. In this section, I
conceptualize general model misspecification by adopting the seminal work of
Huber, (1964).

Suppose that there exists a “true” but unknown parameter
$\bm{\theta}_{*}\in\bm{\Theta}$ whose associated density a categorical random
variable of interest $\bm{Z}$ follows, that is,
$\bm{Z}\sim\bm{p}\left(\bm{\theta}\right)=(p_{\bm{z}}\left(\bm{\theta}_{*}\right))_{\bm{z}\in\bm{\mathcal{Z}}}$.
Subsequently, the primary statistical problem is to estimate the true
$\bm{\theta}_{*}$ based on a sample of realizations of $\bm{Z}$. If all data
points in the sample are indeed generated by the true density
$\bm{p}\left(\bm{\theta}_{*}\right)$, then one can generally construct a
consistent estimator of $\bm{\theta}_{*}$, for instance through maximum
likelihood estimation. However, if some observed data points are not generated
by $\bm{p}\left(\bm{\theta}_{*}\right)$ but some other distribution outside of
model $\\{\bm{p}\left(\bm{\theta}\right):\bm{\theta}\in\bm{\Theta}\\}$, then
one says that the model is misspecified.

Without loss of generality, potential model misspecification can be
conceptualized by means of a mixture density

$\bm{\mathcal{Z}}\ni\bm{z}\mapsto
f_{\varepsilon}(\bm{z})=(1-\varepsilon)p_{\bm{z}}\left(\bm{\theta}_{*}\right)+\varepsilon
h(\bm{z})$ (4)

for misspecification degree $\varepsilon\in[0,1]$ and misspecification type
$h(\cdot)$, which is some density on $\bm{\mathcal{Z}}$ that is not contained
in model $\\{\bm{p}\left(\bm{\theta}\right):\bm{\theta}\in\bm{\Theta}\\}$.
Hence, one samples from the model distribution
$p_{\bm{z}}\left(\bm{\theta}_{*}\right)$ with probability $1-\varepsilon$ and
from some other density $h(\bm{z})$ with probability $\varepsilon$. Neither
$\varepsilon$ nor $h(\cdot)$ are known and are left completely unspecified so
that one makes no assumption whatsoever on the degree, magnitude, or type of
misspecification, which might be absent altogether. Indeed, if
$\varepsilon=0$, the model is correctly specified because the true model
probability and population probability are equal to one another, that is,
$f_{0}(\bm{z})=p_{\bm{z}}\left(\bm{\theta}_{*}\right)$ for all
$\bm{z}\in\bm{\mathcal{Z}}$.

In the context of models for questionnaire responses (like SEMs), leaving
misspecification degree $\varepsilon$ and type $h(\cdot)$ in (4) unspecified
means that the model can be misspecified due to an unlimited variety of
reasons, for instance but not limited to inattentive or careless responding
(e.g., straightlining, pattern responding, random-like responding),
misresponses, item misunderstanding, or accurate responses that are simply not
generated by the assumed model (sampling errors).

Conceptualizing model misspecification in similar fashion as (4) is commonly
done in the robust statistics literature (see Maronna et al.,, 2018, for a
recent overview) and econometric literature on misspecified moment conditions
(e.g. Andrews et al.,, 2017; Armstrong & Kolesár,, 2021; Bugni & Ura,, 2019;
Bonhomme & Weidner,, 2022; Chernozhukov et al.,, 2018, 2022; Kitamura et al.,,
2013; Newey,, 1985).

###  2.4 Why model in the first place?

It is worthwhile to briefly address the need for modeling in categorical data.
First, minimizing the Kullback-Leibler divergence between the empirical
relative frequency of the individual outcomes and some non-parametrized
theoretical outcome probabilities with respect to the theoretical outcome
probabilities yields as estimator the empirical relative frequencies of the
outcomes, which is uninformative. Assuming a model that imposes a
parametrization based on theory is standard in many fields; for instance,
factor models are commonly employed in the analysis of questionnaire responses
with repeated elicitation, which aids in validating, rejecting, or developing
novel theories on human behavior (e.g. Ch. 9 in Mardia et al.,, 1979). Second,
in the context of questionnaire responses, recent work by Bond & Lang, (2019)
has pointed out that without assumptions on the latent trait that is being
measured by the questionnaire, possible location differences between two
groups of respondents are not identified. Yet, Bond & Lang, (2019) consider
the necessary conditions for parametric and nonparametric identification to be
unlikely to be satisfied in practice. In other words, a possible model on a
latent trait is likely to be misspecified. However, this is exactly the type
of situation where robust estimation could be useful: Robust estimators, like
the one presented in this paper, implicitly allow for model misspecification
while preserving estimation accuracy. Hence, if one is worried about
violations of modeling assumptions like those needed for identification, it
might be attractive to employ a robust estimator that is designed to work
reasonably well despite violations of said assumptions.

##  3 Methodology

Throughout this section, suppose one observes a sample
$\\{\bm{Z}_{i}\\}_{i=1}^{N}$ of $N$ independent categorical draws from
population density $f_{\varepsilon}$ in (4), that is, the assumed model
$\\{\bm{p}\left(\bm{\theta}\right):\bm{\theta}\in\bm{\Theta}\\}$ is possibly
misspecified. In the following, I briefly review the commonly employed maximum
likelihood estimation and then proceed by describing the proposed robust
estimator.

###  3.1 Maximum likelihood estimation

Denote by

$N_{\bm{z}}=\sum_{i=1}^{N}\mathds{1}\left\\{\bm{Z}_{i}=\bm{z}\right\\}$

the empirical frequency of an outcome $\bm{z}\in\bm{\mathcal{Z}}$. The maximum
likelihood estimator (MLE) of the true parameter $\bm{\theta}_{*}$ can be
expressed as

$\widehat{\bm{\theta}}_{N}^{\mathrm{\
MLE}}=\arg\max_{\bm{\theta}\in\bm{\Theta}}\Bigg{\\{}\sum_{\bm{z}\in\bm{\mathcal{Z}}}N_{\bm{z}}\log\big{(}p_{\bm{z}}\left(\bm{\theta}\right)\big{)}\Bigg{\\}}.$
(5)

Many models for categorical data are commonly estimated via maximum likelihood
(and extensions thereof), including the famous Rasch model in item response
theory (Ch. 4 in Mair,, 2018), the polychoric correlation model in eq. 3 (see
Olsson,, 1979), and even sample correlation like the Pearson correlation
coefficient can be expressed as MLE.

However, while the MLE is consistent and efficient for a correctly specified
model, it is well-known to be inconsistent for $\bm{\theta}_{*}$ when the
model is misspecified and is generally highly susceptible to already small
degrees of misspecification (Huber & Ronchetti,, 2009; Hampel et al.,, 1986;
Huber,, 1964).

###  3.2 Proposed estimator

In the following, I propose an estimator of $\bm{\theta}_{*}$ that is
asymptotically equivalent to the MLE in the absence of misspecification, and
substantially more accurate than the MLE in the presence of misspecification.
The estimator is based on the following observation. The empirical relative
frequency of a fixed outcome $\bm{z}\in\bm{\mathcal{Z}}$, denoted

$\widehat{f}_{N}(\bm{z})=N_{\bm{z}}\big{/}N=\frac{1}{N}\sum_{i=1}^{N}\mathds{1}\left\\{\bm{Z}_{i}=\bm{z}\right\\},$

is a pointwise strongly consistent nonparameteric estimator of its
corresponding population cell probability

$f_{\varepsilon}(\bm{z})=(1-\varepsilon)p_{\bm{z}}\left(\bm{\theta}_{*}\right)+\varepsilon
h(\bm{z})$

as $N\to\infty$ (see e.g., Chapter 19.2 in Van der Vaart,, 1998). Cosequently,
if there is no misspecification ($\varepsilon=0$), then model
$\\{\bm{p}\left(\bm{\theta}\right):\bm{\theta}\in\bm{\Theta}\\}$ is correctly
specified because

$f_{0}(\bm{z})=p_{\bm{z}}\left(\bm{\theta}_{*}\right),$

whereas if $\varepsilon>0$, then misspecification arises from the model’s cell
probabilities not being equal to the population cell probabilities, that is,

$f_{\varepsilon}(\bm{z})\neq p_{\bm{z}}\left(\bm{\theta}_{*}\right).$ (6)

Consequently, if the model is correctly specified (and $\bm{\theta}_{*}$ is
identified), then there exists a parameter $\bm{\theta}\in\bm{\Theta}$ whose
associated probability mass function $p_{\bm{z}}\left(\bm{\theta}\right)$ the
nonparametric estimator $\widehat{f}_{N}(\bm{z})$ will converge to, namely the
true $\bm{\theta}=\bm{\theta}_{*}$. Conversely, if there is misspecification,
then there is no $\bm{\theta}\in\bm{\Theta}$ whose model probability
$p_{\bm{z}}\left(\bm{\theta}\right)$ the estimate $\widehat{f}_{N}(\bm{z})$
will converge to. This fact can be exploited in the estimation of the model’s
parameters by minimizing a certain divergence between the empirical relative
outcome frequencies, $\widehat{f}_{N}(\bm{z})$, and theoretical outcome
probabilities $p_{\bm{z}}\left(\bm{\theta}\right)$ returned by the assumed
model (1) at parameter value $\bm{\theta}\in\bm{\Theta}$ to find the most
accurate fit that the assumed model can offer for the observed data.
Specifically, the estimator minimizes with respect to $\bm{\theta}$ the loss

$L\left(\bm{\theta},\widehat{f}_{N}\right)=\sum_{\bm{z}\in\bm{\mathcal{Z}}}{\rho}\Bigg{(}\frac{\widehat{f}_{N}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}\Bigg{)}p_{\bm{z}}\left(\bm{\theta}\right),$
(7)

where $\rho:[0,\infty)\to\mathbb{R}$ is a predefined function that will be
defined momentarily. The proposed estimator $\widehat{\bm{\theta}}_{N}$ is
given by the value minimizing the objective loss over $\bm{\Theta}$,

$\widehat{\bm{\theta}}_{N}=\arg\min_{\bm{\theta}\in\bm{\Theta}}L\Big{(}\bm{\theta},\widehat{f}_{N}\Big{)}.$
(8)

Estimators that minimize a loss function of the type in (7) are called minimum
disparity estimators (Lindsay,, 1994) because they minimize a certain
disparity between empirical probabilities ($\widehat{f}_{N}(\bm{z})$ here) and
theoretical probabilities ($p_{\bm{z}}\left(\bm{\theta}\right)$ here). A small
disparity indicates that the assumed model is able to fit observed data well.
The disparity is governed by the choice of the function $\rho(\cdot)$. Many
estimators can be written as minimum disparity estimators, including the MLE
and various estimators for grouped data (Victoria-Feser & Ronchetti,, 1997;
Cressie & Read,, 1984) through the choice of $\rho(\cdot)$ (Victoria-Feser &
Ronchetti,, 1997; Lindsay,, 1994). In the following, I motivate a specific
choice of $\rho(\cdot)$ that makes the estimator $\widehat{\bm{\theta}}_{N}$
less susceptible to misspecification of the polychoric model.

The fraction $\widehat{f}_{N}(\bm{z})/p_{\bm{z}}\left(\bm{\theta}\right)$ in
(7) is called a Pearson residual333Technically, Lindsay, (1994) defines
Pearson residuals as
$\widehat{f}_{N}(\bm{z})/p_{\bm{z}}\left(\bm{\theta}\right)-1$. I renounce on
subtracting the value 1 because it makes notation simpler in this paper.
(Lindsay,, 1994) and can be interpreted as a goodness-of-fit measure of
outcome $\bm{z}$. Values close to 1 indicate a good fit between data and
assumed model at $\bm{\theta}$, whereas values toward 0 or $+\infty$ indicate
a poor fit. Indeed, outcomes whose observations are primarily generated by a
different distribution than that of the assumed model will generally have a
Pearson residual away from 1 (see eq. 6). Hence, to achieve robustness to
misspecifcation, outcome frequencies that cannot be modeled well by the
assumed model, as indicated by their Pearson residual being away from 1,
should be downweighted in the estimation procedure such that they do not over-
proportionally affect the fit. Such robustness can be achieved by choosing an
appropriate function for $\rho(\cdot)$ in loss function (7). I propose to
choose the following specification, suggested by Ruckstuhl & Welsh, (2001) for
robustly fitting the binomial model,

$\rho(x)=\begin{cases}x\log(x)&\textnormal{ if }x\in[0,c],\\\
x(\log(c)+1)&\textnormal{ if }x>c,\end{cases}$ (9)

where $c\in[1,\infty]$ is a prespecified tuning constant.444Equation (9) is
actually just a special case of a more general formulation in Ruckstuhl &
Welsh, (2001). Figure 1 visualizes this function for the example choice
$c=1.6$. Note that the function $\rho(\cdot)$ is convex, but whether the
estimator’s optimization problem in (8) is also convex depends on the assumed
model.

(a) $\rho(x)$

(b) $\rho^{\prime}(x)$

Figure 1: Visualization of the function $\rho(x)$ in (9) (left panel) and its
derivative (right panel), for $c=1.6$ (vertical dashed blue line).

It is easy to see that for the choice $c=+\infty$ in function $\rho(\cdot)$,
minimizing the loss (7) is equivalent to maximizing the log-likelihood
objective in (5), meaning that the estimator $\widehat{\bm{\theta}}_{N}$ is
equal to $\widehat{\bm{\theta}}_{N}^{\mathrm{\ MLE}}$ for this choice of $c$.
More specifically, if a Pearson residual
$x=\frac{\widehat{f}_{N}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}$ of an
outcome $\bm{z}\in\bm{\mathcal{Z}}$ is such that $x\in[0,c]$ for fixed $c\geq
1$, then the estimation procedure behaves at this cell like in maximum
likelihood estimation. If this Pearson residual $x$ equals 1, then its
associated outcome can be fitted perfectly with the presumed model, so at this
outcome the estimation procedure will behave like maximum likelihood
regardless of the choice of $c\geq 1$. In the absence of misspecification
($\varepsilon=0$),
$\widehat{f}_{N}(\bm{z})\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}p_{\bm{z}}\left(\bm{\theta}_{*}\right)$
as $N\to\infty$ for all $\bm{z}\in\bm{\mathcal{Z}}$, meaning that all Pearson
residuals are asymptotically equal to 1. In other words, if there is no
misspecification, $\widehat{\bm{\theta}}_{N}$ is asymptotically equivalent to
$\widehat{\bm{\theta}}_{N}^{\mathrm{\ MLE}}$ no matter the choice of tuning
constant $c\geq 1$. On the other hand, if an outcome’s Pearson residual is far
away from 1, it cannot be fitted well with the presumed model, which is
typically indicative of the polychoric model being misspecified. In this case,
this particular outcome should not be treated like in maximum likelihood
estimation because maximum likelihood is not consistent under
misspecification. Instead, the outcome’s influence on the final estimate
should be downweighted to avoid that outcomes that cannot be fitted well
dominate the fit, which happens in maximum likelihood. Such downweighting is
employed by function $\rho(x)$ in (9) whenever $x>c\geq 1$, that is, the
Pearson residual is sufficiently far away from the ideal value 1, where the
choice of $c$ governs what is deemed “sufficiently far”. Indeed, for values
$x>c$, the function $\rho(x)$ increases only linearly with $x$, as opposed to
the non-linear exponential increase when $x\leq c$. The notion of requiring
nonlinear effects for the bulk of the data and linear effects in its tails is
similar to classic robust estimation as in Huber, (1964).

It is shown in Figure 1 how $\rho(x)$ transitions from exponential growth to
linear growth at $x=c$, as well as the boundedness of its first derivative
when $c$ is finite. Hence, if $c$ is finite, any Pearson residual can only
have a bounded effect on the final estimator, as opposed to unbounded effects
in maximum likelihood estimation where $c=+\infty$. Thus, I achieve robustness
against misspecification through the choice of $c$. The closer to 1 one
chooses $c$, the more robust the estimator becomes. However, there is a well-
known tradeoff between robustness and efficiency for robust estimators: the
more robust an estimator, the more estimation variance is introduced (e.g.
Huber & Ronchetti,, 2009). Therefore, by choosing $c$, one is effectively
choosing between robustness and efficiency. A characterization of the this
tradeoff is work in progress.

With the proposed choice of $\rho(\cdot)$, I stress that estimator
$\widehat{\bm{\theta}}_{N}$ in (8) has the same time complexity as maximum
likelihood, that is, $O\big{(}\\#\bm{\mathcal{Z}}\big{)}$, since one needs to
calculate the Pearson residual of all possible outcomes for any given
candidate parameter value. Consequently, the proposed estimator does not incur
any additional computational cost compared to maximum likelihood, and
therefore robustness can be achieved without having to pay a computational
price.

###  3.3 Estimand

Considering the potential presence of model misspecification, it is worth
studying what quantity the proposed estimator $\widehat{\bm{\theta}}_{N}$ in
(8) estimates. In population, its estimand is given by

$\bm{\theta}_{0}=\arg\min_{\bm{\theta}\in\bm{\Theta}}L\Big{(}\bm{\theta},f_{\varepsilon}\Big{)}.$
(10)

This minimization problem is simply the population analogue to the
minimization problem in (7) that the sample-based estimator
$\widehat{\bm{\theta}}_{N}$ solves, since the $\widehat{f}_{N}(\bm{z})$ are
the population analogues to the $f_{\varepsilon}(\bm{z})$. In the absence of
misspecification ($\varepsilon=0$), estimand $\bm{\theta}_{0}$ equals the true
parameter $\bm{\theta}_{*}$. In the presence of misspecification
($\varepsilon>0$) it is generally different from $\bm{\theta}_{*}$. How much
different it is depends on the degree and type of of misspecification as well
as the choice of tuning constant $c$ in $\rho(\cdot)$. In general, the larger
$\varepsilon$ (more severe misspecification) and $c$ (less downweighting of
hard-to-fit cells), the further $\bm{\theta}_{0}$ is away from
$\bm{\theta}_{*}$. Hence, for fixed misspecification degree $\varepsilon$, the
MLE ($c=+\infty$) will estimate a parameter that is farther or equally far
away from the true $\bm{\theta}_{*}$ than for finite choices of $c$.
Correspondingly, finite choices of $c$ lead to an estimator that is at least
as accurate as the MLE, and more accurate under misspecification of the
presumed model.

###  3.4 Assumptions

In the following, I list a set of assumptions that will be entertained in the
asymptotic analysis of the proposed estimator $\widehat{\bm{\theta}}_{N}$
computed on random sample $\\{\bm{Z}_{i}\\}_{i=1}^{N}$.

######  Assumption Set A.

Suppose that the following assumptions hold true.

  1. A.1

$c\in[1,+\infty]$

  2. A.2

$\bm{\Theta}\subset\mathbb{R}^{d}$ is compact,

  3. A.3

$\bm{\theta}_{0}=\arg\min_{\bm{\theta}\in\bm{\Theta}}L(\bm{\theta},f_{\varepsilon})$
is a unique global minimum, and $\bm{\theta}_{0}$ is an interior point of
$\bm{\Theta}$,

  4. A.4

$p_{\bm{z}}\left(\bm{\theta}\right)$ is continuously differentiable with
respect to $\bm{\theta}\in\bm{\Theta}$ and twice differentiable at
$\bm{\theta}_{0}$ for all cells $\bm{z}\in\bm{\mathcal{Z}}$,

  5. A.5

$\left\|\frac{\partial
p_{\bm{z}}\left(\bm{\theta}\right)}{\partial\bm{\theta}}\right\|<\infty$ for
all $\bm{\theta}\in\bm{\Theta},\bm{z}\in\bm{\mathcal{Z}}$,

  6. A.6

$p_{\bm{z}}\left(\bm{\theta}\right)>0$ for all $\bm{\theta}\in\bm{\Theta}$,
$\bm{z}\in\bm{\mathcal{Z}}$.

  7. A.7

$\\#\big{\\{}\bm{z}\in\bm{\mathcal{Z}}:f_{\varepsilon}(\bm{z})>0\big{\\}}>d$,

  8. A.8

$L(\bm{\theta},f_{\varepsilon})$ is convex in a neighborhood of
$\bm{\theta}_{0}$,

  9. A.9

$\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\neq c$
for any $\bm{z}\in\bm{\mathcal{Z}}$.

Assumption A.1 ensures that function $\rho(\cdot)$ exhibits meaningful
behavior when evaluated at Pearson residuals, such as the ideal residual value
1 being included in the interval $[0,c]$. Compactness of the parameter space
(Assumption A.2) is primitive and required for a technicality when proving
consistency of $\widehat{\bm{\theta}}_{N}$.555This assumption can possibly be
modified to $\bm{\Theta}$ being open by equipping it with a specific
topological structure. Uniqueness of a global minimum in the parameter space’s
interior (Assumption A.3) is a common assumption in the literature on $M$-type
estimators (e.g., Chapter 5.2 in Van der Vaart,, 1998) like the one presented
in this paper. Assumptions A.4 and A.5 pertain to the presumed model and
require smoothness and first-order boundedness, as well as second-order
differentiability at estimand $\bm{\theta}_{0}$. These assumptions are
standard in the literature on minimum-disparity-type estimators (e.g., Cressie
& Read,, 1984; Ruckstuhl & Welsh,, 2001; Victoria-Feser & Ronchetti,, 1997).
The assumption of strictly positive model probabilities (A.6) is also standard
in this literature, and rules out that one divides by zero when calculating
Pearson residuals. Related is Assumption A.7, which imposes that the number of
positive cell population probabilities is strictly larger than the number of
model parameters. In other words, there must be more sources of variation
(populated cells) than parameters. One may view this assumption as a rank
condition that ensures invertibility of the Hessian matrix of the minimization
problem (7), and is required to prevent rank deficiency of the asymptotic
covariance matrix of estimator $\widehat{\bm{\theta}}_{N}$. Since
$\widehat{f}_{N}(\bm{z})>0$ implies $f_{\varepsilon}(\bm{z})>0$, this
assumption can be empirically verified. Assumption A.8 is a local convexity
assumption that only becomes relevant when the population minimization problem
(10) is not convex, which can happen for some models. Specifically, this
assumption rules out that the gradient of population loss
$L(\bm{\theta},f_{\varepsilon})$ is flat at value zero in a region around
$\bm{\theta}_{0}$, because such a situation would precipitate an
identification problem: If the gradient is zero-valued in a neighborhood of
estimand $\bm{\theta}_{0}$ instead of being zero at $\bm{\theta}_{0}$ only,
then the estimand $\bm{\theta}_{0}$ is not uniquely determined. As such, this
assumption refines the assumption of $\bm{\theta}_{0}$ being a global minimum
of the population loss (Assumption A.3) by requiring well-separatedness of the
global minimum, which is a common assumption in the asymptotic analysis of
$M$-type estimators (see Ch. 5.2 in Van der Vaart,, 1998, for a discussion).
Finally, Assumption A.9 imposes that the Pearson residual at the global
minimum is not equal to tuning constant $c$. This is a primitive condition
that is required for the loss to be twice differentiable at $\bm{\theta}_{0}$
(in combination with Assumption A.4), which is a requirement for the existence
of the estimator’s asymptotic covariance matrix.

I emphasize that no assumption in Assumption Set A restricts the type, source,
or magnitude of potential misspecification of the considered model. In fact,
Assumptions A.2–A.8 are also required for consistency and asymptotic normality
of the MLE. Only assumptions A.1 and A.9 are specific to the proposed
estimator because they pertain to tuning constant $c$.

###  3.5 Asymptotic analysis

The following theorem establishes strong consistency of
$\widehat{\bm{\theta}}_{N}$ for $\bm{\theta}_{0}$. The theorem’s proof and
that of all subsequent mathematical statements are given in Appendix D.

######  Theorem 1 (Consistency).

Under Assumptions A.1–A.6, it holds true that

$\widehat{\bm{\theta}}_{N}\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}\bm{\theta}_{0},$

as $N\to\infty$.

I now study the limit distribution of the estimator. Doing so necessitates
additional notation. For fixed tuning constant $c\geq 1$, let

$w(x)=\mathds{1}\left\\{x\in[0,c]\right\\}+c\mathds{1}\left\\{x>c\right\\}/x\quad\mathrm{for}\
x\geq 0,$

with first derivative

$w^{\prime}(x)=0\mathds{1}\left\\{x\in[0,c]\right\\}-c\mathds{1}\left\\{x>c\right\\}/x^{2},$

and further define the $d$-dimensional gradient of
$\log\left(p_{\bm{z}}\left(\bm{\theta}\right)\right)$ for cell
$\bm{z}\in\bm{\mathcal{Z}}$ at parameter $\bm{\theta}\in\bm{\Theta}$ as

$\bm{s}_{\bm{z}}\left(\bm{\theta}\right)=\frac{1}{p_{\bm{z}}\left(\bm{\theta}\right)}\left(\frac{\partial}{\partial\bm{\theta}}\
p_{\bm{z}}\left(\bm{\theta}\right)\right),$

as well as the $d\times d$ Hessian matrix of
$\log(p_{\bm{z}}\left(\bm{\theta}\right))$ as

$\bm{Q}_{\bm{z}}(\bm{\theta})=\frac{1}{p_{\bm{z}}\left(\bm{\theta}\right)}\left(\frac{\partial^{2}}{\partial{\bm{\theta}}\partial{\bm{\theta}}^{\top}}p_{\bm{z}}\left(\bm{\theta}\right)\right)-\bm{s}_{\bm{z}}\left(\bm{\theta}\right)\bm{s}_{\bm{z}}\left(\bm{\theta}\right)^{\top}.$

In addition, define $d$-dimensional vectors

$\bm{w}_{\bm{z}}\left(\bm{\theta}\right)=\bm{s}_{\bm{z}}\left(\bm{\theta}\right)\mathds{1}\left\\{\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}\in[0,c]\right\\},$

and a $d\times m$ matrix

$\displaystyle\bm{W}(\bm{\theta})=\bigg{(}\bm{w}_{\bm{z}_{1}}\left(\bm{\theta}\right),\bm{w}_{\bm{z}_{2}}\left(\bm{\theta}\right),\cdots,\bm{w}_{\bm{z}_{m}}\left(\bm{\theta}\right)\bigg{)}$

that row-binds all $m$ vectors $\bm{w}_{\bm{z}}\left(\bm{\theta}\right)$ for
fixed $\bm{\theta}$. The $j$-th column of this matrix equals the gradient
$\bm{s}_{\bm{z}_{j}}(\bm{\theta})$ if the Pearson residual of the associated
cell $\bm{z}_{j},j=1,\dots,m,$ is in the MLE-part of the function
$\rho(\cdot)$, and is otherwise equal to a $d$-dimensional zero. In similar
fashion, define the $m$-dimensional vector

$\bm{f}_{\varepsilon}=\Big{(}f_{\varepsilon}(\bm{z}_{1}),f_{\varepsilon}(\bm{z}_{2}),\dots,f_{\varepsilon}(\bm{z}_{m})\Big{)}^{\top}$

that holds all $m$ evaluations of the population density $f_{\varepsilon}$,
and put

$\bm{\Omega}=\mathrm{diag}(\bm{f}_{\varepsilon})-\bm{f}_{\varepsilon}\bm{f}_{\varepsilon}^{\top}.$

With this notation, I can establish root-$N$ consistency and asymptotic
normality of estimator $\widehat{\bm{\theta}}_{N}$ in the following theorem.

######  Theorem 2 (Asymptotic normality).

Grant the assumptions of Assumption Set A. Then

$\sqrt{N}\left(\widehat{\bm{\theta}}_{N}-\bm{\theta}_{0}\right)\stackrel{{\scriptstyle\mathrm{d}}}{{\longrightarrow}}\textnormal{N}_{d}\Big{(}\bm{0},\bm{\Sigma}\left({\bm{\theta}_{0}}\right)\Big{)},$

as $N\to\infty$, where, as a function of $\bm{\theta}\in\bm{\Theta}$, the
invertible asymptotic covariance matrix is given by

$\bm{\Sigma}\left({\bm{\theta}}\right)=\bm{M}\left(\bm{\theta}\right)^{-1}\bm{U}\left(\bm{\theta}\right)\bm{M}\left(\bm{\theta}\right)^{-1},$

where

$\displaystyle\bm{U}\left(\bm{\theta}\right)$
$\displaystyle=\bm{W}\left(\bm{\theta}\right)\bm{\Omega}\bm{W}\left(\bm{\theta}\right)^{\top}\qquad\text{and}$
$\displaystyle\bm{M}\left(\bm{\theta}\right)$
$\displaystyle=\sum_{\bm{z}\in\bm{\mathcal{Z}}}f_{\varepsilon}(\bm{z})\Bigg{(}{w^{\prime}}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}\Bigg{)}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}\bm{s}_{\bm{z}}\left(\bm{\theta}\right)\bm{s}_{\bm{z}}\left(\bm{\theta}\right)^{\top}-{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}\Bigg{)}\bm{Q}_{\bm{z}}(\bm{\theta})\Bigg{)}$

are $d\times d$ symmetric and invertible matrices.

A strongly consistent estimator of the unobserved asymptotic covariance matrix
$\bm{\Sigma}\left({\bm{\theta}_{0}}\right)$ can be constructed as follows.
Replace all population class probabilities $f_{\varepsilon}(\bm{z})$ by their
corresponding empirical counterparts $\widehat{f}_{N}(\bm{z})$ in matrices
$\bm{W}\left(\bm{\theta}\right),\bm{M}\left(\bm{\theta}\right)$, and
$\bm{\Omega}$. Then exploit the plug-in principle and evaluate
$\bm{U}\left(\bm{\theta}\right)$ and $\bm{M}\left(\bm{\theta}\right)$ at the
point estimate $\widehat{\bm{\theta}}_{N}.$ Denote the ensuing plug-in
estimator by $\bm{\Sigma}\left({\widehat{\bm{\theta}}_{N}}\right)$, which is
strongly consistent for $\bm{\Sigma}\left({\bm{\theta}_{0}}\right)$ by Theorem
1 and the continuous mapping theorem.

Moreover, it can be shown that in the absence of misspecification, the
asymptotic covariance matrix of the proposed estimator,
$\bm{\Sigma}(\bm{\theta}_{0})$ from Theorem 2, is equal to the asymptotic
covariance matrix of the MLE whenever $c>1$ (see Lemma B.5 in the Appendix).
Hence, if the model is correctly specified, both estimators are first and
second order equivalent.

###  3.6 Goodness-of-fit test

Suppose one wishes to test the null hypothesis that an individual cell in a
$k$-way contingency table can be fitted well by the presumed model. Rejecting
this null hypothesis is indicative of the model’s misspecification, at least
for that cell. This notion can be conceptualized by means of Pearson
residuals. Recall that a Pearson residual of value 1 indicates that the
corresponding cell can be fitted well, whereas a Pearson residual
significantly larger than 1 indicates poor fit. For a given cell
$\bm{z}\in\bm{\mathcal{Z}}$, this translates into the natural null hypothesis
with one-sided alternative

$\textnormal{H}_{0}:\frac{\widehat{f}_{N}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}=1\qquad\text{vs.}\qquad\textnormal{H}_{1}:\frac{\widehat{f}_{N}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}>1,$

which is equivalent to

$\textnormal{H}_{0}:p_{\bm{z}}\left(\bm{\theta}_{0}\right)=\widehat{f}_{N}(\bm{z})\qquad\text{vs.}\qquad\textnormal{H}_{1}:p_{\bm{z}}\left(\bm{\theta}_{0}\right)<\widehat{f}_{N}(\bm{z}).$
(11)

Ideally, a test for such a hypothesis will reject $\textnormal{H}_{0}$ if the
presumed model is misspecified for that cell, and sustain $\textnormal{H}_{0}$
if it is correctly specified for that cell. It turns out that a test statistic
that satisfies these two desirable properties is given by

$T_{N}(\bm{z})=\frac{p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)-\widehat{f}_{N}(\bm{z})}{\sqrt{\sigma_{\bm{z}}^{2}\left(\bm{\theta}_{0}\right)\big{/}N}},$
(12)

where

$\sigma_{\bm{z}}^{2}\left(\bm{\theta}\right)=\bm{g}_{\bm{z}}(\bm{\theta})^{\top}\bm{\Sigma}\left({\bm{\theta}}\right)\bm{g}_{\bm{z}}(\bm{\theta})$

for gradient

$\bm{g}_{\bm{z}}(\bm{\theta})=\frac{\partial
p_{\bm{z}}\left(\bm{\theta}\right)}{\partial\bm{\theta}}.$

The intuition behind using test statistic $T_{N}(\bm{z})$ to test
$\textnormal{H}_{0}:p_{\bm{z}}\left(\bm{\theta}_{0}\right)=\widehat{f}_{N}(\bm{z})$
is as follows. If the presumed model is correctly specified ($\varepsilon=0$
so that $\bm{\theta}_{0}=\bm{\theta}_{*}$ ), it holds true that
$\widehat{f}_{N}(\bm{z})\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}p_{\bm{z}}\left(\bm{\theta}_{0}\right)$
as well as
$p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}p_{\bm{z}}\left(\bm{\theta}_{0}\right)$
by Theorem 1 and the continuous mapping theorem. It follows that the
difference between $p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)$ and
$\widehat{f}_{N}(\bm{z})$ vanishes as $N$ grows to infinity. Indeed, the
following corollary (Corollary 3.1) shows that this difference, when scaled
appropriately to equal $T_{N}(\bm{z})$, converges in distribution to a zero-
mean Gaussian distribution, provided that the model is correctly specified.
However, if the model is misspecified ($\varepsilon>0$), then
$\widehat{f}_{N}(\bm{z})\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}f_{\varepsilon}(\bm{z})$
and
$p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)\stackrel{{\scriptstyle\mathrm{a.s.}}}{{\longrightarrow}}p_{\bm{z}}\left(\bm{\theta}_{0}\right)$,
but $f_{\varepsilon}(\bm{z})\neq p_{\bm{z}}\left(\bm{\theta}_{0}\right)$
because of the misspecification. Consequently, the difference between
$p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)$ and
$\widehat{f}_{N}(\bm{z})$ does not converge to zero (and need not converge at
all). Hence, if $T_{N}(\bm{z})$ is statistically significantly different from
a zero-mean Gaussian sequence, there is evidence that the presumed model is
misspecified, at least at cell $\bm{z}$. The following corollary of Theorem 2
formally establishes the validity of test statistic $T_{N}(\bm{z})$ for
testing the null hypothesis
$\textnormal{H}_{0}:\widehat{f}_{N}(\bm{z})=p_{\bm{z}}\left(\bm{\theta}_{0}\right)$.

######  Corollary 3.1 (Limit distribution of test statistic).

Grant the assumptions of Assumption Set A and consider a given cell
$\bm{z}\in\bm{\mathcal{Z}}$. Then, under the null hypothesis in (11), the test
statistic $T_{N}(\bm{z})$ in (12) possesses the following limit distribution:

$T_{N}(\bm{z})\stackrel{{\scriptstyle\mathrm{d}}}{{\longrightarrow}}\textnormal{N}(0,1),$

as $N\to\infty$.

In practice, the variance term
$\sigma_{\bm{z}}^{2}\left(\bm{\theta}_{0}\right)$ in test statistic
$T_{N}(\bm{z})$ is unobserved, but it can be strongly consistently estimated
by $\sigma_{\bm{z}}^{2}\left(\widehat{\bm{\theta}}_{N}\right)$, which follows
from Theorem 1 and the continuous mapping theorem. Hence, in practice, one
uses the approximate test statistic

$\widehat{T}_{N}(\bm{z})=\frac{p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)-\widehat{f}_{N}(\bm{z})}{\sqrt{\sigma_{\bm{z}}^{2}\left(\widehat{\bm{\theta}}_{N}\right)\big{/}N}}$

for hypothesis testing. In addition, it should be noted that using this test
to test the null hypothesis in (11) for all $m$ outcomes in sample space
$\bm{\mathcal{Z}}$ creates a multiple comparisons problem. I therefore
recommend to adjust the ensuing $p$-values for multiple comparisons, for
instance through the procedure of Benjamini & Hochberg, (1995).

As an additional theoretical result, I derive the influence function of
estimator $\widehat{\bm{\theta}}_{N}$ in Appendix B. It turns out that for
tuning constant choices $c>1$, the influence function equals that of the non-
robust MLE, but remains bounded at categorical data. Similar results have been
derived for minimum power divergence estimators (Victoria-Feser & Ronchetti,,
1997), minimum Hellinger distance estimators (Simpson,, 1987; He & Simpson,,
1993; Lindsay,, 1994), and minimum disparity estimators at the binomial model
(Ruckstuhl & Welsh,, 2001). However, I argue that the influence function is
not a very informative measure when working with categorical data because the
influence of a single categorical observation is bounded by construction due
to the bounded nature of categorical data. A detailed discussion is provided
in Appendix B.

###  3.7 Computational and practical aspects

The minimization problem of the estimator in (7) can be solved with standard
gradient descent methods, which allows one to impose constraints on the model
parameters. In the spirit of open and reproducible science and to enhance
accessibility and adoption by empirical researchers, an R implementation of
the proposed methodology is freely available in the package robcat (Welz,,
2024, for “ROBust CATegorical data analysis”). To maximize speed and
performance, the package is predominantly developed in C++ and integrated to R
via Rcpp (Eddelbuettel,, 2013).

As for the tuning constant $c\in[1,+\infty]$ in function $\rho(\cdot)$ in (9),
the choice $c=1.6$ yielded a good compromise between robustness and efficiency
in numerous simulation experiments. I therefore recommend this choice for
practical use, but a detailed investigation on the optimal choice of $c$ is
work in progress.

##  4 Monte Carlo experiments

In order to verify the theoretical guarantees of the proposed estimator and
demonstrate its performance in practice, I employ a series of Monte Carlo
simulation experiments. In these experiments, I emulate a SEM on questionnaire
responses where some respondents only give inattentive responses. As described
in Section 2.2, polychoric correlation coefficients estimated with my proposed
robust estimator is the workhorse behind robustifying SEMs against inattentive
responding (among other sources of misspecification). Therefore, the
simulation design is twofold: In the first design, I study estimation of the
pairwise polychoric model in (3), while the second design focuses on fitting
SEMs.

###  4.1 Polychoric correlation

Consider a pairwise polychoric correlation model for two rating items (eq. 3).
Let there be $J_{x}=J_{y}=5$ response categories for each of the two rating
variables $\bm{Z}=(X,Y)^{\top}$ and define the true thresholds in the
polychoric model as

$a_{*,1}=b_{*,1}=-1.5,\quad a_{*,2}=b_{*,2}=-0.5,\quad
a_{*,3}=b_{*,3}=0.5,\quad a_{*,4}=b_{*,4}=1.5,$ (13)

and let the true polychoric correlation coefficient be $\varrho_{*}=0.5$.
Then, for a bivariate standard normally distributed latent random variable
with correlation coefficient $\varrho_{*}$, denoted $(\xi,\eta)$, observations
of rating responses $(X,Y)$ are generated by

$X=\begin{cases}1&\text{if }\xi<a_{*,1},\\\ 2&\text{if
}a_{*,1}\leq\xi<a_{*,2},\\\ 3&\text{if }a_{*,2}\leq\xi<a_{*,3},\\\ 4&\text{if
}a_{*,3}\leq\xi<a_{*,4},\\\ 5&\text{if }a_{*,4}\leq\xi,\\\ \end{cases}$ and
$Y=\begin{cases}1&\text{if }\eta<b_{1},\\\ 2&\text{if
}b_{*,1}\leq\eta<b_{*,2},\\\ 3&\text{if }b_{*,2}\leq\eta<b_{*,3},\\\
4&\text{if }b_{*,3}\leq\eta<b_{*,4},\\\ 5&\text{if }b_{*,4}\leq\eta.\\\
\end{cases}$ (14)

To emulate misspecification of the polychoric model, let a fraction
$\varepsilon$ of the latent $(\xi,\eta)$ be generated by a bivariate normal
distribution with mean $(2,-2)^{\top}$, variances $(0.2,0.2)^{\top}$, and zero
covariance (and therefore zero correlation) that causes the polychoric model
(3) to be misspecified. Hence, in simulated data with nonzero misspecification
$\varepsilon$, the empirical frequency of cells
$(x,y)\in\\{(5,1),(4,3),(5,2)\\}$ will be inflated in the sense that they have
a higher realization probability than under the true polychoric model
distribution. The data points causing these three cells to be inflated are
instances of negative leverage points. Here, such leverage points drag
correlational estimates away from a positive value towards zero or, if there
are sufficiently many of them, even a negative value.

For misspecification degrees $\varepsilon\in\\{0,0.1,0.2\\}$, I sample
$N=1,000$ rating responses from this data generating process and estimate the
true parameter
$\bm{\theta}_{*}=\left(\varrho_{*},a_{*,1},\dots,a_{*,4},b_{*,1},\dots,b_{*,4}\right)^{\top}$
by means of the Pearson sample correlation coefficient (for $\varrho_{*}$),
the MLE, as well as the proposed estimator with tuning constant set to
$c=1.6$, since this choice yielded a good compromise between robustness and
efficiency in further simulation studies. This procedure is repeated for 1,000
simulated datasets. As performance measures, I calculate the average bias,
standard deviation across repetitions, coverage, and length of confidence
intervals at significance level $\alpha=0.05$. The coverage is defined as
proportion of $(1-\alpha)$-th confidence intervals $[\widehat{\varrho}_{N}\mp
q_{1-\alpha/2}\cdot\text{SE}(\widehat{\varrho}_{N})]$ that contain the true
$\varrho_{*}$, where $q_{1-\alpha/2}$ is the $(1-\alpha/2)$-th quantile of the
standard normal distribution and $\text{SE}(\widehat{\varrho}_{N})$ is the
standard error of $\widehat{\varrho}_{N}$, which is constructed using the
limit theory developed in Theorem 2. The length of a confidence interval is
given by $2\cdot q_{1-\alpha/2}\cdot\text{SE}(\widehat{\varrho}_{N})$.

Figure 2: Boxplot visualization of the bias ($\widehat{\varrho}_{N}-\varrho_{*}$) of the three estimators, Pearson sample correlation, the MLE, and the proposed robust estimator with $c=1.6$, for various degrees of misspecification across 1,000 simulated datasets. Misspecification | Estimator | $\widehat{\varrho}_{N}$ | Bias | StDev | Coverage | CI length
---|---|---|---|---|---|---
$\varepsilon=0$ | Polycor robust | 0.504 | 0.004 | 0.027 | 0.930 | 0.104
Polycor ML | 0.500 | 0.000 | 0.026 | 0.943 | 0.102
Sample cor | 0.457 | $-0.043$ | 0.025 | 0.702 | 0.110
$\varepsilon=0.1$ | Polycor robust | 0.466 | $-0.034$ | 0.038 | 0.911 | 0.152
Polycor ML | 0.097 | $-0.403$ | 0.029 | 0.000 | 0.134
Sample cor | 0.084 | $-0.416$ | 0.026 | 0.000 | 0.124
$\varepsilon=0.2$ | Polycor robust | 0.439 | $-0.061$ | 0.051 | 0.951 | 0.220
Polycor ML | $-0.172$ | $-0.672$ | 0.028 | 0.000 | 0.133
Sample cor | $-0.133$ | $-0.633$ | 0.026 | 0.000 | 0.123
Table 1: Performance measures of of the three estimators, Pearson sample
correlation, the MLE, and the proposed robust estimator with $c=1.6$, for
various degrees of misspecification across 1,000 simulated datasets. The true
polychoric correlation coefficient is $\rho_{*}=0.5$. The performance measures
are the average point estimate of the polychoric correlation coefficient,
$\widehat{\varrho}_{N}$, average bias ($\widehat{\varrho}_{N}-\varrho_{*}$),
the standard deviation of the $\widehat{\varrho}_{N}$ (“StDev”), the
estimator’s coverage with respect to the true $\varrho_{*}$ at significance
level $\alpha=0.05$, as well as the length of the estimator’s confidence
interval, again at level $\alpha=0.05$.

Figure 2 visualizes the bias of each estimator with respect to the true
polychoric correlation $\varrho_{*}$ across the 1,000 simulated datasets.
Analogous plots for the whole parameter $\bm{\theta}_{*}$ can be found in
Appendix C; the results are similar to those of $\varrho_{*}$. For all
considered misspecification degrees, the estimates of the MLE and sample
correlation are similar, which is expected because these two estimators are
known to yield similar similar results when there are five or more rating
options (cf. Rhemtulla et al.,, 2012). In the absence of misspecification,
both MLE and the robust estimator yield accurate estimates. Both estimates are
nearly equivalent to one another in the sense that their point estimates,
standard deviation, and coverage at significance level $\alpha=0.05$ are very
similar (Table 1). However, when misspecification is introduced, MLE, sample
correlation, and robust estimator yield noticeably different results. At
misspecification degree $\varepsilon=0.1$, MLE and sample correlation are
substantially biased with average estimates of $0.097$ and $0.084$,
corresponding to biases of $-0.403$ and $-0.416$, respectively, as well as
zero coverage. In contrast, the robust estimator maintains accuracy with an
average estimate of $0.466$, which corresponds to only a minor bias of
$-0.034$ and a good coverage of $0.911$ (Table 1). When the misspecification
is increased to $\varepsilon=0.2$, the contrast between the estimators becomes
even stronger. While the robust estimator is still remarkably close to the
truth with a small bias of $-0.061$, MLE and sample correlation produce
estimates that are not only severely biased (biases of $-0.672$ and $-0.633$),
but also sign-flipped: While the true correlation is strongly positive (0.5),
both estimates are considerably negative ($-0.172$ and $-0.133$). It is worth
noting that in the presence of misspecification, the confidence intervals of
the robust estimator are wider than those of the MLE (see Table 1). This is
expected because of the well-known trade-off between robustness and
efficiency: An estimator that is designed to reduce bias, like a robust
estimator, will inevitably have a larger estimation variance (e.g. Huber &
Ronchetti,, 2009). These wider confidence intervals furthermore explain why
the robust estimator improves its coverage in Table 1 when the degree of
misspecification is increased from 0.1 to 0.2.

This first simulation study demonstrated that already a small degree of
misspecificaion of the polychoric model can render the commonly employed MLE
and sample correlation unreliable, while the proposed robust estimator retains
good accuracy even in the presence of considerable misspecificaion. On the
other hand, when the model is correctly specified, MLE and robust estimator
produce equivalent results. With these affirmative results in mind, I proceed
to the second simulation study, where a structural equation model is fitted to
a robustly estimated polychoric correlation matrix that was computed on data
for which the model is misspecified due to inattentive responding.

###  4.2 Structural Equation Modeling

This second simulation design emulates a popular type of SEM, namely a
confirmatory factor analysis (CFA). In a CFA, a scale compprising multiple
rating items is tested for validity and reliability. That is, does the scale
indeed measure the latent construct it is supposed to measure, and does it so
reliably?

Suppose there are $q=4$ rating items that jointly measure the same
unidimensional ($r=1$) latent construct in a reliable manner by having a
$q\times r$ factor loadings matrix

$\bm{\Lambda}=(0.75,0.75,0.75,0.75)^{\top}$

and diagonal $q\times q$ noise covariance matrix $\bm{\Psi}$ with diagonal
elements all being equal to 0.4375. Then, by factor model (2), the ensuing
$q\times q$ population correlation matrix $\bm{\Sigma}$ is characterized by
having the same pairwise correlation of $0.5625$ between all distinct items.
It follows that this scale is indeed reliable for the latent construct with a
population Cronbach-$\alpha$ of 0.84.666Cronbach’s $\alpha$ (Cronbach,, 1951)
is a lower bound on the reliability of a scale, and is computed from the
scale’s correlation matrix. In general, a scale with Cronbach-$\alpha$ of 0.75
or higher is considered reliable (Robinson,, 2018).

I then generate $N=1,000$ responses to the $q=4$ rating items with $J=5$
response options each by drawing from a zero-mean normal distribution with
correlation matrix $\bm{\Sigma}$ and use these draws to generate rating
responses according to discretization process (14), where the same thresholds
as in (13) are used for each variable. To emulate inattentive responding, I
replace a fraction $\varepsilon\in\\{0,0.1,0.2\\}$ of the generated rating
observations with the negative leverage point $(1,5,1,5)^{\top}$. Such a
response vector corresponds to a respondent who alternates between the most
extreme response points, regardless of item content.777This response vector
may alternatively be interpreted as the responses of a straightliner at the
first option, after recoding negatively worded items (which would here be the
second and fourth item). Based on the resulting $N\times q$ data matrix of
ratings, I estimate the correlation matrix $\bm{\Sigma}$ by means of the
standard Pearson sample correlation matrix as well as pairwise polychoric
correlation coefficients estimated robustly (with choice $c=1.6$) and via MLE.
Then, based on the estimated correlation matrix $\bm{\Sigma}$, the factor
model (2) is fitted with $r=1$ factor by means of maximum
likelihood888Alternatively, one could also fit a factor model via principal
components or least-squares approaches. to obtain estimates of loadings matrix
$\bm{\Lambda}$ and error covariance matrix $\bm{\Psi}$. In addition, I
estimate the scale’s reliability by calculating an estimate of
Cronbach’s-$\alpha$ based on the estimated correlation matrix. This procedure
is repeated 1,000 times. As performance measures, I calculate the root mean
squared error (RMSE) between estimate and true value of each element in
matrices $\bm{\Lambda}$ and $\bm{\Sigma}$ as well as Cronbach’s $\alpha$.

Figure 3: Boxplot visualization of the RMSE of loadings matrix $\bm{\Lambda}$,
correlation matrix $\bm{\Sigma}$, and Cronbach’s $\alpha$ for the three
estimators for various degrees of misspecification across 1,000 simulated
datasets.

Figure 3 visualizes the simulation results by means of boxplots. Like in the
simulation design for pairwise polychoric correlations (Figure 2), MLE-based
estimation of the polychoric correlation matrix yields similar results as
estimation based on the sample correlation matrix, which is in line with
expectations (cf. Rhemtulla et al.,, 2012). In the absence of inattentive
responding ($\varepsilon=0$), both polychoric estimates (MLE and robust) are
equivalent and accurate with respect to all three performance measures.
However, once inattentive responding is introduced, the MLE and sample
correlation-based estimates display a large bias: At misspecification degree
$\varepsilon=0.1$, both have an RMSE of about 0.2 and 0.35 for the loadings
matrix and correlation matrix, respectively, as well as an inaccurate estimate
of Cronbach’s $\alpha$: While the true value of Cronbach’s $\alpha$ of 0.84
indicates a reliable scale, the two non-robust estimators would suggest an
rather unreliable scale with an estimated reliability of only about 0.65. In
contrast, the robust estimator remains accurate with respect to all three
performance measures. In addition, while the performance of the two non-robust
estimators further deteriorates with increasing prevalence of inattentive
responding $(\varepsilon=0.2)$, that of the robust estimator remains stable
and almost unaffected.

To conclude this section, my Monte Carlo experiments demonstrated the
performance and practical usefulness of the proposed robust estimator when
applied to a SEM: While hitherto estimation approaches are highly sensitive to
the adverse effects of inattentive responding, the robust estimator remains
almost unaffected by their presence although it makes no assumption about the
inattention’s magnitude or type. On the other hand, in the absence of
inattentive responding, the robust estimator yields equivalent results to
hitherto estimators. These simulation results suggest that a SEM can be
robustified against inattentive responding by using the novel robust estimator
proposed in this paper.

##  5 Empirical application

###  5.1 Background and study design

In this section, I demonstrate the proposed methodology on empirical data by
using a subset of the 100 unipolar markers of the Big-5 personality traits
(Goldberg,, 1992).999The Big-5 factor model is a fundamental model in
personality psychology. It assumes that human personality can be described by
five latent variables (traits), namely openness, conscientiousness,
extraversion, agreeableness, and neuroticism. Each marker is a questionnaire
item comprising a single English adjective (such as “bold” or “timid”) asking
respondents to indicate how accurately the adjective describes their
personality using a 5-point Likert-type rating scale (very inaccurate,
moderately inaccurate, neither accurate nor inaccurate, moderately accurate,
and very accurate). Here, each Big-5 personality trait is measured with six
pairs of adjectives that are polar opposites to one another (such as
“talkative” vs. “silent”), that is, twelve items in total for each trait. It
seems implausible that an attentive respondent would choose to agree (or
disagree) to both items in a pair of polar opposite adjectives. Consequently,
one expects a strongly negative correlation between polar adjectives if all
respondents respond attentively (Arias et al.,, 2020). Hence, this dataset is
well-suited to test if my robust estimator can correctly identify and account
for inattentive respondents, which are presumably those who gave inconsistent
responses to polar opposite adjective items.

Arias et al., (2020) collect measurements of three Big-5 traits in this way,
namely extroversion, neuroticism, and conscientiousness.101010Arias et al.,
(2020) synonymously refer to neuroticism as emotional stability. Furthermore,
in addition to the three listed traits, Arias et al., (2020) collect
measurements of the trait dispositional optimism by using a different
instrument, and gather data for another scale that is designed to not measure
any construct. For the sake of brevity, I do not consider these scales in this
empirical demonstration. The sample that I shall use, Sample 1 in Arias et
al., (2020), consists of $N=725$ online respondents who are all U.S. citizens,
native English speakers, and tend to have relatively high levels of reported
education (about 90% report to hold an undergraduate or higher degree).
Concerned about respondent inattention in their data, Arias et al., (2020)
construct a factor mixture model for detecting inattentive participants. Their
model crucially relies on response inconsistencies to polar opposite
adjectives and is designed to primarily detect inattentive straightlining
responding. They find that inattentive responding is a sizable problem in
their data. Their model estimates that the proportion of inattentive
participants amounts to 4.7% in the conscientiousness, 6% in the neuroticism,
and 7.3% in the extroversion scale. After some further analyses, the authors
conclude that if unaccounted for, inattentive responses can substantially
deteriorate the fit of theoretical models, produce spurious variance, and
overall jeopardize the validity of research results.

Due to the suspected presence of inattentive respodents, I employ my robust
estimator with tuning constant choice $c=1.6$111111The results remain
qualitatively similar for different finite choices of $c$. to estimate the
polychoric correlation matrix of the twelve items in the neuroticism scale.
The results of the remaining two scales are qualitatively similar. I then
conduct a confirmatory factor analysis (which is a special case of SEM) on
this estimated correlation matrix. I repeat this exercise for a non-robustly
estimated correlation matrix, namely the Pearson sample correlation matrix,
which is the standard way of conducting factor analyses with questionnaire
data. Conducting the analysis with a polychoric correlation matrix estimated
via MLE yields similar results as with a sample correlation matrix; I
therefore omit the former from the reported results.

###  5.2 Results of correlation matrix estimation

Figure 4: Difference between absolute estimates for the robustly estimated
polychoric correlation matrix and the Pearson sample correlation matrix of the
neuroticism scale, using the data of Arias et al., (2020). The numerical
estimates of each method are provided in Table C.1 in Appendix C. The
individual adjective items are “calm” (N1_P), “angry” (N1_N), “relaxed”
(N2_P), “tense” (N2_N), “at ease” (N3_P), “nervous” (N3_N), “not envious”
(N4_P), “envious” (N4_N), “stable” (N5_P), “unstable” (N5_N), “contented”
(N6_P), and “discontented” (N6_N). For the item naming given in parentheses,
items with identical identifier (the integer after the first “N”) are polar
opposites, where a last character “P” refers to the positive opposite and “N”
to the negative opposite.

Figure 4 visualizes the absolute differences between the robustly estimated
polychoric correlation matrix and sample correlation matrix of the neuroticism
scale. For all individual unique item pairs, the robust method estimates a
stronger correlation coefficient than sample correlation. The differences in
absolute estimates on average amount to 0.13, ranging from only marginally
larger than zero to a substantive 0.36. For correlations between polar
opposite adjectives, the average absolute difference between the robust method
and sample correlation is 0.20. The fact that a robust method consistently
yields stronger correlation estimates than a non-robust method, particularly
between polar opposite adjectives, is indicative of the presence of negative
leverage points, which drag negative correlational estimates towards zero,
that is, they attenuate the estimated strength of correlation. Here, such
negative leverage points could be the responses of inattentive participants
who report agreement or disagreement to both items in item pairs that are
designed to be negatively correlated. For instance, recall that it is
implausible that an attentive respondent would choose to agree (or disagree)
to both adjectives in the pair “envious” and “not envious” (cf. Arias et al.,,
2020). If sufficiently many such respondents exist, then the presumably
strongly negative correlation between these two opposite adjectives will be
estimated to be weaker than it actually is.

| Sample cor | MLE | Robust
---|---|---|---
Parameter | Estimate | SE | Estimate | SE | Estimate | SE
$\varrho$ | $-0.562$ | 0.031 | $-0.618$ | 0.025 | $-0.925$ | 0.062
$a_{1}$ |  |  | $-1.370$ | 0.061 | $-1.570$ | 0.276
$a_{2}$ |  |  | $-0.476$ | 0.043 | $-0.560$ | 0.203
$a_{3}$ |  |  | $0.121$ | 0.042 | $0.109$ | 0.187
$a_{4}$ |  |  | $1.060$ | 0.054 | $1.080$ | 0.105
$b_{1}$ |  |  | $-0.857$ | 0.049 | $-0.905$ | 0.073
$b_{2}$ |  |  | $-0.004$ | 0.041 | $-0.040$ | 0.091
$b_{3}$ |  |  | $0.608$ | 0.045 | $0.640$ | 0.364
$b_{4}$ |  |  | $1.580$ | 0.071 | $1.171$ | 0.811
Table 2: Parameter estimates with standard errors (SEs) for the correlation
between the neuroticism adjective pair “envious” and “not envious” in the data
of Arias et al., (2020), using sample correlation, polychoric correlation
fitted by maximum likelihood, and polychoric correlation fitted by my robust
estimator with tuning constant $c=1.6$. Each adjective item has five Likert-
type answer categories. The sample correlation coefficient does not model
thresholds, hence no estimates can be reported for them.

To further investigate the presence of inattentive respondents who attenuate
correlational estimates, I study in detail the adjective pair “not envious”
and “envious”, which featured the largest discrepancy between the non-robust
and robust estimates in Figure 4, with an absolute difference of 0.36. The
results are summarized in Table 2. The maximum likelihood estimate of $-0.618$
and sample correlation estimate of $-0.562$ for the correlation between these
two items seem remarkably weak considering that the two adjectives in question
are polar opposites. In contrast, its robust estimate estimate is given by
$-0.925$, which seems much more in line with what one would expect if all
participants responded accurately and attentively (cf. Arias et al.,, 2020).

Figure 5: Dot plot of cells for the neuroticism item adjective pair “envious”
and “not envious” in the data of Arias et al., (2020), where each item has
five Likert-type response options, anchored by “very inaccurate” (= 1) and
“very accurate” (= 5). Each dot’s size is proportional to the relative
empirical frequency of its associated cell, $\widehat{f}_{N}(\bm{z})$, whereas
its color varies by the value of the cell’s Pearson residual,
$\widehat{f}_{N}(\bm{z})/p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)$, at
robust parameter estimate with tuning constant $c=1.6$: The darker in blue a
dot, the larger the value of the Pearson residual of its associated cell. The
color of cells that could not be fitted well is fixed to red, where I deem a
fit poor if the Pearson residual exceeds value 5.48 (which is the value of the
smallest Pearson residual that is substantially larger than ideal value 1).
The data visualized here can be found in Table C.2 in Appendix C.

To study the potential presence of inattentive responses in each response cell
$\bm{z}=(x,y)\in\\{1,\dots,5\\}^{2}$ of the item pair “envious” and “not
envious”, Figure 5 visualizes the empirical relative frequencies,
$\widehat{f}_{N}(\bm{z})$, through dot size, as well as the associated Pearson
residual at the robust estimate,
$\widehat{f}_{N}(\bm{z})\Big{/}p_{\bm{z}}\left(\widehat{\bm{\theta}}_{N}\right)$,
through dot color (the darker the blue shade, the larger). Importantly, the
color of cells whose Pearson residual exceed 5.48 has been fixed to
red.121212The truncation of the color gradient in Figure 5 prevents that the
color gradient is dominated by single cells with extreme Pearson residuals,
which would blur the distinction between well fitted and poorly fitted cells.
This truncation value is equal to the value of the smallest Pearson residual
that is substantially larger than the ideal value 1. I consider cells whose
Pearson residual exceeds this truncation value to have a poor fit at the
polychoric model. This applies to a total of 12 cells, some of which have
enormous Pearson residuals. The Pearson residuals of the remaining 13 cells
are reasonably close to ideal value 1, ranging from 0.65 to 1.42 with average
0.89. The Pearson residuals as well as relative empirical frequencies of all
cells can be found in Table C.2 in the Appendix. It stands out that all poorly
fitted cells are those whose responses might be viewed as inconsistent.
Indeed, response cells $(x,y)=(1,1),(1,2),(2,1),(1,2)$ indicate that a
participant reports that neither “envious” nor “not envious” characterizes
them accurately, which are mutually contradicting responses, while for
response cells $(x,y)=(4,4),(4,5),(5,4),(5,5)$ both adjectives characterize
them accurately, which is again contradicting. As discussed previously, such
responses are likely due to inattentiveness. The robust estimator suggests
that such responses cannot be fitted well by the polychoric model and
subsequently downweighs their influence in the estimation procedure by mapping
their Pearson residual with the linear part of the $\rho(\cdot)$ function in
(9). Notably, also cells $(x,y)=(1,3),(3,1),(3,5),(5,3)$ are classified as
poorly fitted. These responses report (dis)agreement to one opposite
adjective, while being neutral about the other opposite. It is beyond the
scope of this paper to assess whether such response patterns are indicative of
inattentive responding, but the robust estimator suggests that such responses
at least cannot be fitted well by the polychoric model with the data of Arias
et al., (2020).

|  | “Envious”
---|---|---
|  $X$ \$Y$ | 1 | 2 | 3 | 4 | 5
“Not Envious” | 1 | $<0.0001$ | $0.0002$ | 0.8429 | 0.8274 | 0.8429
2 | $<0.0001$ | 0.0033 | 0.8429 | 0.9146 | 0.8274
3 | 0.8274 | 0.9896 | 0.8429 | 0.8429 | 0.8429
4 | 0.8429 | 0.9079 | 0.8863 | 0.5251 | $0.0003$
5 | 0.8429 | 0.8429 | 0.5251 | $<0.0001$ | $<0.0001$
Table 3: $p$-values, adjusted for multiple comparisons by the procedure of
Benjamini & Hochberg, (1995), of the cellwise goodness-of-fit test in
Corollary 3.1 for the neuroticism adjective pair “envious” and “not envious”
in the data of Arias et al., (2020), where each item has five Likert-type
response options, anchored by “1 = very inaccurate” and “5 = very accurate”.
The test statistics were computed with robust estimates using tuning constant
$c=1.6$. Cells in boldface are those for which the null hypothesis of unit
Pearson residual is rejected at significance level $\alpha=0.001$ in favor of
the alternative of it being larger than one.

Next, I perform the goodness-of-fit test derived in Corollary 3.1 for each
response cell to assess for which cells the polychoric model achieves a
statistically significantly poor fit in the “not envious”–“envious” item pair.
Table 3 presents the $p$-values for the hypothesis test in (11), adjusted for
multiple comparisons via the procedure of Benjamini & Hochberg, (1995). Values
for which the null hypothesis is rejected at significance level $\alpha=0.001$
are in boldface. This choice of significance level is deliberately extremely
conservative because the literature on inattentive responding recommends
overwhelming evidence in favor of inattention before one should label
responses as such (cf. Huang et al.,, 2012). At this significance level, the
null hypothesis of a good fit is rejected for six cells, namely
$(x,y)=(1,1),(2,1),(1,2),(5,4),(4,5),(5,5)$. These six cells comprise 5.52% of
the entire sample. As discussed in the previous section, it seems likely that
these responses are due to inattention because of inconsistent and
contradictory responding. Either way, my test offers strong empirical evidence
that these cells are outlying in the sense that they cannot be fitted well by
the polychoric model and therefore lead to deteriorated model fit. This is
consistent with Arias et al., (2020), who find that even a relatively small
proportion of inconsistent responses can drastically reduce a model’s fit. In
their analyses, they estimate that 6% of all respondents in the neuroticism
scale have been inattentive. Yet, albeit similar, I emphasize that the
estimate of 5.52% can be, if at all, understood as a lower bound for the
proportion of inattentive responding because of the extremely conservative
significance level I chose for my analyses. For instance, the null hypothesis
of good model fit was not rejected for the seemingly inconsistent response
cell $(x,y)=(2,2)$ (relative empirical frequency of about 4%) with a $p$-value
of approximately 0.003, but would have been rejected at a slightly more
liberal significance level. In addition, it is worth noting that the null
hypothesis was also not rejected for one more seemingly inconsistent response
cell, namely $(x,y)=(4,4)$, despite a relatively large Pearson residual of
12.66. This non-rejection is likely due to low statistical power stemming from
a small empirical frequency of this cell, since it only counted 14 responses
(out of 725). Similar reasoning applies to the remaining four cells that were
highlighted in red in Figure 4 but for which the null hypothesis of good fit
was not rejected, namely those who indicate (dis)agreement to one adjective,
while being neutral about its opposite. These four cells,
$(x,y)=(1,3),(3,1),(3,5),(5,3)$, only count empirical frequencies of 2, 4, 2,
and 4, respectively.

Overall, leveraging the proposed robust estimator, I find strong evidence for
the presence of inattentive respondents in the data of Arias et al., (2020).
While they substantially affect the correlational estimate of Pearson sample
correlation, amounting to about $-0.56$, which is much weaker than one would
expect for polar opposite items, my robust estimator can withstand their
influence with an estimate of about $-0.93$ and also identify them by means of
the proposed test. Similar conclusions follow by repeating this analysis for
different item pairs.

###  5.3 Results of structural equation modeling

In this section, I perform confirmatory factor analyses on the data of Arias
et al., (2020) in the neuroticism scale. That is, I fit a factor model as in
(2) with a single factor to each of the two estimated correlation matrices.
Table 4 presents the estimates factor loadings for each estimator. The
estimated loadings of the robust estimator are consistently higher than those
of the non-robust sample correlation estimator, indicating greater internal
consistency. Indeed, the neuroticism factor of the robust analysis can explain
53% of the variation in the data, whereas the factor of the non-robust
analysis can only explain 40%. In addition, the robust analysis yields higher
reliability than the non-robust analysis, as can be seen by Cronbach-$\alpha$
estimates of 0.89 and 0.93, respectively.

Item | Sample corr | Robust polycor
---|---|---
N1_P | 0.70 | 0.80
N1_N | 0.56 | 0.66
N2_P | 0.76 | 0.86
N2_N | 0.68 | 0.78
N3_P | 0.77 | 0.88
N3_N | 0.66 | 0.74
N4_P | 0.35 | 0.46
N4_N | 0.46 | 0.54
N5_P | 0.69 | 0.77
N5_N | 0.67 | 0.73
N6_P | 0.57 | 0.66
N6_N | 0.64 | 0.71
Proportion variance | 0.40 | 0.53
Cronbach’s $\alpha$ | 0.89 | 0.93
Table 4: Factor loadings estimates in the unidimensional confirmatory factor
analysis of the neuroticism scale, using the data of Arias et al., (2020).

The fact that a robust estimator yields a more internally consistent factor
structure than a non-robust estimator is indicative of the presence of
irregular data points, namely the inconsistent (and presumably inattentive)
responses the robust estimator identified in the previous section and whose
adverse influence it can withstand by down-weighting them in the estimation
procedure.

This empirical application demonstrated how the proposed estimator can be used
to robustify SEMs against inattentive respondents. I find compelling evidence
for the presence of inattentive responding, and the proposed robust estimator
can not only identify them, but also account for their presence to obtain a
good model fit.

##  6 Conclusion

In this paper, I develop a novel estimator for models of categorical variables
that is designed to be robust to misspecification of such models, which is the
first of its kind and can be thought of as an analogue to robust
$M$-estimation for non-continuous variables. The estimator is shown to be
consistent and asymptotically normally distributed, and possesses attractive
properties with respect to robustness and computation. Crucially, the
estimator makes no assumption whatsoever on the degree, magnitude, or type of
misspecification. If misspecification is absent, the estimator is
asymptotically equivalent to maximum likelihood estimation (MLE), but more
robust than MLE in the presence of misspecification. In addition, I develop a
novel diagnostic test that can test if a given categorical observation can be
fitted well by the presumed model, allowing one to trace back potential
sources of model misspecification. The methodology proposed in this paper is
implemented in the free open source package robcat (Welz,, 2024) in the
statistical programming environment R, although it is primarily developed in
C++ to maximize speed and computational performance.

I verify the enhanced robustness and theoretical properties of the novel
estimator in simulation studies and demonstrate its practical usefulness in an
empirical application on structural equation modeling of questionnaire
responses to a Big-5 administration. I find compelling evidence for the
presence of inattentive respondents. For instance, in a rating item pair with
polar opposite content where a strong negative correlation is expected, the
robust estimator yields a correlational estimate of $-0.93$, whereas non-
robust estimators yield only $-0.56$ to $-0.62$; it follows that the robust
estimate is more in line with the literature on the corresponding scale.
Utilizing the proposed diagnostic test, I argue that the lower-than-expected
estimates of the non-robust method are likely due to a few possibly
inattentive participants who gave mutually contradictory responses, while the
robust estimator can resist their influence. In addition, conducting a
confirmatory factor analyses with the robust estimator yielded more core
consistent factor loadings, higher explained variance, and greater scale
reliability than commonly employed estimation methods.

Although the proposed estimator is particularly useful in models for
questionnaire responses, I stress that it can be applied to any model of
categorical responses. Examples include models from item response theory,
counting processes, or categorical regression, thereby potentially giving rise
to a new research line. I leave these exciting avenues to further research.

## References

  * Agresti, (2010) Agresti, A. (2010).  Analysis of ordinal categorical data (2nd ed.).  John Wiley & Sons.
  * Agresti, (2012) Agresti, A. (2012).  Categorical data analysis (3rd ed.).  John Wiley & Sons.
  * Alfons & Welz, (2024) Alfons, A. & Welz, M. (2024).  Open science perspectives on machine learning for the identification of careless responding: A new hope or phantom menace?  Social and Personality Psychology Compass.  https://doi.org/10.1111/spc3.12941.  In press.
  * Almlund et al., (2011) Almlund, M., Duckworth, A. L., Heckman, J., & Kautz, T. (2011).  Personality psychology and economics.  Handbook of The Economics of Education, volume 4 of Handbook of the Economics of Education, 1–181. Elsevier.  https://doi.org/10.1016/B978-0-444-53444-6.00001-8
  * Andersen, (1972) Andersen, E. B. (1972).  The numerical solution of a set of conditional estimation equations.  Journal of the Royal Statistical Society: Series B (Methodological), 34(1), 42–54.  https://doi.org/10.1111/j.2517-6161.1972.tb00887.x
  * Andrews et al., (2017) Andrews, I., Gentzkow, M., & Shapiro, J. M. (2017).  Measuring the Sensitivity of Parameter Estimates to Estimation Moments*.  Quarterly Journal of Economics, 132(4), 1553–1592.  https://doi.org/10.1093/qje/qjx023
  * Arias et al., (2020) Arias, V. B., Garrido, L., Jenaro, C., Martínez-Molina, A., & Arias, B. (2020).  A little garbage in, lots of garbage out: Assessing the impact of careless responding in personality survey data.  Behavior Research Methods, 52(6), 2489–2505.  https://doi.org/https://doi.org/10.3758/s13428-020-01401-8
  * Armstrong & Kolesár, (2021) Armstrong, T. B. & Kolesár, M. (2021).  Sensitivity analysis using approximate moment condition models.  Quantitative Economics, 12(1), 77–108.  https://doi.org/https://doi.org/10.3982/QE1609
  * Benjamini & Hochberg, (1995) Benjamini, Y. & Hochberg, Y. (1995).  Controlling the false discovery rate: A practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society: Series B (Methodological), 57(1), 289–300.  https://doi.org/https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
  * Bond & Lang, (2019) Bond, T. N. & Lang, K. (2019).  The sad truth about happiness scales.  Journal of Political Economy, 127(4), 1629–1640.  https://doi.org/10.1086/701679
  * Bonhomme & Weidner, (2022) Bonhomme, S. & Weidner, M. (2022).  Minimizing sensitivity to model misspecification.  Quantitative Economics, 13(3), 907–954.  https://doi.org/https://doi.org/10.3982/QE1930
  * Borghans et al., (2008) Borghans, L., Duckworth, A. L., Heckman, J., & ter Weel, B. (2008).  The economics and psychology of personality traits.  Journal of Human Resources, 43(4), 972–1059.  https://doi.org/10.3368/jhr.43.4.972
  * Bowling et al., (2016) Bowling, N. A., Huang, J. L., Bragg, C. B., Khazon, S., Liu, M., & Blackmore, C. E. (2016).  Who cares and who is careless? Insufficient effort responding as a reflection of respondent personality.  Journal of Personality and Social Psychology, 111(2), 218–229.  https://psycnet.apa.org/doi/10.1037/pspp0000085
  * Bowling et al., (2023) Bowling, N. A., Huang, J. L., Brower, C. K., & Bragg, C. B. (2023).  The quick and the careless: The construct validity of page time as a measure of insufficient effort responding to surveys.  Organizational Research Methods, 26(2), 323–352.  https://doi.org/10.1177/10944281211056520
  * Bugni & Ura, (2019) Bugni, F. A. & Ura, T. (2019).  Inference in dynamic discrete choice problems under local misspecification.  Quantitative Economics, 10(1), 67–103.  https://doi.org/https://doi.org/10.3982/QE917
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. K., & Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  Econometrics Journal, 21(1), C1–C68.  https://doi.org/10.1111/ectj.12097
  * Chernozhukov et al., (2022) Chernozhukov, V., Escanciano, J. C., Ichimura, H., Newey, W. K., & Robins, J. M. (2022).  Locally robust semiparametric estimation.  Econometrica, 90(4), 1501–1535.  https://doi.org/https://doi.org/10.3982/ECTA16294
  * Credé, (2010) Credé, M. (2010).  Random responding as a threat to the validity of effect size estimates in correlational research.  Educational and Psychological Measurement, 70(4), 596–612.  https://doi.org/10.1177/0013164410366686
  * Cressie & Read, (1984) Cressie, N. & Read, T. R. (1984).  Multinomial goodness-of-fit tests.  Journal of the Royal Statistical Society: Series B (Methodological), 46(3), 440–464.
  * Cronbach, (1951) Cronbach, L. J. (1951).  Coefficient alpha and the internal structure of tests.  Psychometrika, 16, 297–334.  https://doi.org/10.1007/BF02310555
  * Curran, (2016) Curran, P. G. (2016).  Methods for the detection of carelessly invalid responses in survey data.  Journal of Experimental Social Psychology, 66, 4–19.  https://doi.org/10.1016/j.jesp.2015.07.006
  * Daley & Vere-Jones, (2003) Daley, D. J. & Vere-Jones, D. (2003).  An Introduction to the Theory of Point Processes (2nd ed.).  Probability and Its Applications. Springer.  https://doi.org/10.1007/b97277
  * Drasgow et al., (1985) Drasgow, F., Levine, M. V., & Williams, E. A. (1985).  Appropriateness measurement with polychotomous item response models and standardized indices.  British Journal of Mathematical and Statistical Psychology, 38(1), 67–86.  https://doi.org/10.1111/j.2044-8317.1985.tb00817.x
  * Eddelbuettel, (2013) Eddelbuettel, D. (2013).  Seamless R and C++ Integration with Rcpp.  Springer.  https://doi.org/10.1007/978-1-4614-6868-4.  ISBN 978-1-4614-6867-7
  * Foldnes & Grønneberg, (2022) Foldnes, N. & Grønneberg, S. (2022).  The sensitivity of structural equation modeling with ordinal data to underlying non-normality and observed distributional forms.  Psychological Methods, 27(4), 541–567.  https://doi.org/10.1037/met0000385
  * Gadermann et al., (2012) Gadermann, A. M., Guhn, M., & Zumbo, B. D. (2012).  Estimating ordinal reliability for likert-type and ordinal item response data: A conceptual, empirical, and practical guid.  Practical Assessment, Research, and Evaluation, 17(1), 1–13.  https://doi.org/10.7275/n560-j767
  * Garrido et al., (2013) Garrido, L. E., Abad, F. J., & Ponsoda, V. (2013).  A new look at Horn’s parallel analysis with ordinal variables.  Psychological Methods, 18(4), 454–474.  https://doi.org/10.1037/a0030005
  * Goldberg, (1992) Goldberg, L. R. (1992).  The development of markers for the Big-Five factor structure.  Psychological Assessment, 4(1), 26–42.  https://doi.org/https://doi.org/10.1037/1040-3590.4.1.26
  * Götze, (1991) Götze, F. (1991).  On the rate of convergence in the multivariate CLT.  Annals of Probability, 19(2), 724–739.  https://doi.org/10.1214/aop/1176990448
  * Hampel, (1974) Hampel, F. R. (1974).  The influence curve and its role in robust estimation.  Journal of the American Statistical Association, 69(346), 383–393.  https://doi.org/10.1080/01621459.1974.10482962
  * Hampel et al., (1986) Hampel, F. R., Ronchetti, E. M., Rousseeuw, P., & Stahel, W. A. (1986).  Robust statistics: The approach based on influence functions.  Wiley series in probability and mathematical statistics. Wiley.
  * He & Simpson, (1993) He, X. & Simpson, D. G. (1993).  Lower bounds for contamination bias: Globally minimax versus locally linear estimation.  Annals of Statistics, 21(1), 314–337.  https://doi.org/10.1214/aos/1176349028
  * Heckman et al., (2006) Heckman, J. J., Stixrud, J., & Urzua, S. (2006).  The effects of cognitive and noncognitive abilities on labor market outcomes and social behavior.  Journal of Labor Economics, 24(3), 411–482.  https://doi.org/10.1086/504455
  * Holgado-Tello et al., (2010) Holgado-Tello, F. P., Chacón-Moscoso, S., Barbero-García, I., & Vila-Abad, E. (2010).  Polychoric versus pearson correlations in exploratory and confirmatory factor analysis of ordinal variables.  Quality & Quantity, 44, 153–166.  https://doi.org/10.1007/s11135-008-9190-y
  * (35) Huang, J. L., Bowling, N. A., Liu, M., & Li, Y. (2015a).  Detecting insufficient effort responding with an infrequency scale: Evaluating validity and participant reactions.  Journal of Business and Psychology, 30(2), 299–311.  https://doi.org/10.1007/s10869-014-9357-6
  * Huang et al., (2012) Huang, J. L., Curran, P. G., Keeney, J., Poposki, E. M., & DeShon, R. P. (2012).  Detecting and deterring insufficient effort responding to surveys.  Journal of Business and Psychology, 27(1), 99–114.  https://doi.org/10.1007/s10869-011-9231-8
  * (37) Huang, J. L., Liu, M., & Bowling, N. A. (2015b).  Insufficient effort responding: Examining an insidious confound in survey data.  Journal of Applied Psychology, 100(3), 828–845.  https://doi.org/10.1037/a0038510
  * Huber, (1964) Huber, P. J. (1964).  Robust estimation of a location parameter.  Annals of Mathematical Statistics, 35(1), 73–101.  https://doi.org/10.1214/aoms/1177703732
  * Huber & Ronchetti, (2009) Huber, P. J. & Ronchetti, E. M. (2009).  Robust Statistics.  Wiley.  https://doi.org/10.1002/9780470434697
  * Kam & Meyer, (2015) Kam, C. C. S. & Meyer, J. P. (2015).  How careless responding and acquiescence response bias can influence construct dimensionality: The case of job satisfaction.  Organizational Research Methods, 18(3), 512–541.  https://doi.org/10.1177/1094428115571894
  * Kim et al., (2018) Kim, D. S., Reise, S. P., & Bentler, P. M. (2018).  Identifying aberrant data in structural equation models with IRLS-ADF.  Structural Equation Modeling: A Multidisciplinary Journal, 25(3), 343–358.  https://doi.org/10.1080/10705511.2017.1379881
  * Kitamura et al., (2013) Kitamura, Y., Otsu, T., & Evdokimov, K. (2013).  Robustness, infinitesimal neighborhoods, and moment restrictions.  Econometrica, 81(3), 1185–1201.  https://doi.org/https://doi.org/10.3982/ECTA8617
  * Lai & Green, (2016) Lai, K. & Green, S. B. (2016).  The problem with having two watches: Assessment of fit when rmsea and cfi disagree.  Multivariate Behavioral Research, 51(2-3), 220–239.  https://doi.org/10.1080/00273171.2015.1134306
  * Lindsay, (1994) Lindsay, B. G. (1994).  Efficiency versus robustness: The case for minimum Hellinger distance and related methods.  Annals of Statistics, 22(2), 1081–1114.  https://doi.org/10.1214/aos/1176325512
  * Mair, (2018) Mair, P. (2018).  Modern Psychometrics with R.  Springer.  https://doi.org/10.1007/978-3-319-93177-7
  * Maniaci & Rogge, (2014) Maniaci, M. R. & Rogge, R. D. (2014).  Caring about carelessness: Participant inattention and its effects on research.  Journal of Research in Personality, 48, 61–83.  https://doi.org/10.1016/j.jrp.2013.09.008
  * Mardia et al., (1979) Mardia, K. V., Kent, J. T., & Bibby, J. M. (1979).  Multivariate Analysis (1st ed.).  Series in Probability and Mathematical Statistics. Academic Press.
  * Maronna et al., (2018) Maronna, R. A., Martin, R. D., Yohai, V. J., & Salibián-Barrera, M. (2018).  Robust Statistics: Theory and Methods (2nd ed.).  John Wiley & Sons.
  * McDonald, (1970) McDonald, R. P. (1970).  The theoretical foundations of principal factor analysis, canonical factor analysis, and alpha factor analysis.  British Journal of Mathematical and Statistical Psychology, 23(1), 1–21.  https://doi.org/10.1111/j.2044-8317.1970.tb00432.x
  * McGrath et al., (2010) McGrath, R. E., Mitchell, M., Kim, B. H., & Hough, L. (2010).  Evidence for response bias as a source of error variance in applied assessment.  Psychological Bulletin, 136(3), 450–470.  https://doi.org/10.1037/a0019216
  * Meade & Craig, (2012) Meade, A. W. & Craig, S. B. (2012).  Identifying careless responses in survey data.  Psychological Methods, 17(3), 437–455.  https://psycnet.apa.org/doi/10.1037/a0028085
  * Newey, (1985) Newey, W. K. (1985).  Generalized method of moments specification testing.  Journal of Econometrics, 29(3), 229–256.  https://doi.org/https://doi.org/10.1016/0304-4076(85)90154-X
  * Newey, (1991) Newey, W. K. (1991).  Uniform convergence in probability and stochastic equicontinuity.  Econometrica, 59(4), 1161–1167.  https://doi.org/10.2307/2938179
  * Neyman & Scott, (1948) Neyman, J. & Scott, E. L. (1948).  Consistent estimates based on partially consistent observations.  Econometrica, 16(1), 1–32.  https://doi.org/10.2307/1914288
  * Olsson, (1979) Olsson, U. (1979).  Maximum likelihood estimation of the polychoric correlation coefficient.  Psychometrika, 44(4), 443–460.  https://doi.org/10.1007/BF02296207
  * Osborne-Groves, (2005) Osborne-Groves, M. (2005).  How important is your personality? Labor market returns to personality for women in the US and UK.  Journal of Economic Psychology, 26(6), 827–841.  https://doi.org/10.1016/j.joep.2005.03.001
  * Pearson, (1901) Pearson, K. (1901).  I. mathematical contributions to the theory of evolution, vii: On the correlation of characters not quantitatively measurable.  Philosophical Transactions of the Royal Society of London, Series A, 195(262-273), 1–47.  https://doi.org/10.1098/rsta.1900.0022
  * Pearson & Pearson, (1922) Pearson, K. & Pearson, E. S. (1922).  On polychoric coefficients of correlation.  Biometrika, 14(1–2), 127–156.  https://doi.org/10.1093/biomet/14.1-2.127
  * Pison et al., (2003) Pison, G., Rousseeuw, P. J., Filzmoser, P., & Croux, C. (2003).  Robust factor analysis.  Journal of Multivariate Analysis, 84(1), 145–172.  https://doi.org/10.1016/S0047-259X(02)00007-6
  * Rasch, (1960) Rasch, G. (1960).  Probabilistic models for some intelligence and attainment tests.
  * Raymaekers & Rousseeuw, (2021) Raymaekers, J. & Rousseeuw, P. J. (2021).  Fast robust correlation for high-dimensional data.  Technometrics, 63(2), 184–198.  https://doi.org/10.1080/00401706.2019.1677270
  * Rhemtulla et al., (2012) Rhemtulla, M., Brosseau-Liard, P. É., & Savalei, V. (2012).  When can categorical variables be treated as continuous? a comparison of robust continuous and categorical SEM estimation methods under suboptimal conditions.  Psychological Methods, 17(3), 354–373.  https://doi.org/10.1037/a0029315
  * Robinson, (2018) Robinson, M. A. (2018).  Using multi-item psychometric scales for research and practice in human resource management.  Human Resource Management, 57(3), 739–750.  https://doi.org/10.1002/hrm.21852
  * Rousseeuw, (1985) Rousseeuw, P. J. (1985).  Multivariate estimation with high breakdown point.  Mathematical Statistics and Applications, 283–297. Reidel.
  * Ruckstuhl & Welsh, (2001) Ruckstuhl, A. F. & Welsh, A. H. (2001).  Robust fitting of the binomial model.  Annals of Statistics, 29(4), 1117–1136.  https://doi.org/10.1214/aos/1013699996
  * Rudin, (1976) Rudin, W. (1976).  Principles of Mathematical Analysis (3rd ed.).  McGraw-Hill.
  * Rust et al., (2020) Rust, J., Golombok, S., & Stillwell, D. (2020).  Modern psychometrics: The science of psychological assessment (4th ed.).  Routledge.
  * Schennach, (2022) Schennach, S. (2022).  Measurement systems.  Journal of Economic Literature, 60(4), 1223–1263.  https://doi.org/10.1257/jel.20211355
  * Schmitt & Stults, (1985) Schmitt, N. & Stults, D. M. (1985).  Factors defined by negatively keyed items: The result of careless respondents?  Applied Psychological Measurement, 9, 367–373.  https://doi.org/10.1177/014662168500900405
  * Schroeders et al., (2022) Schroeders, U., Schmidt, C., & Gnambs, T. (2022).  Detecting careless responding in survey data using stochastic gradient boosting.  Educational and Psychological Measurement, 82(1), 29–56.  https://doi.org/10.1177/00131644211004708
  * Simpson, (1987) Simpson, D. G. (1987).  Minimum hellinger distance estimation for the analysis of count data.  Journal of the American Statistical Association, 82(399), 802–807.  https://doi.org/10.1080/01621459.1987.10478501
  * Stantcheva, (2022) Stantcheva, S. (2022).  How to run surveys: A guide to creating your own identifying variation and revealing the invisible.  Working Paper 30527, National Bureau of Economic Research.  https://doi.org/10.3386/w30527
  * Van de Geer, (2000) Van de Geer, S. (2000).  Empirical Processes in $M$-estimation.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Van der Vaart, (1998) Van der Vaart, A. W. (1998).  Asymptotic Statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Vershynin, (2018) Vershynin, R. (2018).  High-dimensional probability: An introduction with applications in data science.  Cambridge University Press.  https://doi.org/10.1017/9781108231596
  * Victoria-Feser & Ronchetti, (1997) Victoria-Feser, M.-P. & Ronchetti, E. (1997).  Robust estimation for grouped data.  Journal of the American Statistical Association, 92(437), 333–340.  https://doi.org/10.1080/01621459.1997.10473631
  * Ward & Meade, (2023) Ward, M. & Meade, A. W. (2023).  Dealing with careless responding in survey data: Prevention, identification, and recommended best practices.  Annual Review of Psychology, 74(1), 577–596.  https://doi.org/10.1146/annurev-psych-040422-045007
  * Welz, (2024) Welz, M. (2024).  robcat: Robust Categorical Data Analysis.  https://github.com/mwelz/robcat.  R package version 0.0.1
  * Welz & Alfons, (2023) Welz, M. & Alfons, A. (2023).  I don’t care anymore: Identifying the onset of careless responding.  https://doi.org/10.48550/arXiv.2303.07167.  arXiv:2303.07167
  * Woods, (2006) Woods, C. M. (2006).  Careless responding to reverse-worded items: Implications for confirmatory factor analysis.  Journal of Psychopathology and Behavioral Assessment, 28(3), 186–191.  https://psycnet.apa.org/doi/10.1007/s10862-005-9004-7

##  Appendix A Popular categorical models

In addition SEM presented in Section 2.2, this section describes two more
commonly used models for categorical models, namely the Rasch model from item
response theory and a Poisson counting process. I stress that these three
models (SEM, Rasch, Poisson) are designed for three different types of
categorical variables, namely ordinal ones, nominal ones, and counting
variables. The proposed robust estimator can be applied to fit all three
models.

######  Example 1 (Rasch model).

The Rasch model (Rasch,, 1960) is a fundamental model to model a test taker’s
probability to answer correctly to a test’s questions. Let the test comprise
$k$ questions and let $\bm{X}=(X_{1},\dots,X_{k})^{\top}$ be a binary vector
whose $j$-th element equals 1 if the $j$-th question was answered correctly,
and zero otherwise. It follows that the sample space is given by
$\bm{\mathcal{X}}=\\{0,1\\}^{k}$. Further, denote by $S=\sum_{j=1}^{k}X_{j}$
the number of correctly answered questions, known as the score. Since the
unconditional Rasch model suffers from the incidental parameter problem,iiiThe
unconditional Rasch model is defined by the individual-specific probability
$\mathbb{P}\left[\bm{X}_{i}=\bm{x}\right]=\frac{\exp\big{(}\sum_{j=1}^{k}x_{j}(\alpha_{i}-\theta_{j})\big{)}}{\prod_{j=1}^{k}(1+\exp(\alpha_{i}-\theta_{j}))},$
where $\theta_{j}$ parametrizes the difficulty of question $j$, whereas
$\alpha_{i}$ parametrizes the ability of test taker $i=1,\dots,N$, through
fixed effects. Joint maximum likelihood estimates of the ability and
difficulty parameters will not be consistent since the number of parameters
grows with the sample size $N$, which is an instance of the incidental
parameter problem (Neyman & Scott,, 1948).  one instead (Andersen,, 1972)
works with the conditional probability for $\bm{x}\in\bm{\mathcal{X}}$ given a
score $s\in\\{0,\dots,k\\}$, being

$p_{\bm{z}}\left(\bm{\theta}\right)=\mathbb{P}_{\bm{\theta}}\left[\bm{X}=\bm{x}\
|\
S=s\right]=\frac{\exp(-\bm{x}^{\top}\bm{\theta})}{\sum_{\bm{y}\in\bm{\mathcal{D}}(s)}\exp(-\bm{y}^{\top}\bm{\theta})},$

where $\bm{z}=(\bm{x}^{\top},s)^{\top}$, $\bm{\theta}\in\mathbb{R}^{k}$ is a
$k$-vector of question difficulties, and
$\bm{\mathcal{D}}(s)=\left\\{\bm{x}\in\bm{\mathcal{X}}:\sum_{j=1}^{k}x_{j}=s\right\\}$
denotes the set of response vector with score $s$.

######  Example 2 (Poisson process).

Counting processes model how often a certain event occurs in a given time
period, which is a discrete outcome. Assume that one has access to counting
observations in $k$ periods, where the $j$-th period is given by
$(a_{j},b_{j}]$ with known finite boundaries $a_{j}<b_{j}\leq a_{j+1}$, and
random variable $Z_{j}\in\mathbb{N}_{0}=\\{0,1,2,\dots\\}$ counts the number
of event in the $j$-th period, $j=1,\dots,k$. Consequently, the
$k$-dimensional random variable $\bm{Z}=(Z_{1},\dots,Z_{k})^{\top}$ holds the
counts per period and takes values in sample space
$\bm{\mathcal{Z}}=\mathbb{N}_{0}^{k}$.iiiiiiIn practice, the empirical support
is bounded, with the largest observed number of counts being the upper bound.
One of the most popular counting processes is the stationary Poisson point
process (e.g. Daley & Vere-Jones,, 2003, equation 2.1.1), which defines the
probability of $\bm{z}=(z_{1},\dots,z_{k})^{\top}\in\bm{\mathcal{Z}}$ as

$p_{\bm{z}}\left(\lambda\right)=\mathbb{P}_{\lambda}\left[\bm{Z}=\bm{z}\right]=\prod_{j=1}^{k}\frac{\big{(}\lambda(b_{j}-a_{j})\big{)}^{z_{j}}}{z_{j}!}\exp\big{(}-\lambda(b_{j}-a_{j})\big{)},$

where $\lambda>0$ is an intensity parameter.

##  Appendix B Additional theoretical results

The influence function (Hampel,, 1974) measures how an infinitesimally small
degree of misspecification affects an estimator. To construct an influence
function, one adopts a functional notation where an estimator
$\widehat{\bm{\theta}}_{N}$ is interpreted as a function of the distribution
of the data that it is calculated from. Hence, if the observed data are
distributed according to a density $f$, one writes
$\widehat{\bm{\theta}}_{N}=\widehat{\bm{\theta}}_{N}(f)$. Expressing an
estimator as a statistical functional of a distribution allows one to study
how it changes with respect to that distribution. This notion is formalized by
means of the influence function, which is a functional (Gâteaux) derivative of
the estimator with respect to the degree of model misspecification.

######  Definition 1 (Influence function; adapted from Hampel,, 1974).

Let $\widehat{\bm{\theta}}_{N}$ be an estimator that estimates a model
$\\{\bm{p}\left(\bm{\theta}\right):\bm{\theta}\in\bm{\Theta}\\}$ with finite
support $\bm{\mathcal{Z}}$. Evaluated at model density
$p_{\bm{z}}\left(\bm{\theta}\right),\theta\in\bm{\theta}$, the estimator’s
influence function at a data point $\bm{z}\in\bm{\mathcal{Z}}$ is given by

$\mathrm{IF}\left(\bm{z},\widehat{\bm{\theta}}_{N},\bm{p}\left(\bm{\theta}\right)\right)=\lim_{\varepsilon\downarrow
0}\frac{\widehat{\bm{\theta}}_{N}\big{(}(1-\varepsilon)p_{\bm{z}}\left(\bm{\theta}\right)+\varepsilon\Delta_{\bm{z}}\big{)}}{\varepsilon}=\frac{\partial}{\partial\varepsilon}\widehat{\bm{\theta}}_{N}\big{(}(1-\varepsilon)p_{\bm{z}}\left(\bm{\theta}\right)+\varepsilon\Delta_{\bm{z}}\big{)}\bigg{|}_{\varepsilon=0},$

where
$\bm{\mathcal{Z}}\ni\bm{y}\mapsto\Delta_{\bm{z}}(\bm{y})=\mathds{1}\left\\{\bm{z}=\bm{y}\right\\}$
is the point mass density at point $\bm{z}$.

Essential theory on $M$-estimation (e.g. Huber & Ronchetti,, 2009, Section
3.2) reveals the following influence function of the MLE (eq. 5) at point
$\bm{z}\in\bm{\mathcal{Z}}$ and density $\bm{p}\left(\bm{\theta}\right)$ is
given by

$\mathrm{IF}\left(\bm{z},\widehat{\bm{\theta}}_{N}^{\mathrm{\
MLE}},\bm{p}\left(\bm{\theta}\right)\right)=\bm{J}\left(\bm{\theta}\right)^{-1}\bm{s}_{\bm{z}}\left(\bm{\theta}\right),$
(B.1)

where

$\bm{J}\left(\bm{\theta}\right)=-\sum_{\bm{z}\in\bm{\mathcal{Z}}}\bm{Q}_{\bm{z}}(\bm{\theta})p_{\bm{z}}\left(\bm{\theta}\right)$

denotes the Fisher information matrix at $\bm{\theta}\in\bm{\Theta}$. The
following theorem derives the influence function of the proposed estimator
$\widehat{\bm{\theta}}_{N}$ in (8).

######  Theorem 3 (Influence function).

Grant Assumption Set A. Then, the influence function of estimator
$\widehat{\bm{\theta}}_{N}$ at cell $\bm{z}\in\bm{\mathcal{Z}}$ and true
density $\bm{p}\left(\bm{\theta}_{*}\right)$ is given by

$\displaystyle\mathrm{IF}\left(\bm{z},\widehat{\bm{\theta}}_{N},\bm{p}\left(\bm{\theta}_{*}\right)\right)=\begin{cases}\mathrm{IF}\left(\bm{z},\widehat{\bm{\theta}}_{N}^{\mathrm{\
MLE}},\bm{p}\left(\bm{\theta}_{*}\right)\right)&\text{ if }c>1,\\\
\bigg{[}\bm{J}\left(\bm{\theta}_{*}\right)-p_{\bm{z}}\left(\bm{\theta}_{*}\right)\bm{s}_{\bm{z}}\left(\bm{\theta}_{*}\right)\bm{s}_{\bm{z}}\left(\bm{\theta}_{*}\right)^{\top}\bigg{]}^{-1}\bm{s}_{\bm{z}}\left(\bm{\theta}_{*}\right)p_{\bm{z}}\left(\bm{\theta}_{*}\right)&\text{
if }c=1.\end{cases}$

This theorem’s proof is provided in Appendix D.4. As an example, Figure B.1
visualizes this influence function for a polychoric correlation model.

(a) $c>1$

(b) $c=1$

Figure B.1: Influence functions of polychoric correlation coefficient
$\varrho_{*}=0$ in the polychoric model with five Likert-type answer
categories and thresholds like in (13). The top plot visualizes the influence
function for tuning constant $c>1$, and the bottom plot for the most robust
choice, $c=1$; see Theorem 3.

The following section contains a number of technical lemmas that will be
useful in the proofs of the main theoretical results.

###  B.1 Some useful lemmas

######  Lemma B.1.

Under the Assumptions of Theorem 1, the sequence
$\left\\{{L}\left({\bm{\theta},\widehat{f}_{N}}\right)\right\\}_{N}$ is
equicontinuous on the parameter space $\bm{\Theta}\ni\bm{\theta}$.

Proof. I verify equicontinuity by its definition (e.g. Definition 7.22 in
Rudin,, 1976). Put
${L_{N}}\left({\bm{\theta}}\right)={L}\left({\bm{\theta},\widehat{f}_{N}}\right)$
and let

$\bm{g}_{N}(\bm{\theta})=\frac{\partial}{\partial\bm{\theta}}{L_{N}}\left({\bm{\theta}}\right)=-\sum_{\bm{z}\in\bm{\mathcal{Z}}}\bm{s}_{\bm{z}}\left(\bm{\theta}\right)f_{\varepsilon}(\bm{z}){w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}\Bigg{)}$

denote its gradient at $\bm{\theta}\in\bm{\Theta}$. By Assumption A.5, there
exists a universal constant $C>0$ such that $\|\bm{g}_{N}(\bm{\theta})\|<C$
for all $\bm{\theta}\in\bm{\Theta}$. Fix $\varepsilon>0$ and put
$\delta=\varepsilon/C$. Let $\bm{\theta}_{a},\bm{\theta}_{b}\in\bm{\Theta}$ be
vectors such that $\|\bm{\theta}_{a}-\bm{\theta}_{b}\|<\delta$ and the entire
line segment between these two vectors is contained in $\bm{\Theta}$. Then, by
the multivariate mean value theorem, there exists a $\bm{\theta}_{c}$ in this
line segment such that

$L_{N}(\bm{\theta}_{b})-L_{N}(\bm{\theta}_{a})=\bm{g}_{N}(\bm{\theta})^{\top}(\bm{\theta}_{b}-\bm{\theta}_{a}).$

Thus,

$\displaystyle|L_{N}(\bm{\theta}_{b})-L_{N}(\bm{\theta}_{a})|$
$\displaystyle=|\bm{g}_{N}(\bm{\theta})^{\top}(\bm{\theta}_{b}-\bm{\theta}_{a})|$
$\displaystyle\leq\|\bm{g}_{N}(\bm{\theta})\|\|\bm{\theta}_{b}-\bm{\theta}_{a}\|$
$\displaystyle\leq C\|\bm{\theta}_{b}-\bm{\theta}_{a}\|$
$\displaystyle<C\delta$ $\displaystyle=\varepsilon,$

where I have applied the Cauchy-Schwarz inequality in the second line. It
follows from the definition of equicontinuity that the sequence
$\left\\{{L_{N}}\left({\bm{\theta}}\right)\right\\}_{N}$ is equicontinuous on
$\bm{\Theta}$. $\blacksquare$

The following lemma can be seen as a multivariate special case of the classic
Berry-Esseen theorem (e.g. Theorem 2.1.3 in Vershynin,, 2018).

######  Lemma B.2.

For $i=1,\dots,N$, put
$\bm{J}_{i}=\big{(}J_{i,\bm{z}_{1}},J_{i,\bm{z}_{2}},\dots,J_{i,\bm{z}_{m}}\big{)}^{\top}=\big{(}\mathds{1}\left\\{\bm{Z}_{i}=\bm{z}_{1}\right\\},\mathds{1}\left\\{\bm{Z}_{i}=\bm{z}_{2}\right\\},\dots,\mathds{1}\left\\{\bm{Z}_{i}=\bm{z}_{m}\right\\}\big{)}^{\top}$
and
$\bm{X}_{i}=-\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\bm{W}(\bm{\theta}_{N}(\bm{t}))\left(\bm{J}_{i}-\bm{f}_{\varepsilon}\right)$,
where $\bm{\theta}_{N}(\bm{t})=\bm{\theta}_{0}+N^{-1/2}\bm{V}\bm{t}$ for a
fixed vector $\bm{t}\in\mathbb{R}^{d}$ and symmetric invertible matrix
$\bm{V}\in\mathbb{R}^{d\times d}$. Then, under the assumptions of Theorem 2,
the moment $\mathbb{E}\left\|\bm{X}_{i}\right\|^{3}$ is finite, and, in
addition, there exists a constant $C_{d}>0$ only depending on $d$ such that

$\sup_{\bm{x}\in\mathbb{R}^{d}}\Bigg{|}\Phi_{d}(\bm{x})-\mathbb{P}\left[N^{-1/2}\sum_{i=1}^{N}\bm{X}_{i}\leq\bm{x}\right]\Bigg{|}<C_{d}\mathbb{E}\left\|\bm{X}_{i}\right\|^{3}\big{/}\sqrt{N}+o\left(N^{-1/2}\right).$

Before I turn to its proof, note that this lemma implies that

$\sup_{\bm{x}\in\mathbb{R}^{d}}\Bigg{|}\Phi_{d}(\bm{x})-\mathbb{P}\left[N^{-1/2}\sum_{i=1}^{N}\bm{X}_{i}\leq\bm{x}\right]\Bigg{|}=o\left(N^{-1/2}\right)+o\left(N^{-1/2}\right).$

Proof. For the first assertion, note that
$J_{i,\bm{z}}=\mathds{1}\left\\{\bm{Z}_{i}=\bm{z}\right\\}$ follows a
Bernoulli distribution with probability parameter $f_{\varepsilon}(\bm{z})$
for any $\bm{z}\in\bm{\mathcal{Z}}$. An easy calculation reveals that

$\mathbb{E}\big{|}J_{i,\bm{z}}-f_{\varepsilon}(\bm{z})\big{|}^{3}=f_{\varepsilon}(\bm{z})(1-f_{\varepsilon}(\bm{z}))\big{(}(1-f_{\varepsilon}(\bm{z}))^{2}+f_{\varepsilon}(\bm{z})^{2}\big{)}.$
(B.2)

Furthermore, I have that

$\displaystyle\mathbb{E}\left\|\bm{W}(\bm{\theta}_{N}(\bm{t}))\big{(}\bm{J}_{i}-\bm{f}_{\varepsilon}\big{)}\right\|^{3}$
$\displaystyle=\mathbb{E}\left\|\sum_{\bm{z}\in\bm{\mathcal{Z}}}\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\mathds{1}\left\\{\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}\in[0,c]\right\\}\big{(}J_{i,\bm{z}}-f_{\varepsilon}(\bm{z})\big{)}\right\|^{3}$
$\displaystyle\leq\mathbb{E}\left\|\sum_{\bm{z}\in\bm{z}}\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\big{(}J_{i,\bm{z}}-f_{\varepsilon}(\bm{z})\big{)}\right\|^{3}$
$\displaystyle(\textit{Minkowski's inequality})\quad$
$\displaystyle\leq\Bigg{[}\sum_{\bm{z}\in\bm{\mathcal{Z}}}\Big{(}\|\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\|^{3}\mathbb{E}\big{|}J_{i,\bm{z}}-f_{\varepsilon}(\bm{z})\big{|}^{3}\Big{)}^{1/3}\Bigg{]}^{3}$
$\displaystyle(\textit{equation \eqref{eq:momentbernoulli}})\quad$
$\displaystyle=\Bigg{[}\sum_{\bm{z}\in\bm{\mathcal{Z}}}\Big{(}\|\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\|^{3}f_{\varepsilon}(\bm{z})(1-f_{\varepsilon}(\bm{z}))\big{(}(1-f_{\varepsilon}(\bm{z}))^{2}+f_{\varepsilon}(\bm{z})^{2}\big{)}\Big{)}^{1/3}\Bigg{]}^{3}$
$\displaystyle=:\widetilde{\mu}^{3},$

where the first line follows by definition of
$\bm{\theta}\mapsto\bm{W}(\bm{\theta})$. By Assumption A.5,
$\|\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\|^{3}<\infty$, and,
clearly,
$f_{\varepsilon}(\bm{z})(1-f_{\varepsilon}(\bm{z}))\big{(}(1-f_{\varepsilon}(\bm{z}))^{2}+f_{\varepsilon}(\bm{z})^{2}\big{)}<\infty$.
Subsequently, $\widetilde{\mu}^{3}$ is finite. It follows from the definition
of $\bm{X}_{i}$ that

$\begin{split}\mathbb{E}\left\|\bm{X}_{i}\right\|^{3}&\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{3}\mathbb{E}\left\|\bm{W}(\bm{\theta}_{N}(\bm{t}))\big{(}\bm{J}_{i}-\bm{f}_{\varepsilon}\big{)}\right\|^{3}\\\
&\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{3}\widetilde{\mu}^{3}\\\
&<\infty,\end{split}$ (B.3)

where the second line follows from the previous display, and the third line
from $\widetilde{\mu}^{3}<\infty$ and the fact that
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}<\infty$,
which is implied by Assumption A.5. This proves the first assertion.

For the second assertion, note that $\bm{J}_{i}$ obeys a multinomial
distribution with 1 trial, $m$ events, and event probabilities gathered in
vector $\bm{f}_{\varepsilon}$. Hence,
$\mathbb{E}\bm{J}_{i}=\bm{f}_{\varepsilon}$ and
$\mathbb{V}\mathrm{ar}\left[\bm{J}_{i}\right]=\mathrm{diag}(\bm{f}_{\varepsilon})-\bm{f}_{\varepsilon}\bm{f}_{\varepsilon}^{\top}=\bm{\Omega}$.
Therefore, $\mathbb{E}\bm{X}_{i}=\bm{0}$, as well as

$\displaystyle\mathbb{V}\mathrm{ar}\left[\bm{X}_{i}\right]$
$\displaystyle=\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\bm{W}(\bm{\theta}_{N}(\bm{t}))\bm{\Lambda}\bm{W}(\bm{\theta}_{N}(\bm{t}))^{\top}\left(\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\right)^{\top}$
$\displaystyle=\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\bm{U}\left(\bm{\theta}_{N}(\bm{t})\right)\left(\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\right)^{\top}$
$\displaystyle=\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\bm{U}\left(\bm{\theta}_{0}\right)\left(\bm{U}\left(\bm{\theta}_{0}\right)^{-1/2}\right)^{\top}+o\left(N^{-1/2}\right)$
$\displaystyle=\bm{I}_{d}+o\left(N^{-1/2}\right),$

where $\bm{I}_{d}$ denotes the the $d\times d$ identity matrix. The third line
follows from the fact that
$\bm{\theta}_{N}(\bm{t})=\bm{\theta}_{0}+N^{-1/2}\bm{V}\bm{t}=\bm{\theta}_{0}+o\left(N^{-1/2}\right)$
and the continuity of $\bm{\theta}\mapsto\bm{U}(\bm{\theta})$, where the
latter is implied by Assumption A.4. It follows that $\bm{X}_{i}$ has mean
zero and identity covariance matrix when $N$ is large enough. This fact
together with equation (B.3) allows me to apply Theorem 1.3 in Götze, (1991,
equation 1.5), and the second assertion follows. This completes the proof.
$\blacksquare$

######  Lemma B.3.

Grant the assumptions of Theorem 2 and put
$\bm{\theta}_{N}(\bm{t})=\bm{\theta}_{0}+N^{-1/2}\bm{V}\bm{t}$ for a fixed
vector $\bm{t}\in\mathbb{R}^{d}$ and symmetric invertible matrix
$\bm{V}\in\mathbb{R}^{d\times d}$. Then the equality

$\displaystyle-\sqrt{N}\sum_{\bm{z}\in\bm{\mathcal{Z}}}f_{\varepsilon}(\bm{z})\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\Bigg{(}{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}\Bigg{)}-{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}\Bigg{)}$
$\displaystyle\qquad=\sum_{\bm{z}\in\bm{\mathcal{Z}}}{w^{\prime}}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}f_{\varepsilon}(\bm{z})\bm{s}_{\bm{z}}\left(\bm{\theta}_{0}\right)\bm{s}_{\bm{z}}\left(\bm{\theta}_{0}\right)^{\top}\bm{V}\bm{t}+o\left(N^{-1/2}\right)+o\left(N^{-1/2}\right)+o\left(N^{-1/2}\right)$

holds true.

Proof. I can write

$\displaystyle-\sqrt{N}\sum_{\bm{z}\in\bm{\mathcal{Z}}}f_{\varepsilon}(\bm{z})\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)\Bigg{(}{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}\Bigg{)}-{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}\Bigg{)}$
$\displaystyle\quad=-\sqrt{N}\sum_{\bm{z}\in\bm{\mathcal{Z}}}\underbrace{\frac{{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}\Bigg{)}-{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}}{\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}-\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}}}_{=Ⓐ}\underbrace{\left(\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}-\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\right)}_{=Ⓑ}\underbrace{f_{\varepsilon}(\bm{z})\bm{s}_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}_{=Ⓒ}$

For term $Ⓐ$, put $h_{N}=N^{-1/2}\bm{V}\bm{t}$ such that
$\bm{\theta}_{N}(\bm{t})=\bm{\theta}_{0}+h_{N}$ and note that by the
definition of the derivative, it holds true that

$\lim_{h_{N}\to
0}\frac{{w}\Bigg{(}\frac{f_{\varepsilon}(k)}{p_{\bm{z}}\left(\bm{\theta}_{0}+h_{N}\right)}\Bigg{)}-{w}\Bigg{(}\frac{f_{\varepsilon}(k)}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}}{\frac{f_{\varepsilon}(k)}{p_{\bm{z}}\left(\bm{\theta}_{0}+h_{N}\right)}-\frac{f_{\varepsilon}(k)}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}}={w^{\prime}}\Bigg{(}\frac{f_{\varepsilon}(k)}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)},$

where the right hand side exists by Assumption A.9. By definition,
$h_{N}=o\left(N^{-1/2}\right)$. Combined with the previous display, one
obtains

$Ⓐ=\frac{{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}\Bigg{)}-{w}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}}{\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{N}(\bm{t})\right)}-\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}}={w^{\prime}}\Bigg{(}\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}_{0}\right)}\Bigg{)}+o\left(N^{-1/2}\right).$

For term $Ⓑ$, put
$g(\bm{\theta})=\frac{f_{\varepsilon}(\bm{z})}{p_{\bm{z}}\left(\bm{\theta}\right)}$
and note that