A mixed beta regression approach for analyzing bounded sensory data in incomplete block designs

ABSTRACT

Sensory analysis plays a crucial role in the food industry, facilitating innovation and enhancing products. This study includes a sample of individuals, both trained and untrained, who evaluate a product using a hedonic scale or scoring system within a balanced incomplete block (BIB) design. In this context, integrating sensory analysis with robust statistical methods that account for the nature of the response variables is essential to the success of the experimental study. Techniques for analyzing sensory data include response surface models and categorical models. This article proposes using beta regression as a viable alternative to the proportional odds model to address convergence issues, particularly given the number of parameters involved. Furthermore, the beta distribution offers flexibility for modeling data with heteroscedasticity and skewness. To support this approach, we conducted simulation studies that demonstrated agreement rates in product selection using both models. Additionally, we present a compelling case study to guide the selection of grape juice formulations. In this application, the findings from the mixed beta regression model reinforced the choices made using the proportional odds mixed model.

Keywords:
likelihood procedure; simulation studies; formulation selection; sensory attributes; random effects

Introduction

Sensory analysis is essential for assessing the acceptability of food products, covering evaluations from raw materials to final formulations. This analysis examines various attributes, including color, flavor, acidity, and aroma, through the senses of sight, touch, smell, and taste (Duan et al., 2024; Kottaridi et al., 2023). Nevertheless, challenges such as panelist fatigue and high costs can arise. These issues can be addressed by employing incomplete block designs, which enhance efficiency and accuracy (Sipos et al., 2021).

In this scenario, hedonic scales are commonly employed in sensory evaluation, producing ordinal data that may be influenced by panelists’ memory biases, potentially leading to variability and errors (Sugumar and Guha, 2022). It is crucial to utilize appropriate statistical methods to reduce these sources of error and enhance the reliability of the conclusions drawn from the data. Several studies in the literature treat these scores as continuous variables and apply analysis of variance, linear regression, or multivariate techniques, assuming that the underlying assumptions are robust (Smithson and Verkuilen, 2006). However, this approach is often incorrect and can result in misleading conclusions. While nonparametric tests, such as Friedman's test, are frequently used, they tend to be less powerful than parametric tests. Conversely, the generalized logit model offers an efficient alternative, although it may encounter challenges such as convergence issues and the necessity for category grouping (Agresti, 2010).

To address these limitations, this study explores the use of beta regression models with random effects for analyzing sensory data. Given their inherent flexibility, beta regression models are well-suited to handling skewed and heteroscedastic data constrained to bounded intervals (Kubinec, 2023). Therefore, these models bypass the assumptions of normality and homogeneity of variance (Junaid et al., 2024), thereby minimizing complications associated with convergence, overparameterization, and the necessity for grouping in logit models. Additionally, we present a simulation study that compares the performance of cumulative logit models and beta regression, followed by an application to real sensory data.

Materials and Methods

The dataset used in this study was obtained from a sensory evaluation conducted in Piracicaba, São Paulo State, Brazil (22°42’30" S, 47°37’30" W, altitude 546 m) in 2025. The experimental design employed a balanced incomplete block (BIB) design, where each of the 98 untrained panelists assessed four of the eight grape juice formulations. The tested products included eight commercial grape juice brands, divided equally into two categories: four whole juices (100 % juice or not from concentrate; formulations F₈₇₃, F₆₆₁, F₄₁₉, and F₅₇₁) and four reconstituted juices (from concentrate; formulations F₇₁₅, F₁₇₉, F₇₃₂, and F₃₁₈).

Each panelist rated the four samples on a five-point hedonic scale, with higher scores indicating more favorable perceptions. The evaluated sensory attributes included acidity, aroma, color, sweetness, and flavor.

Review of the beta regression model with a random effect

The beta regression model is well-suited for response variables that are strictly bounded within the open interval (0, 1), such as rates and proportions. Beta regression models that incorporate random effects can be viewed as extensions of generalized linear mixed models (GLMMs) (Figueroa-Zúñiga et al., 2013). In GLMMs, the linear predictor integrates fixed and random effects to address correlation structures that typically arise from repeated measurements within subjects (Salinas-Ruíz et al., 2023).

Consider a longitudinal study where $y_i={\left(y_{i1},\dots ,y_{i{n}_i}\right)}^{\text{T}}$ is a vector of n_i observations for the i-th individual, with y_ij denoting the j-th response and i = 1, 2, … N. Let x_ij= (x_ij1, …, x_ijp)^T be the (n_i × p) matrix associated with the vector β = (β₁, …, β_p)^T of unknown regression coefficients related to fixed effects, and z_ij = (z_ij1, …, z_ijq)^T be the (n_i × q) matrix related to the random effects vector u_ij = (u_i1, …, u_iq)^T. Assuming Y_ij | u_i are conditionally independent random variables and follow a beta distribution (μ_ij, ϕ); u₁, …, u_q independent random variables with distribution $u_i\mid {\mathrm{\sigma }}_u^2\sim N_q\left(0,{\mathrm{\sigma }}_u^2\text{I}\right)$, the mixed beta regression model is defined as:

