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 nonrandom 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 [1,55]. 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 [69]. 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 [11,31]. 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 [80]. In this sense, it reflects the subjective interplay of decisional and emotional components which contribute to the final rating response [44]. 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 [68]). Indeed, it is widely recognized that ratings data often suffer from lack of accuracy, for instance because of social desirability [32], faking behaviors [48], personality [56], response styles [27], and violations of rating rules [43,64,66]. 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 [5,21,84]), metrology (e.g., see [60,61]), and risk analysis (e.g., see [72]).
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) [50], beta regression models [29,54,59], and combination of uniform and shifted binomial (CUBE) models [35,62,63]. These models typically represent mean and dispersion components as a linear or nonlinear 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 nonrandom 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 [28] 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 [19,45]. In this setting, maximum likelihood estimation and inference are carried out through the Expectation-Maximization algorithm adapted for the case of fuzzy data [24,74]. 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 [2,13,85]. Similarly, despite the fact that many formal theories have been proposed to deal with subjective uncertainty (e.g., soft sets, rough sets [46,47]), 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 [10,15,19,33,36]).
The reminder 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, Section 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/antcalcagni/fuzzyratingbeta.

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 setÃ of a universal set A is defined by means of its characteristic function ξ A : 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Ã c = {x ∈ A : ξ A (x) = max y∈A ξ A (y)}. In the case max x∈A ξ A (x) = 1 thenÃ is a normal fuzzy set. IfÃ is a normal and convex subset of R thenÃ is a fuzzy number [10]. Broadly speaking, fuzzy sets can be conceived as subsets of R where their Boolean characteristic function ξ A (x), ∀x ∈ A ⊂ R, has been generalized to the real interval [0, 1]. The class of all normal fuzzy numbers is denoted by F (R). Fuzzy numbers can conveniently be represented using parametric models, through which ξ A 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 [10]. Relevant classes of parametric fuzzy numbers are the so-called LR-fuzzy numbers [25] and their generalizations [12]. 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 [16]. A broader class that encompasses a wide range of fuzzy numbers is the so-called beta fuzzy number [3,6,73]: 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 ∈ R and spread/precision s ∈ R + 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 ξ A is still a normal fuzzy set: Interestingly, the re-parameterized Beta fuzzy number resembles the shape of Beta density distribution written using the PERT representation [82].
When a probability space is defined over the reals, the probability of a fuzzy set P(Ã) 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 [83], conditional probability of prior information [18,71], imprecise probability [26], fuzzy numbers [40], and likelihood induced by random events [14]. Following the findings of [24], in this contribution we adopt Zadeh's definition of fuzzy probability [83]. In particular, let (R, B, P) be a probability space. Then, P(Ã) is defined as follows: with ξ A being Borel measurable. In this context, two fuzzy setsÃ andB are said independent w.r.t. to P, if P(ÃB) = P(Ã) · P(B), with the fuzzy product being defined as ξÃB(x) = ξÃ · ξB [83]. The conditional probability of two independent fuzzy events is with P(B) > 0. Note that one can also obtain the conditional probability between a crisp set A and a fuzzy setB 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 = (y 1 , . . . , y n ) from Y 1 , . . . , Y n is available, the likelihood of the sample y can be generalized as follows: where the definition ξ y = n i=1 ξ yi has been used for the joint fuzzy set [33]. Further details about fuzzy generalization of likelihood functions, fuzzy random variables, and fuzzy probability space can be found in [14,20,24,33,34].
We interpret fuzzy data in the context of random variables following the epistemic viewpoint on fuzzy set theory [19]. In particular, for a fuzzy setÃ, ξ A (Y = y) is interpreted as the possibility that the crisp event Y = y has to occur. Indeed, ξ A (Y = y) ∈ (0, 1) can be conceived as a graded plausibility about the occurrence of the event Y = y, with ξ A (Y = y) = 1 indicating the fact that Y = y is fully possible. By contrast, ξ A (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 ξ A (Y = y) is thought as being the consequence of a two-step generation process, in which first a realization y is drawn from Y and then a fuzzy set ξ A 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 ratings data
Since the seminal work of [41], fuzzy sets have been extensively used in the context of ratings data (for a review, see [11,49]). 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 [44]. 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. While the first one aims at turning crisp ratings into fuzzy data by means of supervised or unsupervised conversion systems (e.g., see [81]), fuzzy rating tools instead offer an embedded interface with which the unfolding of the rating process is explicitly or implicitly quantified (e.g., see [11,22]). In the latter case, for instance, computerized tracing techniques can be used to measure the underlying rating process and fuzzy sets can be naturally adopted to represent both final rating responses (e.g., using the coreÃ c of the sets) and decision uncertainty (e.g., using some of the α-sets ofÃ).
In this context, beta fuzzy numbers can be linked to the rating process Y ∼ f Y (y; θ) 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 = y) = 1, s is the precision of m such that smaller values indicate larger levels of hesitation in the rating choice, and ξ Y 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 = (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 [29]. 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 R, as follows: where X and Z 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) [50]. 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 t = ln y, ln (1 − y) being 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 y ∼ f Y (y; µ, φ) has been realized. This leads to a situation where the sample y cannot be precisely observed and a collection of fuzzy dataỹ is instead available. When fuzzy data are represented as beta fuzzy numbers, theñ with m and s being n×1 vectors of modes and precisions/spreads for the fuzzy observations. Turning Eq. (3) into (6), the joint density ofỹ can be written as: where the joint fuzzy set ξ y = n i=1 ξ yi 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 m and s of the fuzzy sets enter the model as observed quantities whereas the parameters β and γ still remain non-fuzzy quantities.

