Modeling random and non-random decision uncertainty in ratings data: a fuzzy beta model

Modeling human ratings data subject to raters’ decision uncertainty is an attractive problem in applied statistics. In view of the complex interplay between emotion and decision making in rating processes, final raters’ choices seldom reflect the true underlying raters’ responses. Rather, they are imprecisely observed in the sense that they are subject to a non-random component of uncertainty, namely the decision uncertainty. The purpose of this article is to illustrate a statistical approach to analyse ratings data which integrates both random and non-random components of the rating process. In particular, beta fuzzy numbers are used to model raters’ non-random decision uncertainty and a variable dispersion beta linear model is instead adopted to model the random counterpart of rating responses. The main idea is to quantify characteristics of latent and non-fuzzy rating responses by means of random observations subject to fuzziness. To do so, a fuzzy version of the Expectation–Maximization algorithm is adopted to both estimate model’s parameters and compute their standard errors. Finally, the characteristics of the proposed fuzzy beta model are investigated by means of a simulation study as well as two case studies from behavioral and social contexts.


Introduction
In social and behavioral research, satisfaction surveys, aptitude and personality testing, demographic inquiries, and life quality questionnaires are widespread tools to collect data involving subjective evaluations, agreements, and judgments. In a typical social survey, a set of questions is administered to a sample of participants and they are asked to express the extent of their agreement on bounded discrete or continuous rating scales (Aiken 1996;Miller and Salkind 2002). Traditionally, ratings data have been collected by means of pencil and paper questionnaires although more flexible implementations, such as online questionnaires, have become increasingly popular. More recently, technological advancements have fostered the development of new rating tools which offer an accurate way to trace the rating process from its beginning to the final rating outcome (Schulte-Mecklenbeck et al. 2011). Hence, unlike standard rating tools, these techniques allow researchers to collect a richer variety of participants' data, including temporal unfolding of ratings and decision uncertainty Freeman and Ambady 2010). As several scholars have pointed out, decision uncertainty can be referred to a subject-specific cognitive component of the rating process designed to construct a coherent mental representation of the question being rated (Ülkümen et al. 2016). In this sense, it reflects the subjective interplay of decisional and emotional components which contribute to the final rating response (Kahneman and Tversky 1982). As such, when appropriately utilized, this source of within-subject heterogeneity can reveal more about ratings then standard crisp responses. Interestingly, this type of non-random, systematic uncertainty has been extensively studied especially with regards to its effects on rating responses (e.g., see Saal et al. 1980). Indeed, it is widely recognized that ratings data often suffer from lack of accuracy, for instance because of social desirability (Furnham 1986), faking behaviors (Lombardi et al. 2015), personality (Muthukumarana and Swartz 2014), response styles (Eid and Zickar 2007), and violations of rating rules (Iannario 2015;Preston and Colman 2000;Rabinowitz et al. 2019). These issues have not only been recognized as important by applied statisticians working with ratings data but also by several researchers working in fields like applied econometrics (e.g., see Angel et al. 2019;De Bruin et al. 2011;Zafar 2011), metrology (e.g., see Pendrill and Petersson 2016;Pendrill 2014), and risk analysis (e.g., see Slovic et al. 2004).
Modeling ratings data is a relevant problem in applied statistics. Commonly used methods to handle with discrete or continuous rating data include generalized linear models (GLMs) (McCulloch 2000), beta regression models (Ferrari and Cribari-Neto 2004;Migliorati et al. 2018;Ospina and Ferrari 2012), and combination of uniform and shifted binomial (CUBE) models (Golia 2015;. These models typically represent mean and dispersion components as a linear or non-linear function of some external covariates, which are intended to explain the observed heterogeneity of the ratings data. Although some of these methods also allow for disentangling individual indecision and heterogeneity of responses induced by the presence of subgroups (e.g., CUBE), they are mainly intended to work with data represented as a crisp collection of responses and do not account for non-random decision components of rating process. The same applies with more general approaches to analyse within-subject heterogeneity such as random-effects and errors-in-variables models (Feng et al. 2018) which do not deal with non-random components of uncertainty. As a consequence, decision uncertainty underlying participant's rating process is not formally represented in these models.
In this contribution we propose a novel method for analysing continuous bounded ratings data that are characterized by non-random and systematic decision uncertainty. In particular, we propose a variable dispersion beta linear model which is generalized to cope with data contaminated by subjective uncertainty. We represent decision uncertainty in the framework of fuzzy data modeling, where crisp ratings data are equipped with non-random systematic uncertainty via normalized set functions Kruse and Meyer 1987). In this setting, maximum likelihood estimation and inference are carried out through the Expectation-Maximization algorithm adapted for the case of fuzzy data (Denoeux 2011;Su et al. 2015). It should be stressed that using beta linear models allows for flexibility in modeling and analysing continuous ratings data, while still retaining simplicity in estimates model's parameters (Algamal 2019;Canterle and Bayer 2019;Zeileis et al. 2010). Similarly, despite the fact that many formal theories have been proposed to deal with subjective uncertainty (e.g., soft sets, rough sets. See Lin and Cercone 2012;Liu 2010), fuzzy set theory offers a good compromise in terms of accuracy and computational costs and benefit from a long tradition of works in statistics (e.g., see Buckley 2006;Chukhrova and Johannssen 2018;Gebhardt et al. 1998;González-Rodríguez et al. 2006).
The remainder of the article is organized as follows. Section 2 briefly describes the basic characteristics of fuzzy data together with their interpretation in terms of decision uncertainty. Section 3 exposes the variable dispersion fuzzy beta model, parameters estimation, and model evaluation. Section 4 reports results of a short simulation study performed to evaluate the finite sample properties of the fuzzy beta model and the consequences of neglecting decision uncertainty on parameters estimation. Section 5 describes an application of the new approach to two case studies involving ratings data from risk-taking behaviors (application 1) and customer satisfaction (application 2). Finally, Sect. 6 concludes the article providing final remarks and suggestions for future extensions. All the materials like datasets and R-scripts used throughout the article are available to download at https:// github. com/ antca lcagni/ fuzzy ratin gbeta.