(1) $$g\left(E\left(Y_i\mid u_i\right)\right)=x_i^T\text{β}+Z_i^T{u}_i={\text{η}}_i$$

where g(.) is a known, strictly monotonic, and twice differentiable link function. The logit link function is commonly used (though not mandatory) and transforms Eq. (1) into (Zimprich, 2010):

(2) $$\log \left({\text{μ}}_{ij}/1-{\text{μ}}_{ij}\right)={\text{η}}_{ij}=x_{ij}^T\text{β}+Z_{ij}^T{u}_i$$

Accordingly, the conditional mean μ_ij= E (Y_ij | u_i) by means Eq. (2) is expressed as:

$${\text{μ}}_{ij}=\exp \left({\text{η}}_{ij}\right)/\left(1+\exp \left({\text{η}}_{ij}\right)\right)$$

For further methodological details on beta regression models and their extensions, refer to Verkuilen and Smithson (2016), Stroup (2012), Zimprich (2010), Huang and Oosterlee (2011), and Pereira et al. (2013).

Proposed method

In sensory studies structured under a BIB design, v treatments (or formulations) are evaluated across b blocks (panelists), with each panelist assessing k treatments. Each treatment appears in exactly r blocks, and every pair of treatments co-occurs in precisely ƛ blocks. Let L sensory attributes be evaluated using a hedonic scale with c response levels. The BIB design must satisfy the following conditions (Yu et al., 2024): bk = rv.

Number of variety pairs:

$$(v(v-1)/2)$$

Number of treatment pairs within each block:

$$(k(k-1)/2)$$

Number across the entire experience:

$$(bk(k-1)/2)=(\lambda v(v-1)/2)$$

In the last equality, replacing bk with rv, we have:

$$\lambda (v-1)=r(k-1)$$

These relationships facilitate balanced comparisons across treatments and pairwise comparisons. To assess the adequacy and sensitivity of this design under the specified conditions, a simulation-based power analysis was conducted using the SIMR package in R software (version 4.3.2) (Green and MacLeod, 2016). The simulation replicated the experimental setup in which 98 untrained panelists evaluated four out of eight formulations. Sensory scores were simulated under the assumption of a fixed difference of 0.15 between formulations F₈₇₃ and F₇₁₅, while incorporating random variation among panelists. The response variable was constrained to the interval (0, 1) to reflect the bounded nature of hedonic scores.

A mixed-effects model was fitted using the lmer() function from the lme4 package (Bates et al., 2015), with formulation as a fixed effect and panelist as a random intercept. Power analysis based on 100 Monte Carlo simulations yielded a power estimate of 99 % (95 % confidence interval: 94.6-99.97 %) for detecting the specified effect size with the current sample size of 392 observations. The large number of observations, combined with the balanced BIB design, provided sufficient degrees of freedom, confirming that the design possesses adequate sensitivity to identify meaningful differences between formulations (Morris et al., 2006).

To apply the modeling framework discussed in the section "Review of beta regression model with random effect", the original responses y_i, which were collected on a five-point hedonic scale (ranging from 1 to 5), were rescaled by dividing each score by five. To ensure the values fell strictly within the open interval (0, 1), as required by the beta distribution, we applied the transformation proposed by Smithson and Verkuilen (2006):

(3) $$y_i^{\ast }=\left(y_i(n-1)+0.5\right)/n$$

where n = bk. After this transformation, the scores are treated as realizations of a continuous variable with bounded support (Sarzo et al., 2023), in accordance with the beta regression class proposed by Ferrari and Cribari-Neto (2004).

To evaluate the product's influence on a sensory attribute, three models were considered. The null model is:

(4) $${\text{η}}_0=\log \left({\text{μ}}_i/1-{\text{μ}}_i\right)={\text{η}}_0$$

which assumes that the mean score of the sensory attribute is independent of the formulation. The second model includes the formulation effect:

(5) $${\text{η}}_1=\log \left({\text{μ}}_i/1-{\text{μ}}_i\right)={\alpha }_0+X_i^T\text{β}$$

where $x_i={(x_{i{t}_1},\dots ,x_{i{t}_k})}^T$ with 1 ≤ t₁, …, t_k ≤ v, representing the formulation indices in each block (panelist), and β = (β₁,.., β_v)^T is the associated coefficient vector capturing the effect of formulation. In this setup, the dimensions are: the design matrix x (bk × v), the response vector y (bk × 1), and the coefficient vector β (v × 1). The design matrix includes dummy variables for the formulations. The third model includes random effects to account for the incomplete block structure and the uncontrolled panelist variability:

(6) $${\text{η}}_2=\log \left({\text{μ}}_i/1-{\text{μ}}_i\right)={\text{α}}_0+X_i^T\text{β}+u_i$$

where $u_i~N\left(0,{\mathrm{\sigma }}_u^{2}I\right)$ corresponding to a specific case of the general Eq. (2) for a given sensory attribute L. For all models (Eq. 4, 5, and 6), the conditional mean is obtained by:

(7) $${\widehat{\text{μ}}}_i=\exp \left({\text{η}}_i\right)/1+\exp \left({\text{η}}_i\right)$$

