2.1 General framework

Several model specifications have been proposed for different data structures and designs(Gwet, Reference Gwet2008; Shrout & Fleiss, Reference Shrout and Fleiss1979; Ten Hove et al., Reference Ten Hove, Jorgensen and van der Ark2022). One-way designs are preferred when rater differences are typically considered as noise (Martinková et al., Reference Martinková, Bartoš and Brabec2023), whereas two-way designs are usually involved if the rater’s effect needs to be identified (Casabianca et al., Reference Casabianca, Lockwood and Mccaffrey2015; Mignemi et al., Reference Mignemi, Calcagnì, Spoto and Manolopoulou2024). Balanced designs require each subject to be rated by all the raters, while in an unbalanced design each subject is only rated by a generally small subset of them (Ten Hove et al., Reference Ten Hove, Jorgensen and van der Ark2021). Raters might be considered either fixed or random (i.e., drawn from the population) depending on the inference the researcher might be interested in (Koo & Li, Reference Koo and Li2016).

The unbalanced two-way design with random raters is considered a general case to present our model. The reasons for this choice are both theoretical and practical. We aim to provide a comprehensive statistical framework for modeling the dependency of ratings on different categorical predictors (i.e., subjects’ and raters’ identities). This setting is a neat compromise between the one-way design, which implies only one categorical predictor (i.e., subject identity), and more complex dependency structures that involve more than two identities (i.e., several categorical predictors). However, an example of these extensions is given in Appendix A.2. Indeed, our proposal might be alternatively reduced or extended to be suitable for these different levels of complexity. The unbalanced design implies some sparsity in the co-occurrence between subjects and raters and each subject is rated only by a small subset of raters (Papaspiliopoulos et al., Reference Papaspiliopoulos, Roberts and Zanella2019; Papaspiliopoulos et al., Reference Papaspiliopoulos, Stumpf-Fétizon and Zanella2023), as a consequence each rater might score a different number of subjects. This makes the framework general and flexible, it might be seen as an extension of cross-classified models in which uncertainty is modeled also hierarchically. From a practical perspective, our choice is reasonable since many large studies and applications use unbalanced designs to distribute the workload across different raters (Ten Hove et al., Reference Ten Hove, Jorgensen and van der Ark2022).

2.2 Preliminaries on Bayesian nonparametric inference

In this subsection, we briefly review some basic preliminaries on BNP inference providing here a very general framework which is detailed in Sections below (refer to Ghosal & van der Vaart, Reference Ghosal and van der Vaart2017 and Hjort et al., Reference Hjort, Holmes, Müller and Walker2010 for exhaustive treatments).

Suppose $Y_1,\dots , Y_n$ , are observations (e.g., ratings), with each $Y_i$ taking values in a complete and separable metric space $\mathbb {Y}$ . Let $\Pi $ denote a prior probability distribution on the set of all probability measures $\mathbf {P}_{\mathbb {Y}}$ such that

(1) $$ \begin{align} Y_i | p \overset{\mathrm{iid}}{\sim} p, \quad \;\; p \sim \Pi, \end{align} $$

for $i=1,\dots ,n$ . Here p is a random probability measure on $\mathbb {Y}$ and $\Pi $ is its probability distribution and might be interpreted as the prior distribution for Bayesian inference (De Blasi et al., Reference De Blasi, Favaro, Lijoi, Mena, Prünster and Ruggiero2015). The inferential problem is called parametric when $\Pi $ degenerates on a finite-dimensional subspace of $\mathbf {P}_{\mathbb {Y}}$ , and nonparametric when the support of $\Pi $ is infinite-dimensional (Hjort et al., Reference Hjort, Holmes, Müller and Walker2010, chapter 3). To the best of our knowledge, the vast majority of the contributions present in rating models literature (Bartoš & Martinková, Reference Bartoš and Martinková2024; Casabianca et al., Reference Casabianca, Lockwood and Mccaffrey2015; Martinková et al., Reference Martinková, Bartoš and Brabec2023; Martinkova et al., Reference Martinkova, Goldhaber and Erosheva2018; Nelson & Edwards, Reference Nelson and Edwards2010, Reference Nelson and Edwards2015; Ten Hove et al., Reference Ten Hove, Jorgensen and van der Ark2021, Reference Ten Hove, Jorgensen and van der Ark2022; Zupanc & Štrumbelj, Reference Zupanc and Štrumbelj2018) are developed within a parametric framework making use of a prior that assigns probability one to a small subset of $\mathbf {P}_{\mathbb {Y}}$ . Although Mignemi et al. (Reference Mignemi, Calcagnì, Spoto and Manolopoulou2024) recently proposed a Bayesian semi-parametric (BSP) model for analyzing rating data. Even if they relax the normality assumption for the rater effect (i.e., the systematic bias), normality is still assumed for the subject true score distribution. This strong prior assumption is overcome through a BNP approach (Ghosal & van der Vaart, Reference Ghosal and van der Vaart2017) in the present work.

Dirichlet processes. For the present proposal, we assume $\Pi $ to be a discrete nonparametric prior and correspond to a DP which has been widely used in BNP psychometric research (Cremaschi et al., Reference Cremaschi, De Iorio, Seng Chong, Broekman, Meaney and Kee2021; Karabatsos & Walker, Reference Karabatsos and Walker2009; Paganin et al., Reference Paganin, Paciorek, Wehrhahn, Rodríguez, Rabe-Hesketh and de Valpine2023; Yang & Dunson, Reference Yang and Dunson2010). Given $\Pi = DP(\alpha P_0)$ , p is a random measure on $\mathbb {Y}$ following a DP with concentration parameter $\alpha>0$ and base measure $P_0$ . This implies that for every finite measurable partition $\{B_1,\dots ,B_k\}$ of $\mathbb {Y}$ , the joint distribution $(p(B_1),\dots ,p(B_k))$ follows a k-variate Dirichlet distribution with parameters $\alpha P_0(B_1),\dots ,\alpha P_0(B_k)$ :

(2) $$ \begin{align} (p(B_1),\dots,p(B_k)) \sim Dir(\alpha P_0(B_1),\dots,\alpha P_0(B_k)). \end{align} $$

The base measure $P_0$ is our prior guess at p as it is the prior expectation of the DP, i.e., $\mathbf {E}[p]=P_0$ . The parameter $\alpha $ (also termed precision parameter) controls the concentration of the prior for p about $P_0$ . In the limit of $\alpha \rightarrow \infty $ , the probability mass is spread out and p gets closer to $P_0$ ; on the contrary, as $\alpha \rightarrow 0$ , p is less close to $P_0$ and concentrates at a point mass.

Dirichlet process mixtures. Given the discrete nature of the DP, whenever $\mathbb {Y} = \mathbb {R}$ it is not a reasonable prior for the real-valued random variable Y. Nonetheless, it might be involved in density estimation through hierarchical mixture modeling (Ghosal & van der Vaart, Reference Ghosal and van der Vaart2017). Let $f(\cdot ;\tilde {\theta })$ denote a probability density function for $\tilde {\theta } \in \Theta \subseteq \mathbb {R}$ , we modify (1) such that for $i=1,\dots ,n$ :

(3) $$ \begin{align} Y_i | \tilde{\theta}_i \overset{\mathrm{ind}}{\sim} f(\cdot;\tilde{\theta}_i), \quad \tilde{\theta}_i|p \overset{\mathrm{iid}}{\sim} p, \quad \;\; p \sim DP(\alpha P_0). \end{align} $$