Data representation
In this section we briefly review some concepts and terminology related to fuzzy numbers, fuzzy probability, and fuzzy ratings data.

Fuzzy numbers and fuzzy probability
A fuzzy subset Ã of a universal set A is defined by means of its characteristic function Ã ∶ A → [0, 1] . It can be easily described as a collection of crisp subsets called -sets, i.e. Ã = {x ∈ A∶ � A (x) > } with ∈ (0, 1] . If the -sets of Ã are all convex sets then Ã is a convex fuzzy set. The support of Ã is Ã 0 = {x ∈ A∶ � A (x) > 0} and the core is the set of all its maximal points, i.e. Ã c = {x ∈ A∶ � A (x) = max y∈A � A (y)} . In the case max x∈A Ã (x) = 1 then Ã is a normal fuzzy set. If Ã is a normal and convex subset of ℝ then Ã is a fuzzy number (Buckley 2006). Broadly speaking, fuzzy sets can be conceived as subsets of ℝ where their Boolean characteristic function A (x) , ∀x ∈ A ⊂ ℝ , has been generalized to the real interval [0, 1]. The class of all normal fuzzy numbers is denoted by F(ℝ) . Fuzzy numbers can conveniently be represented using parametric models, through which Ã is represented by means of few real parameters. Hence, we can define families of fuzzy numbers indexed by some scalars, such as m (mode) and s (spread/precision), which include a number of shapes like triangular, trapezoidal, gaussian, and exponential (Buckley 2006). Relevant classes of parametric fuzzy numbers are the so-called LR-fuzzy numbers (Dubois and Prade 1978) and their generalizations Dombi and Jónás 2018). In this setting, a set of operators and algebras have also been defined for fuzzy numbers, which extend traditional calculus to fuzzy numbers as well (Chwastyk and Kosiński 2013). A broader class that encompasses a wide range of fuzzy numbers is the so-called beta fuzzy number (Alimi 2003;Baklouti et al. 2018;Stein 1985): where x l , x u , a, b ∈ ℝ , with x l and x u being the lower and upper bounds of the set, and m the mode of the fuzzy set. This type of fuzzy numbers uses Beta functions to approximate many regular shapes such as triangular, trapezoidal or Gaussian. Likewise for LR-fuzzy numbers, beta fuzzy numbers can be defined in terms of mode m ∈ ℝ and spread/precision s ∈ ℝ + parameters. In particular, let x l = 0 and x u = 1 without loss of generality. Then Eq. (1) can be re-arranged as follows: with C being a constant ensuring Ã is still a normal fuzzy set: (1) (2) Interestingly, the re-parameterized Beta fuzzy number resembles the shape of Beta density distribution written using the PERT representation (Vose 2008). It should be noted that, like triangular or trapezoidal fuzzy numbers, also beta fuzzy numbers can be adopted to represent rating data variables. Indeed, the beta-like shape is particularly flexible to modeling some features of bounded rating data, such as asymmetry and tail flatness (e.g., see Migliorati et al. 2018).
When a probability space is defined over the reals, the probability of a fuzzy set ℙ(Ã) can also be defined. Over the years, there have been various attempts to define the probability of a fuzzy set in terms of expected value of its membership function (Zadeh 1968), conditional probability of prior information (Coletti and Scozzafava 2004;Singpurwalla and Booker 2004), imprecise probability (Augustin et al. 2014), fuzzy numbers (Hesamian and Shams 2017), and likelihood induced by random events (Cattaneo 2017). Following the findings of Denoeux (2011), in this contribution we adopt Zadeh's definition of fuzzy probability (Zadeh 1968). In particular, let (ℝ, B, ℙ) be a probability space. Then, ℙ(Ã) is defined as follows: with Ã being Borel measurable. In this context, two fuzzy sets Ã and B are said independent w.r.t. to ℙ , if ℙ(ÃB) = ℙ(Ã) ⋅ ℙ(B) , with the fuzzy product being defined as ÃB(x) =Ã(x) ⋅B(x) (Zadeh 1968). The conditional probability of two independent fuzzy events is with ℙ(B) > 0 . Note that one can also obtain the conditional probability between a crisp set A and a fuzzy set B as a special case of Eq. (4). If a discrete or continuous random variable Y is defined over B , then fuzzy probability can be generalized accordingly. For instance, denoting with f Y (y; ) the probability density of Y, then Similarly, when a sample of n independent observations = (y 1 , … , y n ) from Y 1 , … , Y n is available, the likelihood of the sample can be generalized as follows: where the definition ̃ = ∏ n i=1 ỹ i (y) has been used for the joint fuzzy set (Gebhardt et al. 1998). Further details about fuzzy generalization of likelihood functions, fuzzy random variables, and fuzzy probability space can be found in Cattaneo (2017), , Denoeux (2011), Gebhardt et al. (1998, Gil et al. (2006).
We interpret fuzzy data in the context of random variables following the epistemic viewpoint on fuzzy set theory . In particular, for a fuzzy set Ã , Ã (Y = y) is interpreted as the possibility that the crisp event Y = y has to occur. Indeed, Ã (Y = y) ∈ (0, 1) can be conceived as a graded plausibility about the occurrence of the event Y = y , with Ã (Y = y) = 1 indicating the fact that Y = y is fully possible. By contrast, Ã (Y = y) = 0 indicates that Y = y is not possible at all. Hence, fuzzy sets can intuitively be viewed as graded constraints on crisp random variables. In this way, the randomness due to the data generation process and the fuzziness due to observer's state of knowledge can be analysed simultaneously by means of a common statistical representation. As a remark, note that in this setting Ã (Y = y) is thought as being the consequence of a two-step generation process, in which a realization y is drawn from Y first and then a fuzzy set Ã is used to encapsulate the uncertainty about Y = y in terms of possibility distribution. Hence, only the first stage is a random experiment whereas the second stage is a non-random fuzzification of the outcomes being realized.