Parameter estimation for the mixed beta regression model (Eq. 6) was performed using maximum-likelihood estimation. The likelihood contribution from panelist i is given by:

(8) $$f_i\left(y_i\mid \text{β},{\mathrm{\sigma }}_u^2,\mathrm{\varphi }\right)=\int {\displaystyle \prod }_{j=1}^k{f}_{ij}\left(y_{ij}\mid u_i,\text{β},\mathrm{\varphi }\right)f\left(u_i\mid {\mathrm{\sigma }}_u^2\right)d{u}_i$$

where f_ij(y_ij | u_i, β, ϕ) is the conditional beta density for the observation y_ij (i.e., the j-th evaluation in block i), and $f\left(u_i|{\mathrm{\sigma }}_u^2\right)$ is the density of the random effect $u_i~N\left(0,{\mathrm{\sigma }}_u^{2}I\right)$. The overall likelihood function is:

(9) $$\begin{array}{ll}L\left(\text{β},{\mathrm{\sigma }}_u^2,\mathrm{\varphi }\right)& ={\displaystyle {\prod }_{i=1}^b{f}_i}\left(y_i\mid \text{β},{\mathrm{\sigma }}_u^2,\mathrm{\varphi }\right)\\ & ={\displaystyle {\prod }_{i=1}^b\int }{\displaystyle {\prod }_{j=1}^k{f}_{ij}}\left(y_{ij}\mid u_i,\text{β},\mathrm{\varphi }\right)f\left(u_i\mid {\mathrm{\sigma }}_u^2\right)d{u}_i\end{array}$$

Maximizing Eq. (9) yields estimates of the formulation effects (β), panelist variance (${\mathrm{\sigma }}_u^2$), and the precision parameter (f). Taking the natural logarithm of Eq. (9) yields the log-likelihood function:

(10) $$\ell \left(\text{β},{\mathrm{\sigma }}_u^2,\mathrm{\varphi }\right)={\displaystyle {\sum }_{i=1}^b\log }\int {\displaystyle {\prod }_{j=1}^k{f}_{ij}}\left(y_{ij}\mid u_i,\text{β},\mathrm{\varphi }\right)f\left(u_i\mid {\mathrm{\sigma }}_u^2\right)d{u}_i$$

The primary computational challenge lies in evaluating the b integrals over the q-dimensional random effects u_i, which lack a closed-form solution. To address this issue, the integrals were approximated using the Laplace method. This approach approximates the integral as a weighted sum centered on the integrand's mode. The method expands the integrand around the mode of the random effects through a second-order Taylor series and approximates the integral by computing a Gaussian integral at that mode. It is particularly effective for models with non-Gaussian responses and latent variables (Tierney and Kadane, 1986).

The model was estimated using R software with the glmmTMB package, which is built on the Template Model Builder (TMB) framework. This framework allows for the specification of complex likelihood functions that incorporate random effects and facilitates optimization through the Laplace approximation. The code is available at: https://github.com/GabrielRPalma/MixedBetaRegression/tree/main.

After model fitting, selection was based on the maximized log-likelihood, Akaike Information Criterion (AIC), and the Likelihood Ratio Test (LRT). All comparisons were evaluated at the 5 % significance level. The LRT statistic is defined as:

$$\text{Λ}=-2\left[\log \left({\widehat{\text{θ}}}_0\right)-\log \left({\widehat{\text{θ}}}_1\right)\right]$$

where $\log \left({\widehat{\text{θ}}}_0\right)$ and $\log \left({\widehat{\text{θ}}}_1\right)$ are the maximized log-likelihoods under the reduced and full models, respectively. Under standard conditions, $\mathrm{\Lambda }~{\chi }_{df}^2$, where df is the difference in the number of parameters between the two models. The AIC was computed as:

$$\text{AIC}=-2\log (\widehat{\text{θ}})+2\text{k}$$

where $\log \left(\widehat{\text{θ}}\right)$ is the maximized log-likelihood, and k is the number of estimated parameters in the model.

Simulation studies

A simulation study was conducted to evaluate the performance of the beta regression model with random effects in sensory data analysis. This simulation utilized three formulations, arranged in a balanced complete block design. The objective was not to explore alternative experimental designs but rather to assess the modeling strategy compared with the commonly employed cumulative logit model with proportional odds.

Data were generated using parameters from a fitted mixed cumulative logit model with proportional odds, as defined below:

$${\text{η}}_i=\log \left[{\text{γ}}_i/\left(1-{\text{γ}}_i\right)\right]={\text{α}}_i+{\text{β}}^{\text{T}}x+u_i$$

where α_j is the intercept for the j-th response category corresponding to a specific sensory attribute (such as acidity, color, aroma, sweetness, and flavor), β is the vector of fixed effect coefficients associated with the design x_i, and $u_i~N\left(0,{\mathrm{\sigma }}_u^{2}I\right)$ represents the random effect for the i-th panelists. The parameter values used in the simulations are presented in Table 1. The reference category for the response variable was "1 = dislike extremely", and formulation F₁ served as the reference level.