The realizations of the DP are almost surely (a.s.) discrete which implies a positive probability that $\tilde {\theta }_i = \tilde {\theta }_{i'}$ , for $i \neq i'$ . Indeed, a random sample $(\tilde {\theta }_1,\dots , \tilde {\theta }_n)$ from p features $1\leq K_n \leq n$ different unique values $(\tilde {\theta }^{\ast}_1,\dots , \tilde {\theta }^{\ast}_{K_n})$ and leads to a random partition of $\{1,\dots , n\}$ into $K_n$ blocks such that $ \tilde {\theta }_i \in (\tilde {\theta }^{\ast}_1,\dots , \tilde {\theta }^{\ast}_{K_n})$ for $i=1,\dots ,n$ . This naturally induces a mixture distribution for the observations $Y_1,\dots , Y_n$ with probability density:

(4) $$ \begin{align} f(Y) = \int f(Y;\tilde{\theta})p(d\tilde{\theta}). \end{align} $$

To provide some intuition, by using a DP as a prior for an unknown mixture distribution we mix parametric families nonparametrically (Gelman et al., Reference Gelman, Hwang and Vehtari2014). This model specification introduced by Lo (Reference Lo1984) and termed DPM provides a BNP framework to model rating data.

2.3 Proposed model

Consider a subject $i = 1, \dots , I$ , whose attribute is independently scored by a finite random subset of raters $\mathcal {R}_i \subseteq \{1,\dots , J\}$ on a continuous rating scale. We assume that the observed rating $Y_{ij} \in \mathbb {R}$ depends independently on subject i and rater $j \in \mathcal {R}_i$ . The effect of the former is interpreted as i’s true score and is the rating procedure’s focus. We let the residual part, that is the difference between the true and the observed score, depend on rater j’s effects, i.e., systematic bias and reliability.

Modelling rating $Y_{ij}$ . We specify the following decomposition of rating $Y_{ij}$ :

(5) $$ \begin{align} Y_{ij} = \theta_i + \tau_j +\epsilon_{ij}, \quad i = 1,\dots,I; \;\; j \in \mathcal{R}_i. \end{align} $$

Here $\theta _i$ captures the subject i’s latent “true” score and $ \tau _j +\epsilon _{ij}$ is the difference between the observed and the true score, representing the error of rater j. The true score $\theta _i$ might be regarded as the latent consensus across the raters $\mathcal {R}_i$ on the attribute of subject i.We assume these terms to be mutually independent.

Modelling subject’s true score. For each subject $i = 1,\dots ,I$ we assume that the true score $\theta _i$ is independently distributed following a normal distribution with mean $\mu _i$ and variance $\omega ^2_i$ :

(6) $$ \begin{align} \theta_i|\mu_i,\omega^2_i \ \overset{\mathrm{ind}}{\sim} \ N(\mu_i,\omega^2_i). \end{align} $$

Here $\mu _i$ is the mean of subject i’s true score, $\omega ^2_i$ is its variability and we assume them to be independent. Conditional on the rater’s error, higher values of $\theta _i$ imply higher levels of the subjects’ attribute (e.g., higher student proficiency); on the contrary lower values indicate poor levels of their attribute (e.g., poor student proficiency).

We specify a DP prior with precision parameter $\alpha _1$ and base measure $G_0$ for the pair $(\mu _i,\omega ^2_i)$ , $i=1,\dots ,I$ :

(7) $$ \begin{align} (\mu_i,1/\omega^2_i)|G \overset{\mathrm{iid}}{\sim} G, \quad \; G \sim DP(\alpha_1G_0). \end{align} $$

We choose $G_0=N(\mu _0, S_0) \times Ga(w_0, w_0/W_0)$ , where $\mu _0$ and $S_0$ are the mean and variance of the normal distribution and $w_0$ and $W_0$ are, respectively, the shape and the mean parameters of the gamma. We note that G is a.s. discrete with a non-zero probability of ties, such that different subjects will share the same values of $(\mu _i,1/\omega ^2_i)$ with a probability greater than zero, that is $P[(\mu _i,1/\omega ^2_i)=(\mu _{i'},1/\omega ^2_{i'})]>0$ , for $i \neq i'$ . This discreteness property naturally induces clustering across subjects and leads to a location-scale DPM prior for $\theta _i$ . That is, this formulation can capture clusters of subject abilities. Figure 1 shows the hierarchical dependence of subjects’ true scores.

Figure 1 Graphical representation of the dependencies implied by the model. The boxes indicate replicates, the four outer plates represent, respectively, subjects and raters, and the inner grey plate indicates the observed rating.

Modelling rater’s bias and reliability. For each rater $j=1,\dots , J$ , who scores a subset of subjects $\mathcal {S}_j \subseteq \{1,\dots ,I\}: j \in \mathcal {R}_i$ , the difference between the observed rating $Y_{ij}$ and the subject’s true score $\theta _i$ , $i \in \mathcal {S}_j$ , is decomposed into the rater effects $\tau _j$ and $\epsilon _{ij}$ (5), assuming $\tau _j {\perp \!\!\! \perp } \epsilon _{ij}$ . We model $\tau _j$ to be normally distributed with mean $\eta _j$ and variance $\omega ^2_j$ :

(8) $$ \begin{align} \tau_j|\eta_j,\phi^2_j \overset{\mathrm{ind}}{\sim} N(\eta_j,\phi^2_j), \quad j=1,\dots,J. \end{align} $$

Here $\eta _j$ and $\phi ^2_j$ are the mean and the variance of the rater j’s effect $\tau _j$ . It captures j’s specific systematic bias, i.e., the mean difference between the observed rating $Y_{ij}$ and the subject’s true score $\theta _i$ , $i \in \mathcal {S}_j$ . Given two raters such that $\tau _j<\tau _{j'}$ , j is said to be more strict and expected to give systematically smaller ratings than j on average.

The residual term $\epsilon _{ij}$ is assumed to be i.i.d. for $i \in \mathcal {S}_j$ following a normal distribution with zero mean and variance $\sigma ^2_j$ . We let this parameter vary across raters and assume $1/\sigma ^2_j$ follows a gamma distribution with shape and rate parameters $\gamma _j,\gamma _j/\beta _j$ , respectively:

(9) $$ \begin{align} \epsilon_{ij}|\sigma^2_j & \ \ \overset{\mathrm{iid}}{\sim}\ \ N(0, \sigma^2_j), \quad i \in \mathcal{S}_j, \end{align} $$

(10) $$ \begin{align} 1/\sigma^2_j|\gamma_j,\beta_j & \ \ \overset{\mathrm{ind}}{\sim} \ \ Ga(\gamma_j,\gamma_j/\beta_j), \quad j=1,\dots,J. \end{align} $$

Under this parametrization, $1/\sigma ^2_j$ is the rater j’s specific reliability with mean $\beta _j$ and $\gamma _j$ is the shape parameter. We prefer this parametrization for interpretability purposes, which implies a simpler notation below. Conditional on subjects’ true score $\theta _i$ , $i \in \mathcal {S}_j$ , larger values of $\sigma ^2_j$ imply more variability across the ratings given by j and might be interpreted as a poorly consistent rating behaviour. On the contrary, smaller values of $\sigma ^2_j$ indicate less variability and higher consistency for j across subjects. As a consequence, the parameter $1/\sigma ^2_j$ might be equivalently referred to as the rater-specific precision.

We specify a DP prior with concentration parameter $\alpha _2$ and base measure $H_0$ for the four-dimensional vector $(\eta _j,1/\phi ^2_j,\gamma _j,1/\beta _j)$ , $j=1,\dots ,J$ :

(11) $$ \begin{align} (\eta_j,1/\phi^2_j,\gamma_j,1/\beta_j)|H \overset{\mathrm{iid}}{\sim} H, \quad \; H \sim DP(\alpha_2H_0). \end{align} $$

We assume mutual independence for the elements of the vector and choose $H_0=N(\eta _0, D_0) \times Ga(a_0, a_0/A_0) \times Ga(b_0, b_0/B_0) \times Ga(m_0, m_0/M_0)$ , where $\eta _0$ and $D_0$ are mean and scale parameters, respectively; $a_0, b_0, m_0$ are shape parameters and $A_0, B_0, M_0$ are mean parameters. This formulation induces a DPM prior for raters’ bias and reliability $\tau _j$ and $1/\sigma ^2_j$ . Figure 1 gives a graphical representation of the model. The independence assumption might be relaxed by employing a suitable multivariate base measure accounting for possible dependencies among the four elements of the vector. However, this implies a more complex specification, which is beyond the purpose of this work. Further constraints on raters’ systematic bias $\tau $ are needed for identifiability purposes which are discussed in Section 2.5, after presenting the stick-breaking representation.

2.4 Stick-breaking representation

The random probability measures G and H are assigned discrete priors, as a consequence they might be represented as a weighted sum of point masses:

(12) $$ \begin{align} G &= \sum_{n \geq1} \pi_{1n} \delta_{\xi_n}, \end{align} $$

(13) $$ \begin{align} H &= \sum_{k\geq1} \pi_{2k} \delta_{\zeta_k}, \end{align} $$

where the weights $\{\pi _{1n}\}_{n=1}^{\infty }$ and $\{\pi _{2k}\}_{k=1}^{\infty }$ take values on the infinite probability simplex and $\delta _{x}(\cdot )$ stands for the Dirac measure and denotes a point mass at x. Note that, we index the components of the infinite mixture (12) corresponding to the subjects with $n=1,\dots ,\infty $ , whereas $k=1,\dots ,\infty $ is used for that corresponding to the raters (13). The random vectors $\xi _n=(\mu _n,\omega ^2_n)$ , $n=1,\dots ,\infty $ are i.i.d. from the base measure $G_0$ , $\zeta _k=(\eta _k,\phi ^2_k,\gamma _k,\beta _k)$ , $k=1,\dots ,\infty $ are i.i.d. from the base measure $H_0$ , and both vectors are assumed to be independent of the corresponding weights. This makes clear why the expectations of the $DPs$ are $G_0$ and $H_0$ , respectively, and are said to be our prior guess at G and H (see Section 2.2).

This discreteness property of the $DP$ allows us to define G and H through the stick-breaking representation introduced by Sethuraman (Reference Sethuraman1994):

(14) $$ \begin{align} G = \sum_{n \geq1} \pi_{1n} \delta_{\xi_n}, \quad \pi_{1n} = V_{1n} \prod_{l<n}(1-V_{1l}), \quad V_{1n} \overset{\mathrm{iid}}{\sim} Beta(1,\alpha_1), \quad \xi_n \overset{\mathrm{iid}}{\sim} G_0, \end{align} $$

and

(15) $$ \begin{align} H = \sum_{k \geq1} \pi_{2k} \delta_{\zeta_k}, \quad \pi_{2k} = V_{2k} \prod_{l<k}(1-V_{2l}), \quad V_{2k} \overset{\mathrm{iid}}{\sim} Beta(1,\alpha_2), \quad \zeta_k \overset{\mathrm{iid}}{\sim} H_0. \end{align} $$

This construction of the $DP$ implies that, for each subject $i=1,\dots , I$ , $(\mu _i, \omega ^2_i)=\xi _n$ with probability $ \pi _{1n} = V_{1n} \prod _{l<n}(1-V_{1l})$ . Equivalently, for each rater $j=1,\dots ,J$ , the probability that $(\eta _j,\phi ^2_j,\gamma _j,\beta _j)=\zeta _k$ is given by $\pi _{2k} = V_{2k} \prod _{l<k}(1-V_{2l})$ .

Moments of student latent true score  $\theta _i$ . The mean and the variance of the subject’s true score $\theta _i$ , $i=1,\dots , I$ , under a $DP(\alpha _1G_0)$ prior are:

(16) $$ \begin{align} \mathbf{E}[\theta_i|G] = \mu_{G} =\sum_{n\geq 1} \pi_{1n} \mu_n, \quad \quad \mathbf{Var}[\theta_i|G] = \omega^2_{G} =\sum_{n\geq 1} \pi_{1n} (\mu_n^2 + \omega_n^2) - \mu^2_{G}, \end{align} $$

where $\mu _n$ and $\omega _n^2$ are the mean and the variance of $\theta _i$ for the nth component of the mixture. Here $\mu _{G}$ is the weighted average across components and captures the mean true score across subjects. The parameter $ \omega ^2_{G}$ is the conditional variance of the infinite mixture and indicates the variability of true scores across subjects.

Moments of raters’ bias $\tau _j$ . The mean and the variance of the rater’s bias $\tau _j$ , $j=1,\dots , J$ , under a $DP(\alpha _2H_0)$ prior are:

(17) $$ \begin{align} \mathbf{E}[\tau_j|H] = \eta_{H} =\sum_{k \geq 1} \pi_{2k} \eta_k, \quad \quad \mathbf{Var}[\tau_j|H] = \phi^2_{H} =\sum_{k \geq 1} \pi_{2k} (\eta_k^2 + \phi_k^2) - \eta^2_{H}, \end{align} $$

where $\eta _k$ and $\phi _k^2$ are the mean and the variance of $\tau _j$ for the kth component of the mixture. Here $\eta _{H}$ and $\phi ^2_{H}$ capture the mean and the variance of the systematic bias within the general population of raters.

Moments of raters’ reliability $1/\sigma ^2_j$ . Raters’ residual mean is fixed to zero by the model (9), that is $\mathbf {E}[\epsilon ]=0$ ; mean and variance of raters reliability $1/\sigma ^2_j$ under a $DP(\alpha _2H_0)$ prior are:

(18) $$ \begin{align} \mathbf{E}[1/\sigma^2_j|H] = \beta_{H} =\sum_{k \geq 1} \pi_{2k} \beta_k \quad \quad \mathbf{Var}[1/\sigma^2_j|H] = \psi^2_{H} =\sum_{k \geq 1} \pi_{2k} (\beta_k^2 + \psi_k) - \beta^2_{H}, \end{align} $$

where $\beta _{H}$ captures raters’ weighted average reliability and $\psi ^2_{H}$ indicates the total reliability variance across them. Here $\beta _k$ and $\psi _k=\beta _k^2/\gamma _k$ are, respectively, the mean and the variance of $1/\sigma ^2_j$ for the kth component of the mixture.

Note that we model the independent rater’s features, i.e., bias and reliability, by placing the same $DP(\alpha _2H_0)$ prior. In other terms, $\tau _j$ and $1/\sigma _j$ are two independent elements of the same vector drawn from H.

Finite stick-breaking approximation. The recursive generation defined in (14) and (15) implies a decreasing stochastic order of the weights $\{\pi _{1n}\}_{n=1}^{\infty }$ and $\{\pi _{2k}\}_{k=1}^{\infty }$ as the indices n and k grow. Considering the expectations $\mathbf {E}[V_{1n}]=1/(1+\alpha _1)$ and $\mathbf {E}[V_{2k}]=1/(1+\alpha _2)$ it is clear that the rates of decreasing depend on the concentration parameters $\alpha _1$ and $\alpha _2$ , respectively. Values of these parameters close to zero imply a mass concentration on the first couple of atoms, with the remaining atoms being assigned small probabilities; which is consistent with the general formulation of the $DP$ discussed in Section 2.2. Given this property of the weights, in practical applications the infinite sequences (12) and (13), are truncated at enough large values of $R \in \mathbb {N}$ :

(19) $$ \begin{align} G = \sum_{n =1}^{R} \pi_{1n} \delta_{\xi_n}, \quad \quad H = \sum_{k =1}^{R} \pi_{2k} \delta_{\zeta_k}. \end{align} $$

We use this finite stick-breaking approximation proposed by Ishwaran & James (Reference Ishwaran and James2001) to let $V_{1R}=V_{2R}=1$ , and discard the terms $R+1,\dots ,\infty $ , for G and H.

The moment formulas (16), (17), and (18) are readily modified accordingly to the truncation and computed as finite mixture moments.

Nested versions. Semiparametric nested versions of the BNP model might be specified in which alternatively G or H are degenerate on a single component and $R=1$ for one of them in the finite approximation. That is, subjects or raters are all clustered together. For instance, for very small values of J (i.e., raters’ sample size), raters might not be reasonably considered a representative sample of their population and limited information is available for drawing inference about it. Under these scenarios, raters’ effects might be assumed to be i.i.d. from a normal distribution.

2.5 Semi-centered DPM

Hierarchical models (e.g., GLMM, Linear Latent Factor models), might suffer from identifiability issues, and constraints on the latent variable distributions are needed for consistently identify and interpret model parameters (Bartholomew et al., Reference Bartholomew, Knott and Moustaki2011; Gelman & Hill, Reference Gelman and Hill2006; Yang & Dunson, Reference Yang and Dunson2010). More specifically, under the linear random effects models a standard procedure to achieve model identifiability is to constrain the mean of the random effects to be zero (Agresti, Reference Agresti2015). We aim to consistently involve the same mean constraint for our proposal and allow straightforward and interpretable comparisons between the parametric and the nonparametric models. Similar to Yang & Dunson (Reference Yang and Dunson2010), we encompass a DPM-centered prior such that the expected value of the rater systematic bias is fixed to zero, $\mathbf {E}[\tau _j]=0$ , for $j=1,\dots ,J$ .

Since the rating process focuses on the subjects’ true scores, it might be more reasonable to centre the DPM for the raters’ effects and let the model estimate the mean of the true scores $\mu _G$ . Given that the mean of the raters’ residual is fixed to zero in (9), the mean raters’ bias needs to be fixed. We adapt the centering procedure based on a parameter-expanded approach proposed by Yang et al., Reference Yang, Dunson and Baird2010 and Yang & Dunson, Reference Yang and Dunson2010 to our proposal. We specify a semi-centered DPM (SC-DPM) involving an expansion in raters’ systematic bias $\{\tau ^{\ast}_j\}_1^J$ , such that their mean $\eta ^{\ast}_H=0$ a.s. The expanded-parameters (8) can be expressed as

(20) $$ \begin{align} \tau^{\ast}_j = \tau_j - \eta_H, \quad \quad \tau_j|\eta_j,\phi^2_j \overset{\mathrm{ind}}{\sim} N(\eta_j,\phi^2_j), \quad \quad j=1,\dots, J, \end{align} $$

and the decomposition of rating $Y_{ij}$ (5) becomes

(21) $$ \begin{align} Y_{ij} = \theta_i + \tau^{\ast}_j +\epsilon_{ij}, \quad \quad i = 1,\dots,I; \;\; j \in \mathcal{R}_i. \end{align} $$

Given the location transformation in (20) the expectation of the expanded parameters is zero:

(22) $$ \begin{align} \mathbf{E}[\tau^{\ast}_j|H] = 0. \end{align} $$

It is worth noting that the centering needs only to concern the location of the systematic bias and not its scale as it is in the centered-DPM introduced by Yang & Dunson (Reference Yang and Dunson2010), which explains the term “semi-centring” adopted here to avoid confusion. Accordingly, under the semiparametric specifications, the only location of the parametric distribution needs to be fixed; a zero mean normal distribution might be a suitable solution.

4.1 Identifiability and ICC

The moments of the error $\epsilon _{ij}$ , $i=1\dots ,I$ and $ j \in \mathcal {R}_i$ , are, respectively:

(37) $$ \begin{align} \mathbf{E}[\epsilon_{ij}|H] = \eta_{H} =\sum_{k\geq1} \pi_{2k} \eta_k, \quad \quad \mathbf{Var}[\epsilon_{ij}|H] = \phi^2_{H} =\sum_{k\geq1} \pi_{2k} (\eta_k^2 + \phi_k^2) - \eta^2_{H}. \end{align} $$

The centering strategy detailed in Section 2.5 is here used and a SC-DPM is here placed over $\epsilon _{ij}$ :

(38) $$ \begin{align} \epsilon^{\ast}_{ij} = \epsilon_{ij} - \eta_H, \quad \quad \epsilon_{ij}|\eta_{ij},\phi^2_{ij} \overset{\mathrm{ind}}{\sim} N(\eta_{ij},\phi^2_{ij}), \quad \quad i = 1,\dots,I; \;\; j \in \mathcal{R}_i. \end{align} $$

Under this parameter-expanded specification, the decomposition of rating $Y_{ij}$ (34) becomes

(39) $$ \begin{align} Y_{ij} = \theta_i +\epsilon^{\ast}_{ij}, \quad \quad \quad i = 1,\dots,I; \;\; j \in \mathcal{R}_i. \end{align} $$

Given the location transformation in (38), the expectation of the residuals is zero:

(40) $$ \begin{align} \mathbf{E}[\epsilon^{\ast}_{ij}|H] = 0. \end{align} $$

For the one-way designs, the exact general ICC might be consistently estimated.

Proposition 3. Given a random subject i, $i=1,\dots ,I$ , independently scored by two random raters $j, j' \in \mathcal {R}_i$ , $j \neq j'$ , the conditional correlation between the ratings $Y_{ij}$ and $Y_{ij'}$ is

(41) $$ \begin{align} Corr\left(Y_{ij}, Y_{ij'} |G,H \right) = ICC = \frac{\omega^2_G}{\omega^2_G + \phi^2_H}. \end{align} $$

The proof is given in Appendix B. Conditioning on different clusters of subjects or raters and different ICC formulations lead to possible comparisons among clusters similar to the main model.

5.1 Hyperprior specification

Eliciting the concentrations ‘and base measures’ parameters has a role in controlling the posterior distribution over clustering (Gelman et al., Reference Gelman, Hwang and Vehtari2014). Small values of the variance parameters of the base measures $G_0$ , and $H_0$ favor the clustering of subjects and raters, respectively, to different clusters. On the contrary, larger values of $G_0$ and $H_0$ variances favor the allocation of different subjects and raters, respectively, to the same cluster.

We improve model flexibility by placing a prior on the base measures $G_0$ and $H_0$ , and the concentration parameters $\alpha _1$ and $\alpha _2$ letting them be informed by the data. For the subjects’ true score base measure $G_0=N(\mu _0, S_0) \times Ga(w_0, w_o/W_0)$ the following hyperpriors are specified:

$$ \begin{align*} \mu_0 \sim N(\lambda_{\mu_0},\kappa^2_{\mu_0}), \quad S_0 \sim IGa(q_{S_0},Q_{S_0}), \quad w_0 \sim Ga(q_{w_0},Q_{w_0}), \quad W_0 \sim IGa(q_{W_0},Q_{W_0}). \end{align*} $$

We let $\lambda _{\mu _0}$ be the rating scale’s center value (e.g., $\lambda _{\mu _0}=50$ on a 1-100 rating scale), $\kappa ^2_{\mu _0}=100$ and the parameters $q_{w_0}, Q_{w_0},q_{W_0}, Q_{W_0}$ equal to 0.005. For the raters’ base measure $H_0=N(\eta _0, D_0) \times Ga(a_0, a_0/A_0) \times Ga(b_0, b_0/B_0) \times Ga(m_0, m_0/M_0)$ , the following hyperpriors are specified:

$$ \begin{align*} \eta_0 \sim N(\lambda_{\eta_0},\kappa^2_{\eta_0}), \quad D_0 \sim IGa(q_{D_0},Q_{D_0}), \quad a_0 \sim Ga(q_{a_0},Q_{a_0}), \quad A_0 \sim IGa(q_{A_0},Q_{A_0}), \end{align*} $$

$$ \begin{align*} b_0 \sim Ga(q_{b_0},Q_{b_0}), \quad B_0 \sim IGa(q_{B_0},Q_{B_0}), \quad m_0 \sim Ga(q_{m_0},Q_{m_0}), \quad M_0 \sim IGa(q_{M_0},Q_{M_0}). \end{align*} $$

Where $\lambda _{\eta _0}=0$ , $\kappa ^2_{\eta _0}=100$ , and the other hyperparameters are fixed to 0.005. The concentration parameters $\alpha _1$ and $\alpha _2$ are assumed to follow, respectively, a gamma distribution:

$$ \begin{align*} \alpha_1 \sim Ga(a_1,A_1) \quad \quad \alpha_2 \sim Ga(a_2, A_2). \end{align*} $$

where $a_1,A_1,a_2,A_2$ are fixed to 1. The values we fix for the hyperprior’s parameters are very common in literature and they are consistent with those proposed by many other studies on BNP models (e.g., Gelman et al., Reference Gelman, Hwang and Vehtari2014; Heinzl et al., Reference Heinzl, Fahrmeir and Kneib2012; Mignemi et al., Reference Mignemi, Calcagnì, Spoto and Manolopoulou2024; Paganin et al., Reference Paganin, Paciorek, Wehrhahn, Rodríguez, Rabe-Hesketh and de Valpine2023; Yang & Dunson, Reference Yang and Dunson2010).

5.2 Posterior computation

Since most of the parameters in the model have conjugate prior distributions, a Blocked Gibbs sampling algorithm was used for the posterior sampling (Ishwaran & James, Reference Ishwaran and James2001). No conjugate priors are available for the gamma’s shape parameters (e.g., $\gamma _k$ , $k=1,\dots , R$ , $a_0$ , $b_0$ ), thus we approximate the full conditionals using a derivatives-matching procedure (D-M) which is involved as an additional sampling step within the MCMC. This method has several advantages over other sampling schemes (e.g., adaptive rejection sampling or Metropolis-Hasting) in terms of efficiency, flexibility, and convergence property (Miller, Reference Miller2019). We use the same D-M algorithm introduced by Miller, Reference Miller2019 to approximate the posterior of the gamma shape parameters of the base measures, i.e., $w_0, a_0, b_0, m_0$ and a modified version for the parameters $\gamma _k$ , $k=1,\dots , R$ , since the parametrization (30) is adopted. We detail this adapted version of the D-M algorithm in the paragraph below and provide the complete Gibbs sampling in the Supplementary Material.

The notation on the independent allocation of subjects and rater to the corresponding clusters is introduced here. Let $c_{1i}$ denote the cluster allocation of subject $i=1,\dots , I$ , with $c_{1i}=n$ whenever $\xi _i=\xi _n$ , $n=1,\dots , R$ . Given the finite stick-breaking approximation detailed in Section 2.4, R is the maximum number of clusters. We indicate the set of all the subjects assigned to the nth cluster with $\mathcal {C}_{1n}$ and with $N_{1n}=|\mathcal {C}_{1n}|$ its cardinality. Accordingly, let $c_{2j}$ denote the cluster allocation of rater $j=1,\dots , J$ , such that $c_{2j}=k$ whenever $\zeta ^{\ast}_j=\zeta ^{\ast}_k$ , $k=1,\dots , R$ . The set of all the raters assigned to the kth cluster is denoted by $\mathcal {C}_{2k}$ with $N_{2k}=|\mathcal {C}_{2k}|$ being its cardinality.

Derivatives-matching procedure. Since no conjugate priors are available for the gamma’s shape parameters $\{\gamma _k\}_1^R$ , we involve, for each of these parameters, a D-M procedure to find a gamma distribution that approximates the full conditional distribution of these parameters, when their prior is also a gamma distribution (Miller, Reference Miller2019).

We aim to approximate $p(\gamma _k|\cdot )$ , i.e., the true full conditional density of $\gamma _k$ , $k=1,\dots , R$ , by finding $U_{1k}$ and $U_{2k}$ such that

(42) $$ \begin{align} p(\gamma_k|\cdot) \approx g(\gamma_k|U_{1k},U_{2k}), \quad k=1,\dots,R, \end{align} $$

where $g(\cdot )$ is a gamma density, $U_{1k}$ and $U_{2k}$ are shape and rate parameters, respectively. The algorithm aims to find $U_{1k}$ and $U_{2k}$ such that the first and the second derivatives of the corresponding log densities of $p(\gamma _j|\cdot )$ and $g(\gamma _k|U_{1k}, U_{2k})$ match at a point $\gamma _{k}$ . Miller (Reference Miller2019) suggest to choose $\gamma _{k}$ to be near the mean of $p(\gamma _k|\cdot )$ for computational convenience. The approximation is iteratively refined by matching derivatives at the current $g(\cdot )$ mean as shown by Algorithm 1. We adapt the algorithm to our proposal, more specifically we consider the model involving the shape constraint introduced in equation (30). When this constraint is not imposed, the original algorithm by Miller (Reference Miller2019) may be directly used.

We denote with $X_{1k}$ and $X_{2k}$ the sufficient statistics for $\gamma _k$ corresponding to the kth raters’ mixture component. For the implementation of the Algorithm 1 we set the convergence tolerance $\epsilon _0=10^{-8}$ and the maximum number of iterations $M=10$ . Here $\psi (\cdot )$ and $\psi '(\cdot )$ are the digamma and trigamma functions, respectively.

The parameters $U_{1k}$ and $U_{2k}$ , returned by the algorithm, are used to update $\gamma _k \sim Ga(U_{1k},U_{2k})$ , $k=1,\dots ,R$ , through the MCMC sampling. The derivation of the algorithm is given in the Supplementary Material.

5.3 Post-processing procedures

Semi-centered DPM processes. The sampling scheme detailed in the Supplementary Material provides draws under the noncentered DPM model. However, as discussed in Section 2.5, it is not identifiable, and we need to post-process the MCMC samples to make inferences under the SC-DPM parameter-expanded model Yang & Dunson (Reference Yang and Dunson2010). Since it is a semi-centered model that naturally constrains the raters’ systematic bias $\{\tau _j\}_1^J$ to have zero mean, a few location transformations are needed. After computing $\eta _H$ according to 17 for each iteration, the samples of $\mu _0,\mu _G,\{\theta _i\}_1^I$ , and $\{\tau _j\}_1^J$ are computed:

$$ \begin{align*} \mu_0^{\ast} &= \mu_0 + \eta_H, \\ \mu_G^{\ast} & = \mu_G + \eta_H, \\ \theta_i^{\ast} & = \theta_i + \eta_H, \quad \quad \text{for} \;\; i=1,\dots,I; \\ \tau_j^{\ast} & = \tau_j - \eta_H, \quad \quad \text{for} \;\; j=1,\dots,J. \end{align*} $$

The first three are due to the location transformation of $\tau _j$ and have to be considered for inference purposes under the SC-DPM model.

Posterior densities and clusters point estimates. Each density equipped with a BNP prior might be monitored along the MCMC by a dense grid of equally spaced points (Gelman et al., Reference Gelman, Hwang and Vehtari2014; Mignemi et al., Reference Mignemi, Calcagnì, Spoto and Manolopoulou2024; Yang & Dunson, Reference Yang and Dunson2010). Each point of the grid is evaluated according to the mixture resulting from the finite stick-breaking approximation at each iteration. At the end of the MCMC, for each point of the grid posterior mean and credible interval might be computed, and as a by-product, the pointwise posterior distribution of the density might be represented.

The BNP model provides a posterior over the entire space of subjects’ and raters’ partitions, respectively. However, we can summarize these posteriors and determine the point estimates of these clustering structures by minimizing the respective variation of information (VI) loss functions. We refer to Wade & Ghahramani (Reference Wade and Ghahramani2018) and Meilă (Reference Meilă2007) for further details on VI and point estimates of probabilistic clustering.

As for every parameter of the model, we use the posterior distribution of the subjects’ specific parameters for inference purposes. Point estimates of the subjects’ true scores $\{\theta _i\}_1^I$ , such as the posterior mean (i.e., expected a posteriori, EAP) or the maximum a posteriori (i.e., MAP), might be used as official evaluations (i.e., final grades), and the posterior credible intervals as uncertainty quantification around those values. The $ICC_A$ index (32) can be computed at each iteration of the MCMC to get its posterior distribution, which might be used for inference purposes.

Computational details. In the present work, both for the simulations and the real data analysis, similarly to previous works (e.g., Heinzl et al., Reference Heinzl, Fahrmeir and Kneib2012; Paganin et al., Reference Paganin, Paciorek, Wehrhahn, Rodríguez, Rabe-Hesketh and de Valpine2023), the number of iterations is fixed to 80,000 (with a thin factor of 60 due to memory constraints), discarding the first 20,000 as burn-in. We fix the maximum number of clusters to be $R=25,$ respectively, for subjects’ and raters’ DPM priors (Gelman et al., Reference Gelman, Hwang and Vehtari2014). The package mcclust.ext (Wade, Reference Wade2015) is used for the point estimate of the clustering structures based on the VI loss functions. We graphically check out trace plots for convergence and use the package coda for model diagnostics (De Iorio et al., Reference De Iorio, Favaro, Guglielmi and Ye2023; Plummer et al., Reference Plummer, Best, Cowles and Vines2006). Convergence is also confirmed through multiple runs of the MCMC with different starting valuesFootnote ¹.

6.1 Setting

We generate subjects’ ratings on a continuous scale, $Y_{ij} \in (1,100)$ , the number of subjects $I=500$ and raters $J=100$ are fixed, whereas the number of ratings per subject and the true generative model vary across scenarios.

Generative scenarios. We manipulate the number of ratings per subject to be $|\mathcal {R}_i| \in \{2,4\}$ for $i=1,\dots ,I$ , since in many real contexts (e.g., education, peer review) it is common for the subjects to be rated only by two or few more independent raters (Zupanc & Štrumbelj, Reference Zupanc and Štrumbelj2018).

Data are generated as specified by equations 5, 9, and one of the schemes below, according to the three different scenarios:

• Unimodal: Under this scenario, subjects’ true score and raters’ effects densities are unimodal:

  $$ \begin{align*} \theta_i \overset{\mathrm{iid}}{\sim} N(50,50), \quad \quad (\tau_j, 1/\sigma_j^2) \overset{\mathrm{iid}}{\sim} N(0,25) \;Ga(10,10/0.15), \end{align*} $$
  for $i=1,\dots ,I$ and $j=1,\dots ,J$ . This corresponds to the standard BP model in which subjects’ true scores are assumed to be i.i.d across subjects and raters’ effects are drawn jointly i.i.d. across raters.

• Bimodal: In this scenario, both subjects’ and raters’ populations are composed, respectively, of two different clusters:

  $$ \begin{align*} \theta_i &\overset{\mathrm{iid}}{\sim} 0.7 \cdot N(39,50) + 0.3 \cdot N(75.6,30) , \\ (\tau_j, 1/\sigma_j^2) &\overset{\mathrm{iid}}{\sim} 0.5 \cdot N(-5,10) \;Ga(10,10/0.1) + 0.5 \cdot N(5,5) \;Ga(10,10/0.2), \end{align*} $$
  for $i=1,\dots ,I$ and $j=1,\dots ,J$ .

• Multimodal: Under this scenario, both subjects and raters are assigned, respectively, to three clusters:

  $$ \begin{align*} \theta_i &\overset{\mathrm{iid}}{\sim} 0.2 \cdot N(35,50) + 0.2 \cdot N(45,20) + 0.6 \cdot N(56.6,20) , \\ (\tau_j, 1/\sigma_j^2) & \overset{\mathrm{iid}}{\sim} 0.4 \cdot SN(-5,3.162,-5) \;Ga(10,10/0.15) \\ &+ 0.4 \cdot N(0,10) \;Ga(10,10/0.10) \\ &+ 0.2 \cdot N(10,10) \;Ga(10,10/0.20), \end{align*} $$
  for $i=1,\dots ,I$ and $j=1,\dots ,J$ . Here $SN(\xi ,\omega ,\alpha )$ stands for the skew-normal distribution with location, scale and slant parameters, $\xi $ , $\omega $ and $\alpha $ , respectively.

These scenarios mimic three different levels of heterogeneity. From an interpretative point of view, in the first scenario, all the subjects’ true scores are concentrated around the center of the rating scale, and the raters are quite homogeneous in their severity and reliability. The heterogeneity of the subjects and the raters is only at the individual level since they are not nested with clusters. Under the second scenario, we introduce heterogeneity at the population level as both subjects and raters are assigned to different clusters, respectively. Here, we mimic the case in which subjects are clustered within two different levels of true score (e.g., low versus high proficiency level), and raters are either systematically slightly more lenient and reliable or more severe and less reliable. Under the third scenario, subjects and raters are assigned, respectively, to three poorly separated clusters. This results in a highly negatively skewed distribution for the subjects’ true score and a multimodal distribution for the raters’ systematic bias. Figures 3 and 4, Figure C.1 in Appendix C, and Figures 1, 2, and 3 in the Supplementary Material show the respective true densities and the empirical distributions of the generated ratings.

Figure 3 Average estimated density across 10 independent datasets under different scenarios. The columns indicate the cardinality of $|\mathcal {R}_i|=\{2,4\}$ : left and right, respectively; the rows indicate bimodal or multimodal scenario: first and second row, respectively. The solid red lines indicate the true densities; the solid black line and the shaded grey area indicate, respectively, the point-wise mean and $95\%$ quantile-based credible intervals; the density implied by the BP model (black dotted lines).

Figure 4 Top row: empirical distribution of the data (red solid line) and empirical distribution of replicated data (black solid lines) from the respective BNP and BP posterior distributions (left and right columns, respectively). Middle and bottom row: Test statistics computed on the data (red solid line) and histograms of those computed on replicated data.

Ten independent data sets are generated under the six scenarios resulting from the $2 \times 3$ design, for each data set, the standard parametric (BP), the semi-parametric (BSP), and the nonparametric (BNP) models are fitted.

Model recovery assessment. Parameter recovery performance is assessed through the Root Mean Square Error (RMSE) and the Mean Absolute Error (MAE) computed, respectively, as the root mean square difference and the mean absolute difference between the posterior mean and the true value of the parameters across data sets. For the subject and raters specific parameters, i.e., $\{\theta _i\}_1^I$ , $\{\tau _j, 1/\sigma ^2_j\}_1^J$ , RMSE, and MAE are average both across individuals and data sets.

For the sake of comparison across different scenarios, we report the standardized version of both indices (S-RMSE, S-MAE) for the structural parameters. More precisely, those related to $\mu $ and $\mu _G$ are divided by the mean value of the rating scale, i.e., 50; those regarding $\omega ^2$ , $\omega ^2_G$ , $\phi ^2$ , $\phi ^2$ , $\tilde {\sigma }$ , $\tilde {\sigma }_H$ , and the $ICC_A$ are divided by their true value.

The models’ performance in recovering the density distributions of individuals’ specific parameters is evaluated through visual inspection. We give an example of how different densities might lead to very different conclusions on the data generative process (Gelman et al., Reference Gelman, Carlin, Stern, Dunson, Vehtari and Rubin2013; Paganin & de Valpine, Reference Paganin and de Valpine2024; Steinbakk & Storvik, Reference Steinbakk and Storvik2009). Specifically, we draw new replications from the respective posterior predictive distributions and compare these samples to the original data. If the models capture relevant aspects of the data, they should look similar, and replications should not deviate systematically from the data. We measure discrepancy in central asymmetry through the statistic ${T_1(y,\mu _G)=|y_{.25}-\mu _G|-|y_{.75}-\mu _G|}$ , where $y_{.25}$ and $y_{,75}$ are the first and the third quartile, and in the left tail weight by the statistics $T_2(y)=min(y)$ .

8.1 Categorical modeling

Modeling rating $Y_{ij}$ We assume that the observed ordinal rating $Y_{ij}\in \{1,\dots , K\} \subset \mathbb {N}$ is generated by an underlying unobserved normally distributed variable $Y_{ij}^{\ast}$ (Jöreskog & Moustaki, Reference Jöreskog and Moustaki2001) and that we observe $Y_{ij}=k$ if $\delta _{k-1}<Y_{ij}^{\ast}\leq \delta _k$ ; $\delta _0=-\infty < \delta _1 < \dots , <\delta _K=+\infty $ are ordered thresholds over the underlying response variable distribution and are equal across raters. The underlying variable $Y_{ij}^{\ast}$ might be interpreted as a latent rating or the original continuous rating before the coarsening procedure. The conditional probability that $Y_{ij}=k$ is

(43) $$ \begin{align} \mathbb{P}[Y_{ij}=k|\theta_i, \tau_j,\sigma_j,\delta_{k},\delta_{k+1}] = \Phi \left( \frac{\delta_{k+1} - \theta_i - \tau_j}{ \sigma_j} \right) - \Phi \left( \frac{\delta_{k} - \theta_i - \tau_j}{ \sigma_j} \right), \end{align} $$

for $i = 1,\dots ,I; \;\; j \in \mathcal {R}_i$ . Additional considerations on the interpretation of $\sigma _j$ under this formulation are given in the Supplementary Material.

Identifiability issues. Under this parametrization, we need to put additional constraints for identifiability purposes since the underlying response variables’ mean and variance are freely estimated (DeYoreo & Kottas, Reference DeYoreo and Kottas2018; Kottas et al., Reference Kottas, Müller and Quintana2005). Two thresholds (e.g., $\delta _1$ , $\delta _{K-1}$ as proposed by Song et al., Reference Song, Lu, Cai and Ip2013) have to be fixed in advance, as it is common in multi-group analysis (Lockwood et al., Reference Lockwood, Castellano and Shear2018). From a statistical perspective, we note that each rater might be seen as a group of observations (Papaspiliopoulos et al., Reference Papaspiliopoulos, Stumpf-Fétizon and Zanella2023). Moreover, an SC-DPM prior has to be placed on the subject’s true score $\{\theta _i\}_1^I$ to fix their mean and resolve identifiability issues (Gelman et al., Reference Gelman, Hwang and Vehtari2014), as a by-product under the parameter-expanded specification, equation (43) becomes

(44) $$ \begin{align} \mathbb{P}[Y_{ij}=k|\theta_i^{\ast}, \tau_j^{\ast},\sigma_j,\delta_{k},\delta_{k+1}] = \Phi \left( \frac{\delta_{k+1} - \theta^{\ast}_i - \tau^{\ast}_j}{ \sigma_j} \right) - \Phi \left( \frac{\delta_{k} - \theta^{\ast}_i - \tau^{\ast}_j}{ \sigma_j} \right), \end{align} $$

for $i = 1,\dots ,I; \;\; j \in \mathcal {R}_i$ . Whenever $K=2$ , i.e., dichotomous rating scale, $\{\sigma _j\}_1^J$ can not be identified and need to be fixed in advance, e.g., $\sigma _j=1$ , $j=1,\dots , J$ , which implies assuming raters to be equally reliable (Lockwood et al., Reference Lockwood, Castellano and Shear2018).

Generalized ICCs. Under this model specification, the ICCs computed according to propositions 1 and 2 are generalized ICCs that indicate the polychoric correlation between two latent ratings (Jöreskog, Reference Jöreskog1994; Uebersax, Reference Uebersax1993). For instance, proposition 1 implies here:

(45) $$ \begin{align} ICC_{j,j'}^{\ast} = Corr\left(Y_{ij}^{\ast}, Y_{ij'}^{\ast} |G,H, \sigma_j^2, \sigma_{j'}^2 \right) = \frac{\omega^2_G}{\sqrt{\omega^2_G +\phi^2_H +\sigma_j^2}\sqrt{\omega^2_G +\phi^2_H +\sigma_{j'}^2} }, \end{align} $$

where $ICC_{j,j'}^{\ast}$ indicates the conditional pairwise polychoric correlation between the latent ratings given by raters $j \neq j'$ to subject i. Similar considerations might be extended to propositions 2 and 3. As a by-product, the $ICC_A^{\ast}$ is the lower bound of the expected polychoric correlation between the latent ratings $Y_{ij}^{\ast}$ and $Y_{ij'}^{\ast}$ , with $j \neq j'$ :

(46) $$ \begin{align} ICC_A^{\ast} \leq \mathbf{E}[Corr\left(Y_{ij}^{\ast}, Y_{ij'}^{\ast}|G,H \right)] = \mathbf{E}[ICC^{\ast}|G,H]. \end{align} $$

8.2 Posterior computation

A data augmentation procedure may simulate the underlying response variables (Albert & Chib, Reference Albert and Chib1993). The underlying continuous ratings $Y_{ij}^{\ast}$ , $i=1,\dots , I$ , $j \in \mathcal {R}_i$ are sampled:

$$\begin{align*} Y^{\ast}_{ij}| \cdot \overset{\mathrm{ind}}{\sim} N(\theta^{\ast}_i-\tau_j,\sigma^2_j) \times I(\delta_{k-1} < Y^{\ast}_{ij} \le \delta_k), \quad k=1,\dots,K. \end{align*}$$

Here $I(\cdot )$ is an indicator function. Following Albert & Chib (Reference Albert and Chib1993) the conditional posterior distribution of the $K-3$ freely estimated thresholds,e.g., $\delta _2,\dots ,\delta _{K-2}$ might be seen to be uniform on the respective intervals:

$$ \begin{align*} \delta_k|\cdot &\overset{\mathrm{ind}}{\sim}& U(max\{max\{Y_{ij}^{\ast}: Y_{ij}=k\}, \delta_{k-1}\}, min\{min\{Y_{ij}^{\ast}: Y_{ij}=k+1\}, \delta_{k+1}\} ), \end{align*} $$

here $U(\cdot )$ stands for uniform distribution.

All the other parameters are updated according to the posterior sampling scheme detailed in Section 3.1 of Supplementary Material and the post-process transformation outlined in Section 5.3 needs to take into account the double-centering. After computing $\mu _G$ and $\eta _H$ according to 16 and 17 for each iteration, the samples of $\mu _0,\mu _G,\{\theta _i\}_1^I$ , and $\{\tau _j\}_1^J$ are computed as follows:

$$ \begin{align*} \mu_0^{\ast} &= \mu_0 -\mu_G + \eta_H, \\ \theta_i^{\ast} & = \theta_i - \mu_G + \eta_H, \quad \quad \text{for} \;\; i=1,\dots,I; \\ \tau_j^{\ast} & = \tau_j - \eta_H + \mu_G, \quad \quad \text{for} \;\; j=1,\dots,J. \end{align*} $$

8.3 Generated and real coarsened ratings analysis

In this section, we present the analysis of real and generated coarsened ratings and compare the results with those presented in Sections 6 and 7. For the real data, we deliberately coarsened the original continuous ratings analyzed in Section 7 into $K=4$ ordered categories according to the following cutoffs: $\delta _1=20$ , $\delta _2=30$ , $\delta _3=40$ . The fit of the BP, BSP, and BNP models to the data are compared according to the WAIC for ordered data discussed in the Supplementary Material.

We performed a simulation study to assess the accuracy of the BNP and the BP versions for ordered ratings. More specifically, the same data sets generated under the bimodal scenarios in Section 6 are coarsened and considered for this study. We coarse these ratings into $K=4$ ordered categories according to three consecutive cutoffs: $\delta _1=35$ , $\delta _2=50$ , $\delta _3=75$ . The same parameter recovery assessment procedure detailed in Section 6 is consistently used here.

Table 5 RMSE and MAE of individuals parameters across bimodal scenarios with coarsened ratings

Note: The bold text indicates the average, most accurate estimates.

Table 6 The WAIC is reported for each of the fitted models: BNP, BP, and BSP; the pairwise WAIC difference ( $\Delta WAIC$ ) between the model with the best fit and each other is reported

In real context, the cutoffs of the coarsening procedure are generally known since the continuous rating scale is deliberately broken down into a small number of consecutive intervals and raters are explicitly asked to coarse their ratings accordingly (Peeters, Reference Peeters2015; van Praag et al., Reference van Praag, Hop and Greene2025). For example, on a 1–100 continuous scale, they might be asked to indicate which of the following intervals each subject’s score falls into: (1–25), (25–50), (50,75) or (75,100). On the contrary, when ratings are directly given on an ordinal scale, the categories’ labels are not necessarily associated with any continuous scale intervals (e.g., “poor”, “acceptable”, “good”, “very good”). In these scenarios, we consider the observed ordered ratings as coarsened representations of underlying continuous values according to some unknown consecutive cutoffs. In the first case, this coarsening process is factual; in the second, it is merely assumed. However, since in both cases at least two cutoffs need to be fixed for identification purposes, we decide to fix $\delta _1$ and $\delta _3$ to the true values and let the model estimate $\delta _2$ , both for real and generated data.

Results. The total computational elapsed times for the BP, BSP, and BNP models were roughly similar to those of previous Sections. Upon graphical inspections of the MCMC chains and diagnostics, no convergence or mixing issues emerged for both generated and real data. Table 6 gives the WAIC indices for each fitted model and suggests that the BNP model provides the best fit to the data. Based on this model comparison procedure, we focus on the results from the BNP model. As shown in Table 7, the estimates are equivalent to those obtained under the continuous BNP model presented in Section 7. We note that the only notable difference concerns the point estimate of the subjects’ clustering structure (Figure 7). In this case, they are clustered into two (instead of four) subjects’ groups.

Table 7 Posterior mean and $95\%$ quantile-based credible intervals of the estimated structural parameters of the BNP model are reported

Results from generated data suggest that the BNP model provides more accurate estimates of subjects’ and raters’ specific parameters and overcomes the BP model. The only exception is observed for the rater-specific reliability parameter $1/\sigma ^2$ under the scenario $|\mathcal {R}_i|=4$ ; here, the BP model overcomes our proposal. Under the standard parametric model, we only have two population parameters $\gamma $ and $\beta $ (i.e., $(\gamma _j,\beta _j)=(\gamma ,\beta )$ , for $j=1,\dots , J$ ) and, as a consequence, more information is available for their estimation. This might result in a faster accuracy improvement of this model for this set of parameters as the ratio of students per rater increases. The comparison between the $RMSE$ and the $MAE$ of Tables 1 and 5 suggests that the estimates of both the BP and BNP models degrade with coarse data. The same trend emerged regarding the structural parameters and the densities; we report these results in the Supplementary Material.

Figure 7 The estimated densities of the subject’s true score $\theta $ , rater’s systematic bias $\tau $ and the residual term $\epsilon $ are reported; the black solid lines and the shade grey areas indicate the pointwise posterior mean and $95\%$ quantile-based credible intervals of the respective densities. Bottom-right figure shows the posterior distribution of the $ICC_A$ , the black solid and dotted lines indicate, respectively, the $95\%$ credible interval and the posterior mean. The rugs at the margins of the first three figures indicate the clustering of individuals.