Fuzzy scaling and fuzzy data
Since the seminal work of Hesketh et al. (1988), fuzzy sets have been extensively used in the context of ratings data (for a review, see Lubiano et al. 2016;Calcagnì et al. 2021). Although several formats have been proposed for implementing fuzzy ratings tools (e.g., conversion scales, direct rating scales, implicit rating scales), all of them share the same idea that ratings data cannot be coerced into crisp numbers without a certain loss of information. Except for the case of dichotomous ratings, polytomous responses often show some degrees of imprecision and fuzziness, which is essentially due to participants' rating processes (Kahneman and Tversky 1982). In general, ratings data can be enriched by including information related to participants' response process and fuzzy conversion or fuzzy rating systems can be adopted to this purpose. In what follows, we will briefly review both the approaches to fuzzy rating.
Fuzzy conversion scales aim at turning standard crisp ratings into fuzzy numbers through the adoption of user-defined fuzzy systems (e.g., see Vonglao 2017) or statistically-oriented procedures (e.g., see Yu and Wu 2009). In general, a typical implementation of a user-defined fuzzy conversion scale, a fuzzy system relates a space of crisp responses (input) to a space of fuzzy numbers (output). Next implication rules mapping both input and output sources are finally established. For instance, in the simplest case of a 5-point scale, the crisp response Y = 2 activates the fuzzy sets of the output space via an IF-THEN implication rule. On the basis of the implemented fuzzy system (e.g., Mamdani, Sugeno), the output can result in the activation of one or two fuzzy sets and the final fuzzy response-which is given by the intersection of the activated fuzzy sets-will represent a more certain or less certain response, respectively.
Fuzzy rating scales adopt computerized interfaces through which the rater's response process is directly mapped to fuzzy numbers (e.g., see Hesketh et al. 1988;de Sáa et al. 2014). For instance, in the direct rating based implementation (de Sáa et al. 2014), raters are asked to respond by drawing fuzzy sets according to their perceived uncertainty. In this case, the rating procedure proceeds as follows. First, raters draw an interval on a pseudo-continuous graphical scale, which represents the set of admissible responses compatible with the assessment of the item being rated. Then, a degree of confidence is expressed by drawing another interval around the previous interval response. Finally, both the intervals are interpolated to form the final fuzzy responses (e.g., trapezoidal or triangular). In a similar way, fuzzy indirect scales adopt a computerized graphical interface to get responses from raters. However, in this case, they are not asked to directly draw their responses; rather, fuzzy sets are build from a set of implicit measures associated to the final crisp response (e.g., response time). For instance, the DYFRAT scale ) uses a set of biometric measures associated with the cognitive response process (i.e., response time, computer-mouse trajectories) to derive the rater's fuzzy response. In particular, raters are presented with a pseudo-circular scale with M levels and-similarly to the case of crisp Likert-type scales-they are asked to choose which of these levels is the most appropriate to represent their response for a given item being rated. Meanwhile, the system records the streaming coordinates of the computer mouse cursor (at a fixed sampling rate) and the elapsed time since the beginning of the response trajectory. Finally, both information are used to define a fuzzy response according to the following rationale: (1) computer-mouse trajectories are projected onto a linear scale and the histogram of linearized radians is used to derive the fuzzy set (e.g., triangular, beta) in terms of its support and core, (2) response times are used to intensify or reduce the fuzziness of the sets by means of linguistic quantifiers (e.g., the longer the response time, the higher the fuzziness of the set).
In the context of fuzzy scaling, fuzzy numbers are used as a formal representation for rating responses involving individual-based judgments, attitudes, and opinions. Although fuzzy conversion and rating scales differ in the way they derive fuzzy responses, both aim at providing a model for the fuzziness or imprecision which is present in the rating process Y ∼ f Y (y; ) . As for the more general case of LR fuzzy numbers, the parameters of a beta fuzzy number can be linked to the rating process as follows. First, consider a continuous rating scale bounded on a subset (y l , y u ) of reals. Then, m represents the most plausible final rating choice Ỹ (Y = y) = 1 , s is the precision of m such that smaller values indicate larger levels of hesitation in the rating choice, and Ỹ conveys the overall decision uncertainty in terms of fuzziness (the larger the fuzziness, the highest the decision uncertainty). Note that, ideally, if there was no subjective uncertainty, then the fuzziness would tend to zero and true rating realizations would be precisely observed (i.e., Y = y = m ). In this case, there would not need to represent ratings as fuzzy data.

Variable dispersion beta model for fuzzy ratings
In this section we illustrate our proposal to analyse continuous bounded ratings data in situations with decision uncertainty. Hereafter, ratings data will be considered scaled into the real subset (0, 1) without loss of generality.