Table 1
Fixed parameters of the proportional probability model, used to simulate data in sensory analysis considering 13 scenarios.

Scenarios	Parameters (logits) θj = (αj, β[f2], β[f3])	Scenarios	Parameters (logits) θj = (αj, β[f1], β[f2])
F1 = F2 = F3	θ1 = (5, 5, 5)	F3 = F1 < F2	θ1 = (5, 5, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 0, 0)		θ3 = (0, 0, 6)
	θ4 = (0, 0, 0)		θ4 = (0, 0, 0)
F1 = F2 < F3	θ1 = (5, 5, 0)	F3 < F1 = F2	θ1 = (5, 0, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 0, 0)		θ3 = (0, 6, 6)
	θ4 = (0, 0, 7)		θ4 = (0, 0, 7)
F1 < F2 = F3	θ1 = (5, 0, 0)	F3 < F1 < F2	θ1 = (5, 0, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 6, 6)		θ3 = (0, 6, 0)
	θ4 = (0, 0, 0)		θ4 = (0, 0, 7)
F1 < F2 = F3	θ1 = (5, 0, 0)	F3 < F1 < F2	θ1 = (5, 0, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 6, 6)		θ3 = (0, 6, 0)
	θ4 = (0, 0, 0)		θ4 = (0, 0, 7)
	θj = (αj, β[f2], β[f3])		θj = (αj, β[f3], β[f2])
F1 < F2 < F3	θ1 = (5, 0, 0)	F1 < F3 < F2	θ1 = (5, 0, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 6, 0)		θ3 = (0, 6, 0)
	θ4 = (0, 0, 7)		θ4 = (0, 0, 7)
	θj = (αj, β[f1], β[f3])		θj = (αj, β[f3], β[f1])
F2 < F1 = F3	θ1 = (5, 0, 0)	F2 = F3 < F1	θ1 = (5, 5, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 6, 0)		θ3 = (0, 0, 6)
	θ4 = (0, 0, 0)		θ4 = (0, 0, 0)
F2 < F1 < F3	θ1 = (5, 0, 0)	F2 < F3 < F1	θ1 = (5, 0, 0)
	θ2 = (0, 0, 0)		θ2 = (0, 0, 0)
	θ3 = (0, 6, 0)		θ3 = (0, 6, 0)
	θ4 = (0, 0, 7)		θ4 = (0, 6, 0)
	θj = (αj, β[f2], β[f1])
F3 < F2 < F1	θ1 = (5, 0, 0)
	θ2 = (0, 0, 0)
	θ3 = (0, 6, 0)
	θ4 = (0, 0, 7)

After generating the data, both the cumulative logit models with proportional odds and the beta regression model with random effects were fitted. The proportional odds logit models (Eq. 11) were estimated using the clmm2function from the ordinalpackage (Christensen, 2015), whereas the beta regression model with random effects (Eq. 6) was implemented using the glmmTMBpackage (Brooks et al., 2017), as detailed in before.

For each scenario, 1,000 simulated datasets were generated for two sample sizes (N = 90 and N = 300 panelists). Agreement rates between the two fitted models are reported in Table 2, and the simulation scenarios are summarized in Table 1.

Table 2
Agreement rates between the cumulative logit model with proportional odds and the beta regression model with random effects for 13 scenarios in the simulation study, considering three formulations (F₁, F₂, and F₃), for N = 90 and N = 300.

Scenarios	Agreement rates
	N = 90	N = 300
	---------------- % ----------------------
F1 = F2 = F3	98.2	87.1
F1 = F2 < F3	99.1	99.4
F1 < F2 = F3	91	94.4
F1 < F2 < F3	99.9	100
F2 < F1 = F3	91.6	92.6
F3 = F1 < F2	100	98.8
F3 < F1 = F2	92.3	93.3
F3 < F1 < F2	99.9	100
F1 < F3 < F2	99.9	100
F2 = F3 < F1	99.7	98.2
F2 < F3 < F1	99.7	100
F3 < F2 < F1	99.9	100

The results demonstrate that, for N = 90, the agreement rates between the models defined by Eq. (11) and (6) were at least 91 % across all scenarios. These rates showed a steady increase with N = 300, although a few exceptions were noted (e.g., F₂ = F₃ < F₁, F₁ = F₂ = F₃, and F₃ = F₁< F₂), likely due to stochastic variation and random effects.

Certain configurations, specifically F₁ < F₃ < F₂, yielded agreement rates as high as 99.9 % for N = 90 and 100 % for N = 300. This pattern was consistently observed across several scenarios, indicating strong agreement among the models regardless of sample size.

Overall, the simulation results show that the beta regression model with random effects aligns closely with the cumulative logit model, making it well-suited for analyzing sensory data collected on bounded ordinal scales.

Results

In this study, the original ordinal scale was converted to a continuous scale in the interval (0, 1) using Eq. (3), enabling the use of beta regression models. Additionally, an exploratory analysis was conducted to examine the distribution of responses across sensory attributes.