Parameters estimation
To provide estimates for θ = (β, γ) in the context of fuzzy ratings data, one can maximize the log-likelihood l(θ;ỹ) function, which is obtainable 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 [24]. As for the standard EM algorithm, the fuzzy-EM version at the k-th iteration alternates between the E-step, which involves the computation of the expected complete-data log-likelihood 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 [24]). In the fuzzy-EM variant, the complete-data log-likelihood is that obtained if y was precisely observed (see Eq. 9). Therefore, given the (k − 1)-th estimates θ ′ = θ (k−1) the E-step for the k-th 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 Yi|ỹ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 E θ ′ ln Y i ỹ i can be approximated via Taylor expansion around E θ ′ Y i ỹ i as follows: Similarly, the second expectation E θ ′ 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 denote the filtered data which are computed using Eqs. (14)- (15). Finally, θ (k) is obtained by computing the roots of Eqs. (16)- (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 Iθ from the Hessian of the maximum likelihood estimatesθ obtained solving Eqs. (16)-(17) for β and γ. In the context of EM algorithm, Iθ can be approximated using the empirical observed information matrix [53], as follows: is the score vector for the i-th observation calculated atβ andγ. The standard errors are then calculated as usual: Note that this approximation avoids the computation of Hessian of the complete-data log-likelihood as it solely uses the score equations where the unobserved vector y is replaced by y * (1) and y * (2) . Alternatively, standard errors can also be obtained via non-parametric bootstrap [51]. 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(x iβ )) −1 andφ = exp(z iγ ) whereas ξ yi (μ 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. Finally, likewise for the non-fuzzy case, inference on β and γ can be performed using maximum-likelihood theory [24,51] and, consequently, hypothesis testing on model's parameters can be performed using fuzzy version of likelihood ratio test (e.g., see [7,57]).

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 [21]), 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,16x4 GB Ram whereas computations and analyses have been done in the R framework for statistical analyses. (a) X n l ×j k = [1 n l , X n l ×j k−1 ] and Z n l ×hm = [1 n l , Z n l ×hm−1 ] were drawn from U nif (1, 5) (b) location terms were computed as follows: µ n l ×1 = logit −1 1 + exp(X n l ×j k β 0 j k ×1 ) (c) precision terms were computed via exp function: φ n l ×1 = exp(Z n l ×hm γ 0 hm×1 ) (d) 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 (e) fuzzy data were generated by making y n l imprecise via a two-step data-generation process [65,75]. First, spread components were generated as s 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 (f) parameters β j k and γ hm 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 ξ yi (y) dy and the first-maximum method y * i = sup y∈(0,1) {ξ yi (y)} for i = 1, . . . , n l . The ML procedure implemented in the R library betareg has been used in this case [85].
Outcome measures. For each condition of the simulation design, sample results were evaluated using bias of estimates and root mean square errors.
Results. Tables 1-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 indirect fuzzy rating scales [11]. By contrast, the second application is about the analysis of customer satisfaction data collected using direct fuzzy rating scales [22].

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 [79]. 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 [8,52,70,76]). Several studies have recognized the role of subjective factors like sensation-seeking, normlessness, anxiety, aggressiveness, driving attitudes in determining risky behaviors [8]. Researchers have also assessed the contribution of parenting styles and peer relations to young drivers' intention to take risks [76,77]. 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 considered 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: (i) the Reckless Driving Behavior Scale (RDB) [52] used to assess those behaviors that increase the probability of a vehicle crash due to driving under the influence of substances (drugs), extreme motorsport behaviors (extreme), and speeding/steering behaviors (positioning); (ii) a short version of the Driving Anger Scale  Table 3: Monte Carlo study: overall ratios r between over and under-estimation, percentage p of over-estimation, and overall root mean square error. Note that all the indices were computed over B = 5000 replications.
(DAS) [23], adopted to assess driving angers provoked by someone else's behaviors like slow driving and discourtesy; (iii) a simplified version of the Family Climate Road Safety (FCRS) questionnaire [78], 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 [11], a computerized fuzzy rating scale which adopts the mouse-tracking methodology [31] as a tool to implicitly quantify rating processes. According to this rating technique, participants' rating responses were represented in terms of beta fuzzy numbers, with modes m 1 , . . . , m n representing final ratings responses and precisions s 1 , . . . , s n the decision uncertainty involved during the rating process [11].
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 [39]. Similarly, fcrs was obtained by aggregating 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. 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. Table 4 shows the final estimates along with their standard errors computed using the empirical observed information matrix (see Eqs. [18][19] 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 [7]. Table 4 reports the results for model 2 and model 3. With regards to model 2, the likelihoodratio 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). Overall, results were in line with the literature (e.g., see [8]) and suggested that self-reported reckless-driving behaviors were positively associated to driving anger and substance use. By contrast, family climate acted as a protective factor with satisfactory family climate being negatively associated to risky behaviors. Interestingly, variability of self-reported responses varied as a function of gender, with female drivers showing more variable responses than male drivers.

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 [4,58]. 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 [38]. 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 [4]. 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 [22] 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-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 fuzzy rating scale with trapezoidal fuzzy numbers [22]. Unlike DYFRAT technique [11], participants answer using a two-stage strategy where they first are asked to draw the core of the set, which contains the most plausible rating value. Then, conditioned on the previous choice, they are asked to draw the support of the set, which instead represents the most compatible set for the rating. For the purposes of this application, trapezoidal fuzzy numbers were converted into beta fuzzy numbers adopting a procedure minimizing the information content of the fuzzy sets [17]. Similarly to the first application, modes m 1 , . . . , m n of the beta fuzzy numbers represent final participants' responses whereas their precisions s 1 , . . . , s n model the decision uncertainty occurred during the rating process.
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 5 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 [4,38].

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 m and precisions s 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 E θ [Y |ỹ]. However, although this constitutes an important advance, it should be noticed that fuzziness of the data is 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 (e.g., see [30,37,42]).
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 [62]). 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 [67] Table 5: Case study 2: Variable dispersion fuzzy beta models for service quality in restaurant industry. Note that categorical variables were codified using dummy coding with the following reference levels: type (ref.: self-service), price (ref.: high). The parameters µ and φ were linked to the response variable using logit(.) and log(.) link functions, respectively.