Model
Let = (y 1 , … , y i , … , y n ) ∈ (0, 1) n be a sample of n observations from Beta distributed independent random variables Y 1 , … , Y i , … , Y n with density: with ∈ (0, 1) n being the n × 1 vector of location parameters and ∈ (0, ∞) n the n × 1 vector of precision parameters (Ferrari and Cribari-Neto 2004). The sequence Y 1 , … , Y i , … , Y n models the ratings for each of the n participants, with the convention that {0, 1} represent the lower and upper bounds of the rating domain, respectively. In order to account for heterogeneity and non-constant variance in rating responses, location and dispersion parameters can be non-linearly re-written using monotonic and twice differentiable link functions, mapping the support (0, 1) into ℝ , as follows: where and are n × J and n × H matrices of known continuous or categorical covariates, with and being vectors of appropriate order containing unknown parameters. The functions g 1 (⋅) and g 2 (⋅) can be chosen among a variety of link functions (e.g., logit, probit, log. See McCulloch 2000). Two typical choices are the logit and logarithm functions, which yield to: Under Eq. (8), the log-likelihood function for the variable dispersion beta model is: with (⋅) being the Euler gamma function and = ln , ln ( − ) the sufficient statistics for the inference on = ( , ).
In light of the data representation adopted in this work, decision uncertainty is treated as a systematic and non-random component which occurs after the sampling process ∼ f ( ; , ) has been realized. This leads to a situation where the sample cannot be precisely observed and a collection of fuzzy data ̃ is instead available. When fuzzy data are represented as beta fuzzy numbers, then with and being n × 1 vectors of modes and precisions/spreads for the fuzzy observations. Turning Eqs. (3) into (6), the joint density of ̃ can be written as: Modeling random and non-random decision uncertainty in ratings… where the joint fuzzy set ̃ = ∏ n i=1 ỹ i (y) has been factorized in terms of product. The log-likelihood of the model under fuzzy observations can analogously be obtained using Eq. (10). Note that, in this representation the vectors and of the fuzzy sets enter the model as observed quantities whereas the parameters and still remain non-fuzzy quantities.

Parameter estimation
To provide estimates for = ( , ) in the context of fuzzy ratings data, one can maximize the log-likelihood l( ;̃ ) function, which is obtained by Eq. (10). This would require an iterative procedure, alternating between the numerical computation of the integral and the maximization of the function. However, to avoid the problem of approximating integrals in Eq. (10) and have a way to compute standard errors consistently, we will use a variant of the Expectation-Maximization algorithm generalized for the case of fuzzy data (Denoeux 2011). As for the standard EM algorithm, the fuzzy-EM version at the kth iteration alternates between the E-step, which involves the computation of the expected complete-data loglikelihood using (k−1) , and the M-step, which instead maximizes the expected complete-data log-likelihood w.r.t. to (k) . These steps generate a non-decreasing sequence of lower bounds for the maximization of the observed-data log-likelihood l( ;̃ ) (for formal details, see Denoeux 2011). In the fuzzy-EM variant, the complete-data log-likelihood is that obtained if was precisely observed (see Eq. 9). Therefore, given the (k − 1) th estimates � = (k−1) the E-step for the kth iteration of the algorithm consists in the computation of the following quantity: where C contains all the terms of Eq. (9) that do not involve the random quantities to be filtered. To compute the conditional expectations of the E-step, note that Y i | |ỹi is a random variable conditioned on fuzzy events and its density can be obtained using Eq. (4) simplified for the case where A is a crisp event: Under the case of beta fuzzy numbers and up to a normalization constant, the conditional density f Y i |ỹ i (y; i , i ) corresponds to a beta density with parameters given as a function of fuzzy data and crisp parameters of the complete-data density: Then, the first expectation ′ ln Y i | |ỹi can be approximated via Taylor expansion around ′ Y i | |ỹi as follows: Similarly, the second expectation � ln(1 − Y i ) | |ỹi is obtained by symmetry of the beta function: Once expected values are computed, the M-step of the algorithm involves the maximization of Q , ′ with respect to the elements of and can be performed by plugging-in Eqs. (14)-(15) into the log-likelihood Eq. (11). The simultaneous score equations for M-step are as follows: where and are defined as in Eq. (8) whereas the terms  (17) numerically, for instance using Gauss-Newton method.

Standard errors, diagnostics, and inference
Standard errors for ̂ and ̂ can be computed using the observed information matrix Î from the Hessian of the maximum likelihood estimates ̂ obtained solving Eqs.
(16)-(17) for and . In the context of EM algorithm, Î can be approximated using the empirical observed information matrix (Meilijson 1989), as follows: where is the score vector for the ith observation calculated at ̂ and ̂ . The standard errors are then calculated as usual: Note that this approximation avoids the computation of Hessian of the completedata log-likelihood as it solely uses the score equations where the unobserved vector is replaced by * (1) and * (2) . Alternatively, standard errors can also be obtained via non-parametric bootstrap (McLachlan and Peel 2004). An important quantity commonly used to assess the quality of the estimated model is that involving standardized residuals which, for the case of fuzzy data, can be generalized as follows: where ̂i = (1 + exp( î) ) −1 and ̂= exp( î) whereas � y i (̂i) is the fuzzy membership computed for the predicted quantity ̂i . In general, diagnostics for the model can be performed by plotting, for instance, r 1 , … , r n against the indices of the observations in order to check for particular trends or patterns in the predicted data. Similarly to generalized and beta linear models, also the overall fit of the estimated fuzzy beta model with respect to the observed fuzzy data can be assessed by means of pseudo-R 2 indices, which generalize the standard residuals-based R 2 indices (Ferrari and Cribari-Neto 2004;Veall and Zimmermann 1994). In this context, we resorted in applying a likelihood-ratio based pseudo-R 2 index (Aldrich and Nelson 1984;Veall and Zimmermann 1994) where = 2 l(̂; ) M 1 − l(̂; ) M 0 and = 1 n l(̂; ) M 0 . The likelihood-ratio based pseudo-R 2 is normalized in [0, 1] and approximates the relationship between the likelihood-ratio statistic and the R 2 index in the linear regression setting (Veall and Zimmermann 1994).
Finally, likewise for the non-fuzzy case, inference on and can be performed using maximum-likelihood (ML) theory (Denoeux 2011;McLachlan and Peel 2004) and, consequently, hypothesis testing on model's parameters can be performed using fuzzy version of likelihood ratio test (e.g., see Berkachy and Donzé 2019;Najafi et al. 2010). In this context, as for the maximum-likelihood theory under the EM procedure, inferential results are based on the asymptotic properties of ML-based estimators. For further details, we refer the reader to Denoeux (2011), McLachlan andPeel (2004), and Berkachy and Donzé (2019).