Boxplots illustrating the scaled evaluations (ranging from 0 to 1) for each sensory attribute are shown in Figure 1. Most formulations exhibited asymmetric response distributions across all attributes, except for F₇₁₅ and F₇₃₂ for acidity and flavor, and F₇₁₅ for sweetness and F₇₃₂ for aroma, which showed approximately symmetrical patterns. Outliers were detected in several formulations, notably for color (F₁₇₉, F₃₁₈, F₄₁₉, F₅₇₁), sweetness (F₁₇₉), and aroma (F₃₁₈ and F₆₆₁).

Figure 1
Boxplots of the sensory attributes for the eight juice formulations, based on a sensory study conducted in Piracicaba, São Paulo State, Brazil, in 2025.

Median values were equal to or greater than 0.8 for most formulations. Exceptions included F₇₁₅ and F₇₃₂ for acidity; F₇₁₅, F₇₃₂, and F₃₁₈ for flavor; F₇₁₅ for sweetness; F₃₁₈ and F₇₃₂ for aroma, each with a median of approximately 0.6, and F₇₃₂ for color, with a median of 0.4.

This exploratory analysis offers an initial overview of the data distribution. As outlined in the methodology, beta regression models were applied to all sensory attributes. Three model structures were evaluated: the null model (Eq. 4), which includes only the intercept; Model 1 (Eq. 5), which accounted for the fixed effect of formulation; and Model 2 (Eq. 6), which incorporates both formulation as a fixed effect and panelist as a random effect. To determine the most suitable model, we assessed the estimated dispersion parameter ($\widehat{\mathrm{\varphi }}$), the logarithm of the maximized likelihood function (log(L)), and the AIC. A summary of these metrics is provided in Table 3.

Table 3
Estimated dispersion parameter $\widehat{\mathrm{\varphi }}$, the logarithm of maximized likelihood function (log(L)) and Akaike Information Criterion (AIC) for all sensory attributes, considering the three models: Null (Eq. 4), formulation effect (Eq. 5) and with formulation and random effect (Eq. 6) referring to the study carried out in Piracicaba, São Paulo state, Brazil, in 2025.

Sensory attribute	Models	$\widehat{\mathrm{\varphi }}$	log(L)	AIC
Acidity	Null (intercept)	1.7	190.519	–377.0
	1 (intercept + formulation effect)	1.8	3.6	–369.3
	2 (intercept + formulation effect + random effects)	2.0	199.1	–378.2
Aroma	Null (intercept)	2.0	185.1	–366.2
	1 (intercept + formulation effect)	2.2	203.2	–388.4
	2 (intercept + formulation effect + random effects)	2.5	207.8	–395.7
Sweetness	Null (intercept)	1.7	193.3	–382.7
	1 (intercept + formulation effect)	1.8	200.4	–382.8
	2 (intercept + formulation effect + random effects)	2.0	203.0	–386.1
Flavor	Null (intercept)	1.7	200.0	–396.1
	1 (intercept + formulation effect)	1.8	209.0	–400.0
	2 (intercept + formulation effect + random effects)	2.0	211.1	–402.3
Color	Null (intercept)	1.7	473.2	–870.3
	1 (intercept + formulation effect)	2.0	465.2	–912.4
	2 (intercept + formulation effect + random effects)	2.1	265.9	–911.7

The log-likelihood, estimated dispersion parameter, and goodness-of-fit statistics presented in Table 3 progressively increased from the null model to Model 2. Since these metrics indicate the quality of model fit, higher values signify a better fit. Incorporating a random effect for the panelist as an incomplete block proved relevant for elucidating response heterogeneity, as evidenced by lower AIC values for the mixed model across most attributes. Furthermore, the likelihood ratio test for the inclusion of the random effect (p < 0.01) confirmed its significance for all attributes except color, supporting the models described by Eq. (6) as the more suitable functional structure for the data. The estimated variances of the panelist random effect for each attribute were as follows:

$$\begin{array}{ll}{\widehat{\mathrm{\sigma }}}_{acidity}^2=0.20,{\widehat{\mathrm{\sigma }}}_{aroma}^2=0.15,{\widehat{\mathrm{\sigma }}}_{sweetness}^2=0.13,{\widehat{\mathrm{\sigma }}}_{flavor}^2=0.12,\hfill & \hfill \\ {\widehat{\mathrm{\sigma }}}_{color}^2=0.06\hfill & \hfill \end{array}$$

The magnitude of the variance of the panelist random effects for nearly all attributes not only supports the structure of the chosen design but also demonstrates the mixed beta model's sensitivity in identifying variability in responses influenced by unmeasurable factors represented by these random effects.

For the "color" attribute, the model with only the fixed effect had a slightly lower AIC than the mixed model.

The marginal predicted means obtained from the fitted model are summarized in Table 4. Formulations F₄₁₉, F₁₇₉, and F₅₇₁ achieved the highest overall sensory scores across various attributes. Notably, formulation F₄₁₉ consistently recorded the highest predicted values for aroma (0.8375), flavor (0.7995), and sweetness (0.7833), suggesting greater consumer acceptance. In contrast, formulations F₇₁₅ and F₇₃₂ had the lowest estimated means across multiple attributes, particularly for aroma (0.6450 and 0.6451, respectively) and color (0.7808 and 0.5761, respectively). These findings align with the mixed cumulative logit models; however, for brevity, these details are not included in this study.