Simulation study
The aim of this study is twofold. First, we will evaluate the performances of EM estimators for location and precision parameters for the fuzzy beta linear model. Second, we will assess whether standard methods, such as fixed/random-effects beta linear models, perform as good as the proposed method if applied on defuzzified data. Although the EM algorithm for fuzzy data has been validated elsewhere (e.g., see De Bruin et al. 2011;Su et al. 2015), in the present study we have preferred to evaluate the performances of the fuzzy-EM procedure to further provide converging results. The whole simulation procedure has been performed on a (remote) HPC machine based on 16 cpu Intel Xeon CPU E5-2630L v3 1.80 GHz, 16 × 4 GB Ram whereas computations and analyses have been done in the R framework for statistical analyses.
Design The design involved three factors, namely (1) n ∈ {50, 100, 250, 500} , Modeling random and non-random decision uncertainty in ratings… (a) n l ×j k = [ n l , n l ×j k−1 ] and n l ×h m = [ n l , n l ×h m−1 ] , where n l ×j k−1 and n l ×h m−1 were drawn from Unif (1, 5) (b) location and precision terms were computed as follows: (c) crisp data underlying fuzzy observations were generated according to the variable dispersion beta linear model y i ∼ Beta( i i , i − i i ) , i = 1, … , n l (d) fuzzy data were generated by making n l imprecise via a two-step data-generation process (Quost and Denoeux 2016;Su et al. 2014). First, spread components were generated as n l ∼ Gamma(1.025, 0.001) . Second, modes were generated by m i ∼ Beta(y i s i , s i − s i y i ) , i = 1, … , n l (e) parameters j k and h m were estimated using four methods: • fEM: expectation-maximization estimators for the fuzzy case. • dML: maximum-likelihood estimators on two type of defuzzified data computed using the centroid method y * i = ∫ 1 0 y ỹ i (y) dy and the first-maximum method y * i = sup y∈(0,1) ỹ i (y) for i = 1, … , n l . The ML procedure implemented in the R library betareg has been used in this case (Zeileis et al. 2010). • dREML: restricted maximum-likelihood estimators on defuzzified data obtained by treating fuzzy sets ỹ 1 , … , ỹ n as random effects. In particular, for each observation i = 1, … , n l , extremes of -sets y * i = [inf( c i ), sup( c i )] were used, with c i = {y ∈ (0, 1) ∶ ỹ i (y) ≥ } being the set obtained by cutting the ith fuzzy data at = [0.01, 0.35, 0.7] . In this case, estimates have been performed using the R library glmmTMB for random-effects beta linear models (Brooks et al. 2017).
Outcome measures For each condition of the simulation design, sample results were evaluated using bias of estimates and root mean square error.
Results Tables 1 and 2 show the results of the simulation study with regards to averaged bias and root mean square error. For the sake of clarity, results for the cases J = 2 and J = 4 were reported separately. To better interpret the results, it should be noted that conditions with H = 1 represent simplest cases where the precision term is held constant (the linear term for is a simple intercept model). By contrast, conditions with H = 3 represent those situations showing some levels of heterogeneity in the response variable (in this case the linear term for contains two slopes). We first consider the conditions with J = 2 . With respect to the parameters ̂ , all the methods showed negligible bias in estimating the location terms of the beta linear model both in the cases with H = 1 and H = 3 . However, unlike dML and dREML, the fEM solution achieved lowest RMSE. Considering the parameters ̂ , dML and dREML algorithms showed worse performances when compared to fEM both in the case of low ( H = 1 ) and high ( H = 3 ) model complexity. In particular, they showed larger bias in estimating the precision terms of the beta linear model, with bias being higher with increasing model complexity ( H = 3 ). A similar pattern was also observed for RMSE. In particular, when compared to fEM, dML and dREML showed larger values, with severity increasing over model complexity. Interestingly, although the condition with H = 3 represents a more complex situation with respect to parameters estimation, fEM outperformed the other methods. Results for the conditions with J = 4 largely resemble those obtained with H = 1 . Also in this case, the fEM algorithm showed better performances over dML and dREML. Finally, for each method we also computed overall indices of over/under-estimation r and r , as the ratio between the number of positive and negative bias, and the overall percentage of over-estimation p and p . Table 3 reports the results along with the overall RMSE for both the arrays of parameters. In general, fEM and dML (mode) showed negligible overestimation for 0 whereas dML (mean) and dREML tended to overestimate the true arrays of parameters. On the contrary, when estimating 0 , fEM tended to overestimate the true population parameters whereas dML and dREML tended to underestimate. On the whole, fEM showed less variable and less biased estimates of 0 and 0 than the other procedures adopted for defuzzified data.

Applications
In this section we will illustrate the application of the fuzzy beta model to two empirical studies involving fuzzy ratings data collected using two types of fuzzy rating scales. In particular, the first application concerns the analysis of risk-taking behavioral data collected by means of an indirect fuzzy rating approach . By contrast, the second application is about the analysis of customer satisfaction data collected by means of a direct fuzzy rating approach  (de Sáa et al. 2014). Note that the applications are chosen to review two of the most common methods for fuzzy scaling, namely fuzzy indirect (case study 1) and direct rating scales (case study 2).

Case study 1: Reckless-driving and risk-taking behavior in young drivers
Reckless-driving among young people is one of the major cause of mortality and injuries worldwide (Toroyan 2015). It is a very complex phenomenon involving a number of human and non-human factors, such as personality, cognitive styles, social context (e.g., family, peer group), infrastructures (e.g., roads, light), and cultures (e.g., see Bıçaksız and Özkan 2016;McNally and Bradley 2014;Scott-Parker et al. 2015;Taubman-Ben-Ari 2010). Several studies have recognized the role of subjective factors like sensation-seeking, normlessness, anxiety, aggressiveness, driving attitudes in determining risky behaviors (Bıçaksız and Özkan 2016).
Researchers have also assessed the contribution of parenting styles and peer relations to young drivers' intention to take risks (Taubman-Ben-Ari 2010, 2014). Because of its characteristics, assessing reckless-driving behaviors is a typical situation where self-reported measures can show some levels of decision uncertainty, which cannot appropriately be analysed using final crisp responses only. In this application, we consider a set of models where reckless-driving behaviors (rdb) were linearly predicted by the use of substances (drugs), driving anger (anger), and family climate (fcrs). In particular, we hypothesized that both the use of substances and driving anger would linearly increase the self-reported reckless-driving behaviors whereas family climate would instead acts by decreasing the amount of risky behaviors. Moreover, we also assessed whether the dispersion of fuzzy ratings data varied as a function of participants' characteristics such as gender (sex) and frequency of driving (driving_frequency). Data and measures A questionnaire survey was carried out on n = 69 young drivers in Trentino region (north-est of Italy). Of these, 31% were women with mean age of 18.27 years (SD = 0.56). All participants were young drivers with an average of driving experience of 12 months since receipt of their driver's license. About 73% of them drove frequently during the week, 26% drove once a week. The survey consisted of 24 items from three self-reported questionnaires: (1) the Reckless Driving Behavior Scale (RDB) (McNally and Bradley 2014) used to assess those behaviors that increase the probability of a vehicle crash due to driving under the influence of

3
Modeling random and non-random decision uncertainty in ratings… substances (drugs), extreme motorsport behaviors (extreme), and speeding/steering behaviors (positioning); (2) a short version of the Driving Anger Scale (DAS) (Deffenbacher et al. 1994), adopted to assess driving angers provoked by someone else's behaviors like slow driving and discourtesy; (3) a simplified version of the Family Climate Road Safety (FCRS) questionnaire (Taubman-Ben-Ari and Katz-Ben-Ami 2013), adopted to evaluate the role of parents in teens' safe driving, especially with regards to communication, monitoring, and parents' messages. Questionnaires were administered using DYFRAT , a computerized fuzzy rating scale which adopts the mouse-tracking methodology (Freeman and Ambady 2010) as a tool to implicitly quantify rating processes. For each item of the survey, participants were asked to respond using a pseudo-circular rating scale with six-points (RDB), five-points (DAS), and four-points (FCRS) anchors, respectively. Participants gave their responses by mouse-clicking the chosen level of the scale. Meanwhile, we recorded the streaming x − y coordinates of the computer-mouse needed to reach the chosen anchor. According to the DYFRAT methodology, the linearized histograms of the collected radians were used to compute beta fuzzy sets for each item and for each participant (fuzzy sets were obtained by means of an heuristic optimization procedure which converts histograms into fuzzy sets). Finally, modes 1 , … , n and precisions 1 , … , n were used to represent final responses and uncertainties involved during the rating process. Figure 1 shows an example of observed beta fuzzy numbers for the RDB response variable. Data analysis and results The fuzzy response variable rdb was computed by aggregating the fuzzy variables extreme and positioning of the RDB questionnaire in terms of mean (Hanss 2005). Similarly, fcrs was obtained by aggregating the crisp responses for the variables monitoring, messages, and communications from the FCRS questionnaire. To simplify interpretation of the results, variables drugs and fcrs were made categorical using median split on their crisp responses. This yielded to two new dichotomous variables, namely drugs (non-use/use) and fcrs (bad/satisfactory family climate). By contrast, the variable anger entered the model as a fuzzy variable in terms of mode and precision components. First, we run a model (model 1) where dispersion was held fixed for all participants whereas the mean was modeled using the fuzzy variable anger and the categorical variables fcrs and drugs. shows the final estimates along with their standard errors computed using the empirical observed information matrix (see  and Pearson's residuals. The results suggest that rdb increased as a function of mode components of anger ( ̂1 = 1.351 , ̂1 = 0.356 ) whereas its precision/spread components did not affect the response variable ( ̂2 = 0.001 , ̂2 = 0.002 ). Moreover, participants using drugs showed higher levels of rdb ( ̂4 = 0.458 , ̂4 = 0.187 ) when compared to those who did not use substances. As expected, participants with satisfactory family climate showed lower levels of rdb ( ̂3 = −0.192 , ̂3 = 0.440 ) when compared to participants with a bad family climate. To account for heterogeneity in the response variable rdb, two further models were estimated, one which included the dichotomous variable sex, and the other which also included driving_frequency. In order to evaluate models improvements, both the models were compared in terms of fuzzy likelihood-ratio test (Berkachy and Donzé 2019). Table 4 reports the results for model 2 and model 3. With regards to model 2, the likelihood-ratio test computed against model 1 reveled that sex improved the fit of the model ( 2 7−6=1 = 5.001 , p = 0.025 , AIC model1 = −163.38 , AIC model2 = −166.38 ), with men showing lower levels of heterogeneity in the response variable ( 1 = −1.199 , ̂1 = 0.613 ) when compared to women. To further analyse the variability in reckless-driving behaviors, we asked whether this varied as a function of participants' frequency of driving. To assess this hypothesis, model 3 included driving_frequency as an additional term in the precision equation. The likelihood-ratio test conducted against model 2 revealed that driving frequency did not improve the fit of the model ( 2 8−7=1 = 1.626 , p = 0.2021 , AIC model2 = −166.38 , AIC model3 = −166.01 ). The results were in line with the literature (e.g., see: Bıçaksız and Özkan 2016) and suggested that self-reported reckless-driving behaviors were positively associated to driving anger ( ̂1 = 1.485 , ̂1 = 0.378 ) and substance use ( ̂4 = 0.518 , ̂4 = 0.235 ). By contrast, family climate acted as a protective factor with satisfactory family climate being negatively associated to risky behaviors ( ̂3 = −0.126 , ̂3 = 0.157 ). Interestingly, variability of selfreported responses varied as a function of gender, with female drivers showing more variable responses than male drivers ( ̂1 = −1.999 , ̂1 = 0.613 ). Figure 2 plots the predicted curves against the observed fuzzy data as a function of both continuous (panels A-B) and categorical (colors in both panels) predictors (note that only estimated modes ̂ were plotted for the sake of simplicity). To further investigate the role of fuzziness to predict self-reported reckless-driving behaviors (RDB), we contrasted the results of the final model (model 2) against the same model adapted on mean-based defuzzified response data (i.e., data without the components 1 , … , n for the decision uncertainty). Table 5 shows the final estimates computed using maximum-likelihood as implemented in standard beta regression model (Ferrari and Cribari-Neto 2004). The coefficients ̂ for the mean ̂ component were nearly similar to those obtained for the fuzzy case. By contrast, the coefficients ̂ for the precisions were lower if compared to the fuzzy case. This is in line with the results of the simulation study. However, standard errors of the estimates for defuzzified data were smaller than those obtained for the fuzzy case, especially with regards to precision coefficients ̂ . Overall, this is coherent withthe estimation approach used for fuzzy data as, in this case, standard errors were computed using score Eqs. (16)-(17) which are in turn based on filtered data ( * 1 , … , * n ).  Note that categorical variables were codified as for the fuzzy beta linear models (see Table 4). The parameters and were linked to the response variable using logit(⋅) and log ( (a, b). Fitted curves correspond to model 2 in Table 4. Note that rectangles represent -cuts of the observed fuzzy data with = 0.5 , i.e.