Table 4
Means estimated values and standard deviations (between parentheses) by Eq. (7) according to the sensory attributes and formulations, resulting from the study carried out in Piracicaba, São Paulo State, Brazil, in 2025.

Formulation	Acidity	Aroma	Sweetness	Flavor	Color
F 179	0.7650 (0.04)	0.8137 (0.03)	0.7424 (0.03)	0.7711 (0.02)	0.8682 (0.01)
F 318	0.7504 (0.05)	0.6839 (0.05)	0.7256 (0.04)	0.7118 (0.04)	0.8488 (0.01)
F 419	0.7592 (0.05)	0.8375 (0.03)	0.7833 (0.03)	0.7995 (0.03)	0.8533 (0.01)
F 571	0.7358 (0.06)	0.7684 (0.04)	0.7479 (0.03)	0.7312 (0.03)	0.8402 (0.02)
F 661	0.7198 (0.05)	0.7492 (0.04)	0.7157 (0.04)	0.7365 (0.03)	0.8136 (0.01)
F 715	0.6701 (0.06)	0.6450 (0.05)	0.6215 (0.05)	0.6323 (0.04)	0.7808 (0.01)
F 732	0.7168 (0.05)	0.6451(0.04)	0.7406 (0.03)	0.6766 (0.03)	0.5761 (0.02)
F 873	0.7193 (0.05)	0.7484 (0.03)	0.7537 (0.03)	0.7577 (0.03)	0.8195 (0.01)

Discussion

Most sensory tests used for product evaluation produce ordinal categorical data from hedonic scales. Cumulative logit models are the standard method for this type of data (Agresti, 2010), and their applications in sensory science are well-documented (Gadrich et al., 2022; Fatoretto et al., 2018). Despite their theoretical robustness, these models often suffer from overparameterization and convergence issues, limiting their applicability to complex experimental designs.

Alternatively, beta regression models with random effects are particularly well-suited to analyzing hedonic scores transformed to the interval (0, 1). Our proposal approach addresses data asymmetry and heteroscedasticity while incorporating evaluator-specific variability through random intercepts. Moreover, this modeling approach allows prioritizing of formulations with desirable profiles, even under ingredient constraints. In this context, the approach offers methodological and statistical advantages and benefits researchers by improving product selection.

Simulation results show that beta regression produces estimates that align with those generated by the cumulative logit model. However, it offers superior numerical stability, better convergence behavior, and greater interpretability of both parameters and results. These advantages underscore its potential as a reliable analytical framework for sensometric applications.

Applying the model to real data from a grape juice sensory study demonstrated its effectiveness in capturing both formulation effects and interindividual variation. Formulations with higher juice content (F₄₁₉ and F₅₇₁) received the highest hedonic scores. However, F₁₇₉, despite its lower juice concentration, also received favorable evaluations. These findings suggest that certain sensory attributes, such as flavor balance and sweetness, can offset compositional shortcomings. This aligns with the observations of Rossiter et al. (2000) for Actinidia deliciosa (A. Chev.) C.F. Liang & A.R. Ferguson, in which higher °Brix levels reduced perceived acidity and enhanced consumer preference.

From an applied perspective, these findings underscore that achieving high sensory acceptance can be attained not only through compositional enrichment but also through targeted sensory optimization. The proposed modeling approach enables prioritizing formulations with desirable profiles, even under ingredient constraints. Additionally, its integration with BIB designs has proven advantageous in scenarios where panelists have limited evaluation capacity. BIB designs yield conclusions comparable to those from complete block designs while reducing experimental time by more than 40 % and costs in consumer research (Silva et al., 2014).

Despite its strength, certain limitations persist. The current model assumes that sensory attributes are independent, a premise that may not always hold across contexts. Furthermore, while the logit link is commonly used in beta models, exploring alternative link functions such as probit, complementary log-log, or Cauchy could yield a better fit, depending on the data distribution, and warrants further investigation.

Another limitation is the lack of formal residual diagnostics. Although the model has shown strong performance in both simulations and empirical applications, diagnostic tools are crucial for verifying model assumptions and identifying outliers (Espinheira et al., 2017). There remains a scarcity of specific techniques for beta models with random effects; for instance, the quantile residuals proposed by Pereira (2019) are limited to models without random effects. To enhance diagnostic assessment, it is advisable to utilize tools such as the quantile residuals from the DHARMa package (Hartig, 2024) or half-normal plots with simulated envelopes (Moral et al., 2017).

Future research should explore multivariate modeling techniques, such as Dirichlet regression and Bayesian hierarchical models, to address attribute correlations and improve predictive accuracy. Additionally, there is a need for studies focusing on residual analysis and diagnostics, particularly in the context of random effects. A promising avenue for expanding the application of this approach in applied sensory science is the simulation-based optimization of BIB designs, which can be explored across various panel sizes, treatment structures, and formulation subsets. Furthermore, a unified analysis of multiple sensory attributes could help optimize outcomes and improve interpretability.

• Declaration of use of AI technologies
  The authors declare that they did not use artificial intelligence technologies for text generation, data analysis, or any other stage of the development of this manuscript. All content was produced independently by the authors.

Data availability statement

The data in this study contains sensitive and confidential information, protected by legal regulations; therefore, it cannot be shared publicly.

Acknowledgments

This publication has emanated from research conducted with the financial support of Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) process number (88887.821274/2023-00) and Science Foundation Ireland under Grant number 18/CRT/6049. Special thanks to the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), process number 300155/2025-5

References

• Agresti A. 2010. Analysis of Ordinal Categorical Data. John Wiley, Hoboken, NJ, USA. https://doi.org/10.1002/9780470594001
  » https://doi.org/10.1002/9780470594001
• Bates D, Mächler M, Bolker B, Walker S. 2015. Fitting linear mixed-effects models using lme4. Journal of Statistical Software 67: 1-48. https://doi.org/10.18637/jss.v067.i01
  » https://doi.org/10.18637/jss.v067.i01
• Brooks ME, Kristensen K, van Benthem KJ, Magnusson A, Berg CW, Nielsen A, et al. 2017. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal 9: 378-400. https://doi.org/10.32614/RJ-2017-066
  » https://doi.org/10.32614/RJ-2017-066
• Christensen RHB. 2015. Ordinal: regression models for ordinal data. R package version 2015.6-28. Available at: https://CRAN.R-project.org/package=ordinal [Accessed Jan 20, 2025]
  » https://CRAN.R-project.org/package=ordinal
• Duan P, Feng X, Li A, Chen H, Peng X, Wang H, et al. 2024. Identification and characterization of aroma profiles of Huajiao oil from different geographical origins using instruments and sensory analysis. Journal of Food Composition and Analysis 129: 106074. https://doi.org/10.1016/j.jfca.2024.106074
  » https://doi.org/10.1016/j.jfca.2024.106074
• Espinheira PL, Santos EG, Cribari-Neto F. 2017. On nonlinear beta regression residuals. Biometrical Journal 59: 445-461. https://doi.org/10.1002/bimj.201600136
  » https://doi.org/10.1002/bimj.201600136
• Fatoretto MB, Lara IAR, Loro AC, Spoto MHF. 2018. Sensory evaluation of dehydrated tomatoes using the proportional odds mixed model. Journal of Food Processing and Preservation 42: e13822. https://doi.org/10.1111/jfpp.13822
  » https://doi.org/10.1111/jfpp.13822
• Ferrari S, Cribari-Neto F. 2004. Beta regression for modelling rates and proportions. Journal of Applied Statistics 31: 799-815. https://doi.org/10.1080/0266476042000214501
  » https://doi.org/10.1080/0266476042000214501
• Figueroa-Zúñiga JI, Arellano-Valle RB, Ferrari SLP. 2013. Mixed beta regression: a Bayesian perspective. Computational Statistics & Data Analysis 61: 137-147. https://doi.org/10.1016/j.csda.2012.12.002
  » https://doi.org/10.1016/j.csda.2012.12.002
• Gadrich T, Pennecchi FR, Kuselman I, Hibbert DB, Semenova AA, Cheow PS. 2022. Ordinal analysis of variation of sensory responses in combination with multinomial ordered logistic regression vs. chemical composition: a case study of the quality of a sausage from different producers. Journal of Food Quality 2022: 4181460. https://doi.org/10.1155/2022/4181460
  » https://doi.org/10.1155/2022/4181460
• Green P, MacLeod CJ. 2016. SIMR: an R package for power analysis of generalized linear mixed models by simulation. Methods in Ecology and Evolution 7: 493-498. https://doi.org/10.1111/2041-210X.12504
  » https://doi.org/10.1111/2041-210X.12504
• Hartig F. 2024. DHARMa: Residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.4.7. Available at: https://CRAN.R-project.org/package=DHARMa [Accessed Jan 10, 2025]
  » https://CRAN.R-project.org/package=DHARMa
• Huang X, Oosterlee CW. 2011. Generalized beta regression models for random loss-given-default. The Journal of Credit Risk 7: 45-70. https://doi.org/10.21314/JCR.2011.150
  » https://doi.org/10.21314/JCR.2011.150
• Junaid A, Khan A, Alrumayh A, Alghamdi FM, Hussam E, Aljohani HM, et al. 2024. Modified two parameter ridge estimator for beta regression model. Journal of Radiation Research and Applied Sciences 17: 100905. https://doi.org/10.1016/j.jrras.2024.100905
  » https://doi.org/10.1016/j.jrras.2024.100905
• Kottaridi K, Milionis A, Demopoulos V, Rigakou A, Nikolaidis V. 2023. A regression analysis method for the prediction of olive oil sensory attributes. Journal of Agriculture and Food Research 12: 100555. https://doi.org/10.1016/j.jafr.2023.100555
  » https://doi.org/10.1016/j.jafr.2023.100555
• Kubinec R. 2023. Ordered beta regression: a parsimonious, well-fitting model for continuous data with lower and upper bounds. Political Analysis 31: 519-536. https://doi.org/10.1017/pan.2022.20
  » https://doi.org/10.1017/pan.2022.20