3
Modeling random and non-random decision uncertainty in ratings…

Case study 2: Service quality in restaurant industry
Service quality is an important factor to assess the performances of a company and it is of a great interest for restaurants services as well. Several research have been conducted to understand the overall effect of perceived service quality on customer satisfaction which, in turn, leads to positive consumption behaviors like revisiting and recommending the restaurant (Almohaimmeed 2017; Namkung and Jang 2008). A number of variables influencing dining experience have been suggested, such as food quality, menu variety, food presentation, quality of staff service, internal/external environments (Ha and Jang 2010). Last but not least, variables like prices and restaurant type (e.g., fast-food, fine restaurants) seem to be also involved in the heterogeneity of restaurant quality (Almohaimmeed 2017). In this application, we will consider a simple model for restaurant service quality where perceived quality of food (food) and staff's perception of being courteous (employees) were used to predict perceived service quality (service_quality). In addition, we also evaluated the extent to which heterogeneity in the response variable can be accounted by price levels (prices) and restaurant type (type). Data and measures Data were originally collected by de Sáa et al. (2014) and refers to a survey of 14 items administered to a sample of n = 70 customers of different age, background, and occupation. The questionnaire included two of the most important factors of restaurant quality, namely food/beverage and service quality. We considered only complete cases, i.e. cases with no missing values. This yielded to a subset of n = 49 customers (31% women) with modal age between 25 and 34 years. Informal restaurants were about 67% of the total, 16% of them were fine, 10% fast-food, and 7% were self-service restaurants. About 69% of restaurants showed prices levels on average (about 15 Euro) whereas 31% of them reported higher prices. Ratings measures were collected by the authors using a Likert-type computerized fuzzy rating scale (de Sáa et al. 2014). In this case, participants provided their responses by means of a two-step direct rating procedure. First, they were asked to draw on the graphical pseudo-continuous scale the core of the set containing the most plausible rating values. Then, conditioned on the previous choice, they were asked to draw the support of the set, which instead represents the most compatible set for the rating values. Finally, both the intervals were interpolated to form a trapezoidal fuzzy response. For the purposes of this application, trapezoidal fuzzy numbers were converted into beta fuzzy numbers adopting a procedure minimizing the information content of the original fuzzy sets. In general, several transformation procedures are available to this purpose (Nasibov and Peker 2008). Here, for the sake of simplicity, we resorted to adopt the simplest approximation based on the minimization of the area under the curve between two fuzzy sets. In particular, the parameters {m, s} ∈ (0, 1) × ℝ + of the beta fuzzy set B (y;m, s) that approximates a trapezoidal fuzzy set Ã (y;a, b, c, d) with real parameters a < b < c < d were found by minimizing the function (m, s) = | ∫ A 0 Ã (y;a, b, c, d)dy − ∫ B 0 B (y;m, s)dy| w.r.t. m and s. Similarly to the first application, modes 1 , … , n of the beta fuzzy numbers represent final participants' responses whereas their precisions 1 , … , n model the decision uncertainty occurred during the rating process. Figure 3 shows an example of transformed beta fuzzy numbers along with the associated observed trapezoidal fuzzy sets for the service_quality response variable.
Data analysis and results To take into account decision uncertainty for predictors, in this case fuzzy variables food and employees were defuzzified using the centroid method before entering the model. A first model was run including food and employees as predictors for the mean and type and price for precision components. Table 6 shows the final estimates along with their standard errors. As expected, service_quality increased as a function of food ( ̂1 = 1.870 , ̂1 = 0.856 ) and employees ( ̂2 = 1.501 , ̂2 = 0.718 ). This is in line with previous studies on restaurant quality, suggesting that a higher perceived quality of food and staff predicts the perceived quality of restaurant services. In addition, the dispersion component of the model varied as a function of type, with informal restaurants showing higher heterogeneity in service_quality ( 1 = 1.084 , ̂1 = 1.055 ) then self-service restaurants. The same applied for fine restaurants ( 3 = 0.865 , ̂3 = 1.158 ) whereas fast-food showed decreasing levels of variability in service_quality ( 2 = −0.232 , ̂2 = 1.124 ) when compared to self-service restaurants. Interestingly, variability in service_quality is positively associated to price ( 4 = 1.802 , ̂4 = 1.239 ), with high-priced restaurants being more homogeneous in terms of perceived quality. This indicates that service_quality was not homogeneous over normal-priced restaurants and further covariates like internal/external atmospherics or image of restaurant might instead be needed to further account for such differences (Almohaimmeed 2017; Ha and Jang 2010). Figure 4 plots the predicted curves against the observed fuzzy data as a function of both continuous predictors (note that only estimated modes ̂ were plotted for the sake of simplicity). As for the previous case study, we compared the fuzzy beta model with the beta regression model adapted on mean-based defuzzified data. Table 7 shows the final estimates computed using maximum-likelihood as implemented in standard beta regression model (Ferrari and Cribari-Neto 2004). In line with the results from the simulation study, the estimated coefficients ̂ for the precision component ̂ were lower than those for the fuzzy case and showed smaller standard errors. As stated previously, this can be interpreted in light of the estimation method used for the fuzzy case.   Note that categorical variables were codified as for the fuzzy beta linear models (see Table 6). The parameters and were linked to the response variable using logit(⋅) and log(⋅) link functions, respectively