• Moral RA, Hinde J, Demétrio CG. 2017. Half-normal plots and overdispersed models in R: the hnp package. Journal of Statistical Software 81: 1-23. https://doi.org/10.18637/jss.v081.i10
  » https://doi.org/10.18637/jss.v081.i10
• Morris MD, Moore LM, McKay MD. 2006. Sampling plans based on balanced incomplete block designs for evaluating the importance of computer model inputs. Journal of Statistical Planning and Inference 136: 3203-3220. https://doi.org/10.1016/j.jspi.2005.02.016
  » https://doi.org/10.1016/j.jspi.2005.02.016
• Pereira GHA, Botter DA, Sandoval MC. 2013. A regression model for special proportions. Statistical Modelling 13: 125-151. https://doi.org/10.1177/1471082X13478274
  » https://doi.org/10.1177/1471082X13478274
• Pereira GHA. 2019. On quantile residuals in beta regression. Communications in Statistics: Simulation and Computation 48: 302-316. https://doi.org/10.1080/03610918.2017.1381740
  » https://doi.org/10.1080/03610918.2017.1381740
• Rossiter KL, Young H, Walker SB, Miller M, Dawson DM. 2000. The effects of sugars and acids on consumer acceptability of kiwifruit. Journal of Sensory Studies 15: 241-250. https://doi.org/10.1111/j.1745-459X.2000.tb00269.x
  » https://doi.org/10.1111/j.1745-459X.2000.tb00269.x
• Salinas-Ruíz J, Montesinos López OA, Hernández Ramírez G, Crossa Hiriart J. 2023. Generalized Linear Mixed Models for Repeated Measurements. p. 377-423. In: Cárdenas-Peña D, Franco JD, Escobar LA. eds. Generalized linear mixed models with applications in agriculture and biology. Springer, Cham, Germany. https://doi.org/10.1007/978-3-031-32800-8_9
  » https://doi.org/10.1007/978-3-031-32800-8_9
• Sarzo B, Martínez-Minaya J, Pennino MG, Conesa D, Coll M. 2023. Modelling seabirds biodiversity through Bayesian Spatial Beta regression models: a proxy to inform marine protected areas in the Mediterranean Sea. Marine Environmental Research 185: 105860. https://doi.org/10.1016/j.marenvres.2022.105860
  » https://doi.org/10.1016/j.marenvres.2022.105860
• Silva, RCSN, Minim VPR, Silva AN, Simiqueli AA, Della Lucia SM, Minim LA. 2014. Balanced incomplete block design: an alternative for data collection in the optimized descriptive profile. Food Research International 64: 289-297. https://doi.org/10.1016/j.foodres.2014.06.042
  » https://doi.org/10.1016/j.foodres.2014.06.042
• Sipos L, Nyitrai A, Hitka G, Friedrich LF, Kókai Z. 2021. Sensory panel performance evaluation: comprehensive review of practical approaches. Applied Sciences 11: 11977. https://doi.org/10.3390/app112411977
  » https://doi.org/10.3390/app112411977
• Smithson M, Verkuilen J. 2006. A better lemon squeezer? Maximum-likelihood regression with beta-distributed dependent variables. Psychological Methods 11: 54-71. https://doi.org/10.1037/1082-989x.11.1.54
  » https://doi.org/10.1037/1082-989x.11.1.54
• Stroup WW. 2012. Generalized linear mixed models: modern concepts, methods and applications. CRC Press, Boca Raton, FL, USA. https://doi.org/10.1201/b13151
  » https://doi.org/10.1201/b13151
• Sugumar JK, Guha P. 2022. Comparative study on the hedonic and fuzzy logic-based sensory analysis of formulated soup mix. Future Foods 5: 100115. https://doi.org/10.1016/j.fufo.2022.100115
  » https://doi.org/10.1016/j.fufo.2022.100115
• Tierney L, Kadane JB. 1986. Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association 81: 82-86. https://doi.org/10.1080/01621459.1986.10478240
  » https://doi.org/10.1080/01621459.1986.10478240
• Verkuilen J, Smithson M. 2016. Mixed and mixture regression models for continuous bounded responses using the beta distribution. Journal of Educational and Behavioral Statistics 37: 82-113. https://doi.org/10.3102/1076998610396895
  » https://doi.org/10.3102/1076998610396895
• Yu H, Han B, Rios N, Chen J. 2024. Missing data imputation in balanced construction for incomplete block designs. Mathematics 12: 3419. https://doi.org/10.3390/math12213419
  » https://doi.org/10.3390/math12213419
• Zimprich D. 2010. Modeling change in skewed variables using mixed beta regression models. Research in Human Development 7: 9-26. https://doi.org/10.1080/15427600903578136
  » https://doi.org/10.1080/15427600903578136

Edited by

• Edited by:
  Thomas Kumke https://orcid.org/0000-0001-5492-1612

Publication Dates

• Publication in this collection
  14 Aug 2026
• Date of issue
  2026

History

• Received
  02 Apr 2025
• Accepted
  07 Oct 2025