Conclusions
In this article we developed a statistical approach to deal with bounded continuous ratings data in the case of non-random uncertainty. In particular, beta fuzzy numbers were adopted to represent ratings data subject to decision uncertainty and the beta regression framework was used to model the random counterpart of the overall rating process. The fuzzy component of the data was then used to estimate the latent and non-fuzzy characteristics of the raters population. Parameters estimation was performed by maximum likelihood using a version of the Expectation-Maximization algorithm generalized for the case of fuzzy data. A simulation study and two real applications were used to highlight the characteristics of the proposed approach. The simulation study revealed that the fuzzy beta linear model showed more accurate results over a set of standard methods which can be applied in the case of fuzzy data. The applications showed how the proposed method can be adopted in real cases involving ratings data represented as fuzzy numbers. A nice advantage of the proposed approach is its simplicity and flexibility in dealing with fuzzy data. Indeed, as it encapsulates decision rating uncertainty directly in ̃ , the fuzzy beta linear model does not require the extension of its parametric structure to account for modes and precisions of beta fuzzy data. Indeed, in the current proposal, the model's parameters = { , } are not represented as fuzzy numbers. Consequently, parameters estimation and inference can still be performed using the asymptotic properties of maximum likelihood theory. In this setting, the fuzzy beta model recovers the parameters of the underlying rating process Y ∼ f Y (y; ) by filtering the imprecise data in terms of expectation Y|ỹ . However, although this constitutes an important advance, it should be noticed that fuzziness of the data is Case study 2: Observed fuzzy data for service_quality as a function of the two continuous predictors food (a) and employees (b). Fitted curves correspond to model 1 in Table 6. Note that rectangles represent -cuts of the observed fuzzy data with = 0.5 , i.e. i = min {y ∈ [0, 1] ∶ � y i (y) > 0.5} , max {y ∈ [0, 1] ∶ � y i (y) > 0.5} not propagated to the output of the model. Indeed, as a consequence of the averaging approach upon which the fuzzy-EM is based, the predictions ̂ of the fuzzy beta model are non-fuzzy and they are made at the level of the underlying crisp rating mechanism f Y (y; ) . This may limit the use of this approach in some circumstances, for example when researchers are interested in forecasting fuzziness based on current fuzzy data, or rather when components of fuzzy data play a different role in predicting the variable response (e.g., modes and precisions/spreads of the outcome variable interact with those of the explanatory variables in some way). In all these cases, different statistical approaches may instead be preferred, like those based on full-likelihood approaches for imprecise data (e.g., see Denoeux 2014;Kanjanatarakul et al. 2016) or min-max based estimation methods (e.g., see Guillaume and Dubois 2020). Various possible extensions of our approach can be considered in future works. For instance, fuzzy beta model involving fuzzy data with different shapes simultaneously (e.g., beta, triangular, trapezoidal) would offer a way to deal with more complex scenarios. Another future generalization which might be interesting to investigate is the case where fuzziness in ratings data is coupled with random uncertainty which varies as a function of subgroups in the model (like for CUBE models, see . This would offer the opportunity to further decompose the overall uncertainty in ratings responses in terms of different and possibly interacting components underlying participants' rating processes. Finally, an attractive extension of the current approach would consider the case of multivariate fuzzy beta models where joint fuzzy sets do not obey to the product rule (see Eq. 4). In this case, a fuzzy copula representation may instead be used to formally represent the joint fuzziness information (e.g., see Ranjbar and Hesamian 2017).