---
title: 'Everything, altogether, all at once: addressing data challenges when measuring speech intelligibility through entropy scores'
author:
  - name: 
      given: Jose Manuel
      family: Rivera Espejo
    orcid: 0000-0002-3088-2783
    url: https://www.uantwerpen.be/en/staff/jose-manuel-rivera-espejo_23166/
    email: JoseManuel.RiveraEspejo@uantwerpen.be
    corresponding: true
    roles: [conceptualization, methodology, data curation, analysis, writing, software, investigation, validation, visualization]
    affiliation:
      - name: University of Antwerp
        department: Department of Training and Education Sciences
        group: Edubron
  - name: 
      given: Sven
      family: De Maeyer
      #non-dropping-particle: De
    orcid: 0000-0003-2888-1631
    url: https://www.uantwerpen.be/en/staff/sven-demaeyer/
    email: sven.demaeyer@uantwerpen.be
    corresponding: false
    roles: [conceptualization, methodology, editing, supervision, project administration, funding, investigation, resources]
    affiliation:
      - name: University of Antwerp
        department: Department of Training and Education Sciences
        group: Edubron
  - name: 
      given: Steven
      family: Gillis
    orcid: 
    url: https://www.uantwerpen.be/nl/personeel/steven-gillis/
    email: steven.gillis@uantwerpen.be
    corresponding: false
    roles: [conceptualization, methodology, editing, supervision, project administration, funding, investigation, resources]
    affiliation:
      - name: University of Antwerp
        department: Department of Linguistics
        group: Centre for computational linguistics, psycholinguistics, and sociolinguistics (CLiPS)
funding: 
  statement: "The project was founded through the Research Fund of the University of Antwerp (BOF)."
keywords:
  - Bayesian analysis 
  - speech intelligibility 
  - bounded outcomes 
  - clustering 
  - measurement error
  - outliers
  - heteroscedasticity
  - generalized linear latent and mixed models
  - robust regression models.
abstract: |
  When investigating unobservable, complex traits, data collection and aggregation processes can introduce distinctive features to the data such as boundedness, measurement error, clustering, outliers and heteroscedasticity. Failure to collectively address these features can result in statistical challenges that prevent the investigation of hypotheses regarding these traits. This study aimed to demonstrate the efficacy of the Bayesian Beta-proportion Generalized Linear Latent and Mixed Model (Beta-proportion GLLAMM) [@Rabe_et_al_2004a; @Rabe_et_al_2004b; @Rabe_et_al_2004c; @Skrondal_et_al_2004a] in handling data features when exploring research hypotheses concerning speech intelligibility. To achieve this objective, the study reexamined data from transcriptions of spontaneous speech samples initially collected by @Boonen_et_al_2023. The data were aggregated into entropy scores. The research compared the prediction accuracy of the Beta-proportion GLLAMM with the Normal Linear Mixed Model (LMM) [@Holmes_et_al_2019] and investigated its capacity to estimate a latent intelligibility from entropy scores. The study also illustrated how hypotheses concerning the impact of speaker-related factors on intelligibility can be explored with the proposed model. The Beta-proportion GLLAMM was not free of challenges; its implementation required formulating assumptions about the data-generating process and knowledge of probabilistic programming languages, both central to Bayesian methods. Nevertheless, results indicated the superiority of the model in predicting empirical phenomena over the Normal LMM, and its ability to quantify a latent potential intelligibility. Additionally, the proposed model facilitated the exploration of hypotheses concerning speaker-related factors and intelligibility. Ultimately, this research has implications for researchers and data analysts interested in quantitatively measuring intricate, unobservable constructs while accurately predicting the empirical phenomena.
date: last-modified
bibliography: references.bib
notebook-links: global
citation: 
  type: article-journal
  container-title: "Behavior Research Methods"
  doi: "10.3758/s13428-024-02457-6"
  # url: https://example.com/summarizing-output
execute: 
  cache: true
  # eval: true
  echo: false
  # output: true
  # include: true
  warning: false
  error: false
  message: false
---

```{r}
#| results: hide
#| label: code-package

# load packages
libraries = c('RColorBrewer','stringr','dplyr','tidyverse',
              'rechape2','knitr','kableExtra',
              'rstan','StanHeaders','runjags','rethinking')
sapply(libraries, require, character.only=T)

```


```{r}
#| label: code-usd

# load functions
main_dir = '/home/josema/Desktop/1. Work/1 research/PhD Antwerp/#thesis/paper1'
source( file.path( main_dir, 'paper1_manuscript/code', 'user-defined-functions.R') )

```


# Introduction {#sec-introduction}

Intelligibility is at the core of successful, felicitous communication. Thus, being able to speak intelligibly is a major achievement in language acquisition and development. Moreover, intelligibility is considered to be the most practical index to assess competence in oral communication [@Kent_et_al_1994]. Consequently, it serves as a key indicator for evaluating the effectiveness of various interventions like speech therapy or cochlear implantation [@Chin_et_al_2012].

The notion of speech intelligibility may appear deceptively simple, yet it is an intricate concept filled with inherent challenges in its assessment. Intelligibility refers to the extent to which a listener can accurately recover the words in a speaker’s acoustic signal [@Freeman_et_al_2017; @vanHeuven_2008; @Whitehill_et_al_2004]. Furthermore, achieving intelligible spoken language requires to master all core components of speech perception, cognitive processing, linguistic knowledge, and articulation [@Freeman_et_al_2017]. Hence, it is unsurprising that its accurate measurement faces challenges [@Kent_et_al_1989]. These challenges arise from the interplay of the attributes of the communicative environment such as background noise [@Munro_1998], with features of the speaker like speaking rate [@Munro_et_al_1998] or accent [@Jenkins_2000; @Ockey_et_al_2016], and characteristics of the listener like vocabulary proficiency or hearing ability [@Varonis_et_al_1985].

While several approaches have been proposed to assess intelligibility, they commonly rely on two types of speech samples: read-aloud or imitated, and spontaneous speech samples. Most studies favor read-aloud or imitated speech samples due to the substantial control they offer in selecting stimuli for intelligibility assessment. Additionally, these types of speech facilitate a direct and unambiguous comparison between a defined word target, produced by a speaker, and the listener’s identification of it, as exemplified by multiple studies such as @Castellanos_et_al_2014, @Chin_et_al_2012, @Chin_et_al_2014, @Freeman_et_al_2017, @Khwaileh_et_al_2010, and @Montag_et_al_2014. However, it has been demonstrated that these controlled speech samples exhibit limited efficacy in predicting intelligibility among hearing-impaired individuals [@Cox_et_al_1989; @Ertmer_2011]. In contrast, spontaneous speech samples offer a more ecologically valid approach to assess intelligibility, resembling everyday informal speech more than read-aloud or imitated speech samples [@Boonen_et_al_2023]. However, due to the uncertainty surrounding the speaker's intended word production, it is unfeasible to establish a word target for these samples [@Flipsen_2006; @Lagerberg_et_al_2014]. This renders conventional accuracy metrics from imitated speech, such as the percentage of read or imitated words, impractical [@Boonen_et_al_2023].

Yet, various metrics of intelligibility can still be derived from transcriptions of spontaneous speech samples, including the percentage of (un)intelligible words or syllables [@Flipsen_2006; @Lagerberg_et_al_2014], as well as entropy scores [@Boonen_et_al_2023]. In the latter approach, listeners transcribe orthographically spontaneous speech samples produced by various speakers. These transcriptions are then aggregated into entropy scores, where lower scores indicate a higher degree of agreement among the listeners transcriptions and, consequently, higher intelligibility, while higher scores suggest lower intelligibility due to a lower degree of agreement in the transcriptions [@Boonen_et_al_2023; @Faes_et_al_2022]. Notably, the aggregation procedure assumes that speech samples are considered intelligible if all listeners decode them in the same manner. These scores have been instrumental in examining differences in speakers’ speech intelligibility, particularly between children with normal hearing and those with cochlear implants [@Boonen_et_al_2023]. 

However, despite the entropy scores' potential as a fine-grained metric of intelligibility, as proposed by @Boonen_et_al_2023, they exhibit a statistical complexity that cautions researchers against treating them as straightforward indices of intelligibility. This complexity emerges from the processes of data collection and transcription aggregation, endowing the scores with four distinctive features: boundedness, measurement error, clustering, and the possible presence of outliers and heteroscedasticity. Firstly, entropy scores are confined to the interval between zero and one, a phenomenon known as boundedness. Boundedness refers to the restriction of data values within specific bounds or intervals, beyond which they cannot occur [@Lebl_2022]. Secondly, entropy scores are assumed to be a manifestation of a speaker’s intelligibility, with this intelligibility being the primary factor influencing the observed scores. This issue is commonly referred to as measurement error, signifying the disparity between the observed values of a variable, recorded under similar conditions, and some fixed *true value* which is not directly observable [@Everitt_et_al_2010]. Thirdly, due to the repeated assessment of speakers through multiple speech samples, the scores exhibit clustering. Clustering occurs when outcomes stem from repeated measurements of the same individual, location, or time [@McElreath_2020]. Lastly, driven by speech samples with entropy scores located at the extreme of the bounds, and the presence of more than one population in the data (i.e., normal hearing versus hearing-impaired speakers), the scores may exhibit a potential for outliers and heteroscedasticity. Outliers are observations that markedly deviate from other sample data points where they occur [@Grubbs_1969], while heteroscedasticity occurs when the outcome’s variance depends on the values of another variable [@Everitt_et_al_2010].

Failure to collectively address these data features can result in numerous statistical challenges that might hamper the researcher’s ability to investigate intelligibility. Notably, neglecting boundedness can, at best, lead to underfitting and, at worst, to misspecification. Underfitting occurs when statistical models fail to capture the underlying data patterns, potentially generating predictions outside the data range, thus hindering the model’s ability to generalize when confronted with new data. Conversely, misspecification, which is marked by a poor representation of relevant aspects of the true data in the model's functional form, can lead to inconsistent and less precise parameter estimates [@Everitt_et_al_2010]. Additionally, overlooking issues such as measurement error, clustering, outliers, or heteroscedasticity can lead to biased and less precise parameter estimates [@McElreath_2020], ultimately diminishing the statistical power of models and increasing the likelihood of committing type I or type II errors when addressing research inquiries. Type I error results when the null hypothesis is falsely rejected, while Type II error that results when the null hypothesis is falsely accepted [@Everitt_et_al_2010].

In computational statistics and data analysis, several models have been developed to address some of these data features individually and, at times, collectively. For instance, @Ferrari_et_al_2004 and @Simas_et_al_2010 initially introduced and expanded beta regression models to handle outcomes constrained within the unit interval. Subsequently, @Figueroa-Zuniga_et_al_2013 extended these models to address data clustering. Over time, beta regression models have evolved to accommodate clustering and measurement errors in covariates, as demonstrated by @Carrasco_et_al_2012 and @Figueroa-Zuniga_et_al_2018. Furthermore, robust versions of these models have been proposed to account for other statistical data issues, such as outliers and heteroscedasticity, as seen in @Bayes_et_al_2012 and @Figueroa-Zuniga_et_al_2021. Robust models are a general class of statistical procedures designed to reduce the sensitivity of the parameter estimates to mild or moderate departures of the data from the model's assumptions [@Everitt_et_al_2010]. Ultimately, the work of Rabe-Hesketh and colleagues introduced the Generalized Linear Latent and Mixed Model (GLLAMM) [@Rabe_et_al_2004a; @Rabe_et_al_2004b; @Rabe_et_al_2004c; @Skrondal_et_al_2004a], a unified framework that can simultaneously tackle all of the aforementioned data features.

All of these models have found moderate adoption in various fields, including speech communication [@Boonen_et_al_2023], psychology [@Unlu_et_al_2017], cognition [@Verkuilen_et_al_2012; @Lopes_et_al_2023], education [@Pereira_et_al_2020], health care [@Ghosh_2019; @Kangmennaang_et_al_2023], chemistry [@de_Brito_et_al_2021], and policy analysis [@Dieteren_et_al_2023; @Choi_2023; @Zhang_et_al_2023]. Specifically, in the domain of speech communication, @Boonen_et_al_2023 addressed data clustering within the context of intelligibility research. Conversely, @de_Brito_et_al_2021 and @Kangmennaang_et_al_2023 concentrated on tackling non-normal bounded data with measurement error in covariates, within the context of chemical reactions and health care access, respectively. Remarkably, despite these individual efforts, there is, to the best of the authors' knowledge, no study comprehensively addressing all of these data features in a principled way, while also transparently and systematically documenting the Bayesian estimation of the resulting statistical models. 

This study employed Bayesian procedures for three main reasons. Firstly, prior research have consistently demonstrated the superiority of Bayesian methods over frequentist methods, especially with complex and overparameterized models [@Baker_1998; @Kim_et_al_1999], such as the GLLAMM used in this study. Overparameterized models are those with more parameters than observations for estimation [@Everitt_et_al_2010]. Secondly, the Bayesian approach enabled the incorporation of prior information, thereby constraining certain parameters within specified bounds. This feature addressed issues such as non-convergence or improper parameter estimation common in complex models under frequentist methods [@Martin_et_al_1975; @Seaman_et_al_2011]. An example is the estimation of negative variances for random effects in hierarchical models [@Holmes_et_al_2019], a problem resolved in this study through the utilization of prior distributions. Lastly, Bayesian methods have exhibited proficiency in drawing inferences from small sample sizes [@Baldwin_et_al_2013; @Lambert_et_al_2005; @Depaoli_2014]. This feature of the Bayesian methods holds relevance for this study, as it also grapples with a small sample size, where reliance on the asymptotic properties of frequentist methods may not be justified.


## Research questions {#sec-I-RQ}

Considering the imperative need to comprehensively address all features of the data when investigating unobservable and complex traits, this investigation aimed to demonstrate the efficacy of the Generalized Linear Latent and Mixed Model (GLLAMM) in handling entropy score features when exploring research hypotheses concerning speech intelligibility. To achieve this objective, the study reexamined data originating from transcriptions of spontaneous speech samples, initially collected by @Boonen_et_al_2023. The data was aggregated into entropy scores and subjected to modelling through the Bayesian Beta-proportion GLLAMM.

To address the primary objective, the study posed three key research questions. First, given the importance of accurate predictions in developing useful practical models and testing research hypotheses [@Shmueli_et_al_2012], *Research Question 1 (RQ1)* evaluated whether the Beta-proportion GLLAMM yielded more accurate predictions than the widely used Normal Linear Mixed Model (LMM) [@Holmes_et_al_2019]. Second, acknowledging that intelligibility is an unobservable, intricate concept and a key indicator of oral communication competence [@Kent_et_al_1994], *Research Question 2 (RQ2)* investigated how the proposed model can estimate speakers’ latent intelligibility from manifest entropy scores. Thirdly, recognizing that research involves developing and comparing hypotheses, *Research Question 3 (RQ3)* illustrated how these research hypotheses can be examined within the model’s framework. Specifically, RQ3 assessed the influence of speaker-related factors on the newly estimated latent intelligibility.

Ultimately, this study offers researchers studying speech intelligibility through entropy scores and those in similar or different fields facing analogous data challenges with a statistical tool that improves upon current research models. This tool assess the predictability of empirical phenomena and develops a quantitative measure for the latent variable of interest. This quantitative measure, in turn, facilitates the appropriate comparison of existing hypotheses related to the latent variable, and even encourages the formulation of new ones.


# Methods {#sec-methods}

## Data {#sec-M-D}

The data comprised the transcriptions of spontaneous speech samples originally collected by @Boonen_et_al_2023. The data is not publicly available due to privacy restrictions. Nonetheless, the data can be provided by the corresponding author upon reasonable request.

```{r}
#| label: code-data

data_dir = '/home/josema/Desktop/1. Work/1 research/PhD Antwerp/#thesis/#data/final'
data_nam = "data_H_list.RData"
model_data = file.path(data_dir, data_nam )
load( model_data )
# str(dlist)

var_int = c('bid','cid','uid','wid','HS','A','Am','Hwsib')
data_H = data.frame(dlist[var_int])
# str(data_H)

```


### Speakers {#sec-M-S}

@Boonen_et_al_2023 selected $32$ speakers, comprising $16$ normal hearing children (NH) and $16$ hearing-impaired children with cochlear implants (HI/CI). At the time of the collection of the speech samples, the NH group were between $68$ and $104$ months old ($M=86.3$, $SD=9.0$), while HI/CI group were between $78$ and $98$ months old ($M=86.3$, $SD=6.7$). All children were native speakers of Belgian Dutch.


### Speech samples {#sec-M-SS}

Boonen and colleagues selected speech samples from a large corpus of children’s spontaneously spoken speech recordings. These recordings were made in Belgian Dutch and obtained while the children narrated a story prompted by the picture book "Frog, Where Are You?" [@Mayer_1969] to a caregiver ‘unfamiliar with the story’. Before the actual recording, the children were allowed to skim over the booklet and examine the pictures. Prior to the selection of the samples, the recordings were orthographically transcribed using the CHAT format in the CLAN editor [@MacWhinney_2020]. These transcriptions were exclusively used in the selection of appropriate speech samples. To ensure the quality of the selection, Boonen and colleagues excluded sentences containing syntactically ill-formed or incomplete statements, with background noise, crosstalk, long hesitations, revisions, or non-words. Finally, ten speech samples were randomly chosen for each of the $32$ selected speakers. Each of these samples comprised a single sentence with a length of three to eleven words ($M=7.1$, $SD=1.1$). The process resulted in a total of $320$ selected sentences collectively comprising $2,263$ words.


### Listeners {#sec-M-L}

Boonen and colleagues recruited $105$ students from the University of Antwerp. All participants were native speakers of Belgian Dutch and reported no history of hearing difficulties or prior exposure to the speech of hearing-impaired speakers.


### Transcription task and entropy scores {#sec-M-TS}

@Boonen_et_al_2023 distributed the $320$ speech samples and $105$ listeners into five blocks through random allocation. Each block comprised $21$ listeners and $64$ sentences with no overlap between the blocks. The listeners were tasked with transcribing each sentence, which were presented to them in a random order. This resulted in a total of $47,514$ transcribed words from the original $2,263$ words available in the speech samples. These orthographic transcriptions were automatically aligned with a python script [@Boonen_et_al_2023], at the sentence level in a column-like grid structure like the one presented in @tbl-alignment. This alignment process was repeated for each sentence from every speaker, and the output was manually checked and adjusted (if needed) in order to appropriately align the words. For more details on the random assignment and alignment procedures refer to the original authors.

Next, the aligned transcriptions were aggregated by listener, yielding $2,263$ entropy scores, one score per word for every sentence. The entropy scores were calculated following Shannon’s formula [-@Shannon_1948]:

$$
H_{wsib} = \frac{ -\left[ \sum_{k=1}^{K}  p_{k} \cdot log_{2}(p_{k}) \right] }{ log_{2}(J)}
$$ {#eq-entropy}

where $H_{wsib}$ denotes the entropy scores confined to an interval between zero and one, with $w$ defining the word index, $s$ the sentence index, $i$ the speaker index, and $b$ the block index. In addition, $K$ describes the number of different word types within transcriptions, and $J$ defines the total number of word transcriptions. Notice that by design, the total number of word transcriptions $J$ corresponds with the number of listeners per block, i.e., $21$ listeners. Lastly, $p_{k} = \sum_{j=1}^{J} 1(T_{jk}) / J$ denotes the proportion of word types within transcriptions, with $1(T_{jk})$ describing an indicator function that takes the value of one when the word type $k$ is present in the transcription $j$. See @sec-appA for an example of how entropy scores are computed. 

These entropy scores served as the outcome variable, capturing agreement or disagreement among listeners’ word transcriptions. Lower scores indicated a higher degree of agreement between transcriptions and therefore higher intelligibility, while higher scores indicated lower intelligibility, due to a lower degree of agreement in the transcriptions [@Boonen_et_al_2023; @Faes_et_al_2022]. Furthermore, no score was excluded from the modelling process using univariate procedures, rather, the identification of highly influential observations was performed within the context of the proposed models, as recommended by @McElreath_2020.

+---------------+---------+---------+---------+---------+---------+
| Transcription | Words   |         |         |         |         |
+---------------+---------+---------+---------+---------+---------+
| Number        | 1       | 2       | 3       | 4       | 5       |
+:=============:+:=======:+:=======:+:=======:+:=======:+:=======:+
| 1             | de      | jongen  | ziet    | een     | kikker  |
+---------------+---------+---------+---------+---------+---------+
|               | the     | boy     | sees    | a       | frog    |
+---------------+---------+---------+---------+---------+---------+
| 2             | de      | jongen  | ziet    | de      | [X]     |
+---------------+---------+---------+---------+---------+---------+
|               | the     | boy     | sees    | the     | [X]     |
+---------------+---------+---------+---------+---------+---------+
| 3             | de      | jongen  | zag     | [B]     | kokkin  |
+---------------+---------+---------+---------+---------+---------+
|               | the     | boy     | saw     | [B]     | cook    |
+---------------+---------+---------+---------+---------+---------+
| 4             | de      | jongen  | zag     | geen    | kikkers |
+---------------+---------+---------+---------+---------+---------+
|               | the     | boy     | saw     | no      | frogs   |
+---------------+---------+---------+---------+---------+---------+
| 5             | de      | hond    | zoekt   | een     | [X]     |
+---------------+---------+---------+---------+---------+---------+
|               | the     | dog     | searches| a       | [X]     |
+---------------+---------+---------+---------+---------+---------+
| **Entropy**   | $0$     |$0.3109$ |$0.6555$ |$0.8277$ |$1$      |
+---------------+---------+---------+---------+---------+---------+

: Hypothetical alignment of word transcriptions and entropy scores. **_Note:_** Extracted from Boonen et al. [-@Boonen_et_al_2023], and slightly modified for illustrative purposes. Entropy scores were calculated from words of the first sentence, produced by the first speaker assigned to the first block, and transcribed by five listeners $\left( s=1, i=1, b=1, J=5 \right)$. Transcriptions are in Belgian Dutch followed by their English translation. *[B]* represent a blank space, and *[X]* an unidentifiable speech. {#tbl-alignment .striped .hover}


## Statistical models {#sec-M-SM}

This section articulates the probabilistic formalism of both the Normal LMM and the proposed Beta-proportion GLLAMM. Subsequently, it details the set of fitted models and the estimation procedure, along with the criteria employed to assess the quality of the Bayesian inference results. Lastly, the section outlines the methodology employed for model comparison.


### Normal LMM {#sec-M-SM-NLMM}

The general mathematical formalism of the Normal LMM posits that the likelihood of the (manifest) entropy scores  follow a normal distribution, i.e.

$$
H_{wsib} \sim \text{Normal} \left( \mu_{sib}, \sigma_{i} \right)
$$ {#eq-normal-LMM-likelihood}

where $\mu_{sib}$ represents the average entropy at the word-level and $\sigma_{i}$ denotes the standard deviation of the average entropy at the word-level, varying for each speaker. Given the clustered nature of the data, $\mu_{sib}$ is defined by the linear combination of individual characteristics and several random effects:

$$
\mu_{sib} = \alpha + \alpha_{HS[i]} + \beta_{A, HS[i]} (A_{i} - \bar{A}) + u_{si} + e_{i} + a_{b}
$$ {#eq-normal-LMM-linearpred}

where $HS_{i}$ and $A_{i}$ denote the hearing status and chronological age of speaker $i$, respectively. Additionally, $\alpha$ denotes the general intercept, $\alpha_{HS[i]}$ represents the average entropy for each hearing status group, and $\beta_{A,HS[i]}$ denotes the evolution of the average entropy per unit of chronological age $A_{i}$ for each hearing status group. Furthermore, $u_{si}$ denotes the sentence-speaker random effects measuring the unexplained entropy variability within sentences for each speaker, $e_{i}$ denotes the speaker random effects describing the unexplained entropy variability between speakers, and $a_{b}$ denotes the block random effects assessing the unexplained variability between experimental blocks. 

Several notable features of the Normal LLM can be discerned from the equations. Firstly, @eq-normal-LMM-likelihood indicates that the variability of the average entropy at the word level can differ for each speaker, enhancing the model *robustness* to mild or moderate data departures from the normal distribution assumption, such as in the presence of heteroscedasticity or outliers. Secondly, @eq-normal-LMM-linearpred reveals that the model assumes that no transformation is applied to the relationship between the average entropy and the linear combination of speakers' characteristics. This is commonly known as a direct link function. In addition, the equation indicates that chronological age is *centered* around the minimum chronological age in the sample $\bar{A}$. The *centering* procedure prevents the interpretation of parameters outside the range of chronological ages available in the data [@Everitt_et_al_2010]. Also, the equation implies the model considers separate intercept and separate age slopes for each hearing status group, i.e., $\alpha_{HS[i]}$ and $\beta_{A, HS[i]}$ for NH and HI/CI speakers, respectively. Lastly, the presence of a general intercept $\alpha$ in the equation reveals that the model is overparameterized. Although the estimation of overparameterized models is only possible under Bayesian methods, their estimation does not violate any statistical principle [@McElreath_2020, pp. 345]. In contrast, in this study, the overparameterized model facilitates: (1) the comparison between the specific parameter interpretations of the Normal LMM and the Beta-proportion GLLAMM, with  $\alpha$ serving no particular purpose in the former case, and (2) the assignment of prior distributions.


### Beta-proportion GLLAMM {#sec-M-SM-BGLLAMM}

The general mathematical formalism of the proposed Beta-proportion GLLAMM comprises four components: a response model likelihood, a linear predictor, a link function, and a structural model. The likelihood of the response model posits that entropy scores follow a Beta-proportion distribution,

$$
H_{wsib} \sim \text{BetaProp} \left( \mu_{ib}, M_{i} \right)
$$ {#eq-beta-GLLAMM-likelihood}

where $\mu_{ib}$ denotes the average entropy at the word-level and $M_{i}$ signifies the *dispersion* of the average entropy at the word-level, varying for each speaker. Additionally, $\mu_{ib}$ is defined as, 

$$
\mu_{ib} = \text{logit}^{-1}[ a_{b} - SI_{i} ]
$$ {#eq-beta-GLLAMM-linpred}

where $\text{logit}^{-1}(x) = exp(x) / (1+exp(x))$ is the inverse-logit link function, $a_{b}$ denotes the block random effects, and $SI_{i}$ describes the speaker's latent *potential intelligibility*. Conversely, the structural equation model relates the speakers' latent potential intelligibility to the individual characteristics:

$$
SI_{i} = \alpha + \alpha_{HS[i]} + \beta_{A, HS[i]} (A_{i} - \bar{A}) + e_{i} + u_{i}
$$ {#eq-beta-GLLAMM-structural}

where $\alpha$ defines the general intercept, $\alpha_{HS[i]}$ denotes the potential intelligibility for different hearing status groups, and $\beta_{A,HS[i]}$ indicates the evolution of potential intelligibility per unit of chronological age for each hearing status group. Furthermore, $e_{i}$ represents speakers block effects, describing unexplained potential intelligibility variability between speakers, and $u_{i} = \sum_{s=1}^{S} u_{si}/S$ denotes sentence random effects, assessing the average unexplained potential intelligibility variability within sentences for each speaker, with $S$ denoting the total number of sentences per speaker.

Several features are evident in the probabilistic representation of the model. Firstly, akin to the Normal LMM, @eq-beta-GLLAMM-likelihood reveals that the *dispersion* of average entropy at the word level can differ for each speaker. This enhances the model's robustness to mild or moderate data departures from the beta-proportion distribution assumption. Secondly, in contrast with the Normal LMM, @eq-beta-GLLAMM-linpred shows the potential intelligibility of a speaker has a negative non-linear relationship with the entropy scores. The negative relationship explicitly highlights the inverse relationship between intelligibility and entropy, while the non-linear relationship maps the unbounded linear predictor to the bounded limits of the entropy scores. Thirdly, in contrast with the Normal LMM, @eq-beta-GLLAMM-structural demonstrates that the structural parameters are interpretable in terms of the latent potential intelligibility scores, where the scale of the latent trait is set by the general intercept $\alpha$, as it is required in latent variable models [@Depaoli_2021]. Furthermore, the equation implies the model also considers separate intercept and separate age slopes for each hearing status group, i.e., $\alpha_{HS[i]}$ and $\beta_{A, HS[i]}$ for NH and HI/CI speakers, respectively. In addition, it indicates that chronological age is *centered* around the minimum chronological age in the sample $\bar{A}$. Lastly, the equation also reveals that the intelligibility scores have two sources of unexplained variability. The term $e_{i}$ represents inherent differences in potential intelligibility among different speakers. The term $u_{i}$ assumes that different sentences measure potential intelligibility differently due to variations in word difficulties and their interplay within the sentence.


### Prior distributions {#sec-M-SM-P}

Bayesian procedures require the incorporation of priors. Priors are probability distributions summarizing the information about known or assumed parameters prior to observing any empirical data [@Everitt_et_al_2010]. Upon observing empirical data, these priors undergo updating to posterior distributions following Bayes' rule [@Jeffreys_1998]. In cases requiring greater modelling flexibility, a more refined representation of the parameters’ priors can be defined in terms of hyperparameters and hyperpriors. *Hyperparameters* refer to parameters indexing a family of possible prior distributions for the original parameter, while *hyperpriors* are prior distributions for such hyperparameters [@Everitt_et_al_2010].

This study established priors and hyperpriors for the parameters of both the Normal LMM and the Beta-proportion GLLAMM using prior predictive simulations. This procedure entails the semi-independent simulation of parameters, which are subsequently transformed into simulated data values according to the models' specifications. The procedure aims to establish meaningful priors and comprehend their implications within the context of the model before incorporating any information derived from empirical data [@McElreath_2020]. For reader inspection, the prior predictive simulations are provided in the accompanying digital walk-through document (see @sec-M-SM-OS Open Science Statement).


#### Normal LMM {#sec-M-SM-P-NLMM}

For the parameters of the Normal LMM, non-informative priors and hyperpriors were established to align with analogous model assumptions in frequentist methods. A *non-informative* prior reflects the distributional commitment of a parameter to a wide range of values within a specific parameter space [@Everitt_et_al_2010]. The specified priors were as follows:

$$
\begin{aligned}
r_{S} &\sim \text{Exponential}\left( 2 \right) \\ 
\sigma_{i} &\sim \text{Exponential}\left( r_{S} \right) \\
m_{i} &\sim \text{Normal} \left( 0, 0.05 \right) \\
s_{i} &\sim \text{Exponential} \left( 2 \right) \\
e_{i} &\sim \text{Normal} \left( m_{i}, s_{i} \right) \\
m_{b} &\sim \text{Normal} \left( 0, 0.05 \right) \\
s_{b} &\sim \text{Exponential} \left( 2 \right) \\
a_{b} &\sim \text{Normal} \left( m_{b}, s_{b} \right) \\
\alpha &\sim \text{Normal} \left( 0, 0.05 \right) \\
\alpha_{HS[i]} &\sim \text{Normal} \left( 0, 0.2 \right) \\
\beta_{A,HS[i]} &\sim \text{Normal} \left( 0, 0.1 \right)
\end{aligned} 
$$ {#eq-normal-LMM-priors}

#### Beta-proportion GLLAMM {#sec-M-SM-P-BGLLAMM}

For the parameters of the Beta-proportion GLLAMM, weakly informative priors and hyperpriors were established. *Weakly informative priors* reflect the distributional commitment of a parameter to a weakly constraint range of values within a realistic parameter space [@McElreath_2020]. The specified priors were as follows:

$$
\begin{aligned}
r_{M} &\sim \text{Exponential}\left( 2 \right) \\ 
M_{i} &\sim \text{Exponential}\left( r_{M} \right) \\
m_{i} &\sim \text{Normal} \left( 0, 0.05 \right) \\
s_{i} &\sim \text{Exponential} \left( 2 \right) \\
e_{i} &\sim \text{Normal} \left( m_{i}, s_{i} \right) \\
m_{b} &\sim \text{Normal} \left( 0, 0.05 \right) \\
s_{b} &\sim \text{Exponential} \left( 2 \right) \\
a_{b} &\sim \text{Normal} \left( m_{b}, s_{b} \right) \\
\alpha &\sim \text{Normal} \left( 0, 0.05 \right) \\
\alpha_{HS[i]} &\sim \text{Normal} \left( 0, 0.3 \right) \\
\beta_{A,HS[i]} &\sim \text{Normal} \left( 0, 0.1 \right)
\end{aligned} 
$$ {#eq-beta-GLLAMM-priors}


+-------+---------+--------------+---------+----------------+-------------+-------------------+
|       | Model   | Entropy      | Robust  | Fixed effects  |             |                   |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| Model | type    | distribution | feature | $\beta_{HS[i]}$| $\beta_{A}$ | $\beta_{A,HS[i]}$ |
+:=====:+:=======:+:============:+:=======:+:==============:+:===========:+:=================:+
| 1     | LMM     | Normal       | No      | No             | No          | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 2     | LMM     | Normal       | No      | Yes            | Yes         | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 3     | LMM     | Normal       | No      | Yes            | No          | Yes               |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 4     | LMM     | Normal       | Yes     | No             | No          | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 5     | LMM     | Normal       | Yes     | Yes            | Yes         | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 6     | LMM     | Normal       | Yes     | Yes            | No          | Yes               |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 7     | GLLAMM  | BetaProp     | No      | No             | No          | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 8     | GLLAMM  | BetaProp     | No      | Yes            | Yes         | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 9     | GLLAMM  | BetaProp     | No      | Yes            | No          | Yes               |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 10    | GLLAMM  | BetaProp     | Yes     | No             | No          | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 11    | GLLAMM  | BetaProp     | Yes     | Yes            | Yes         | No                |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
| 12    | GLLAMM  | BetaProp     | Yes     | Yes            | No          | Yes               |
+-------+---------+--------------+---------+----------------+-------------+-------------------+
: Fitted models. **_Note:_** *Yes* indicates the feature or parameter is included in the model. {#tbl-fitted .striped .hover}


### Fitted models {#sec-M-SM-FM}

This study evaluated the comparative predictive capabilities of both the Normal LMM and the Beta-proportion GLLAMM (RQ1) while simultaneously examined various formulations regarding how speaker-related factors influence intelligibility (RQ3). In this context, the predictive capabilities of the models were intricately connected to these formulations. As a result, the study required fitting $12$ different models, each representing a specific manner to investigate one or both research questions. The models comprised six versions of both the Normal LMM and the Beta-proportion GLLAMM. The differences among the models hinged on (1) whether they addressed data clustering in conjunction with measurement error, denoted as the model type, (2) the assumed distribution for the entropy scores, which aimed to handle boundedness, (3) whether the model incorporated a robust feature to address mild or moderate departures of the data from distributional assumptions, and (4) the inclusion or exclusion of speaker-related factors in the models. A detailed overview of the fitted models is available in @tbl-fitted. 


### Estimation and chain quality {#sec-M-SM-CQ}

The models were estimated using `R` version 4.2.2 [@R_2015] and `Stan` version 2.26.1 [@Stan_2020]. Four Markov chains were implemented for each parameter, each with distinct starting values. Each chain underwent $4,000$ iterations, where the first $2,000$ serving as a warm-up phase and the remaining $2,000$ were considered samples from the posterior distribution. Verification of stationarity, convergence, and mixing for the parameter chains involved graphical analysis and diagnostic statistics. Graphical analysis utilized trace, trace-rank, and autocorrelation plots (ACF). Diagnostic statistics included the *potential scale reduction factor statistics* $\widehat{\text{R}}$ with a cut-off value of $1.05$ [@Vehtari_et_al_2021]. Furthermore, to confirm whether the parameters posterior distributions were generated with a sufficient number of uncorrelated sampling points, each posterior distribution density plot was inspected along with their effective sample size statistics $n_{\text{eff}}$ [@Gelman_et_al_2014]. 

In general, both graphical analysis and diagnostic statistics indicated that all chains exhibited low to moderate autocorrelation, explored the parameter space in a seemingly random manner, and converged to a constant mean and variance in their post-warm-up phase. Moreover, the density plots and  statistics collectively confirmed that all posterior distributions were unimodal distributions with values centered around a mean, generated with a satisfactory number of uncorrelated sampling points, making substantive sense compared to the models' prior beliefs. The trace, trace-rank, ACF, and distribution density plots, along with $\widehat{\text{R}}$ and $n_{\text{eff}}$ statistics, are provided in the accompanying digital walk-through document for reader inspection (see @sec-M-SM-OS Open Science Statement).


### Model comparison {#sec-M-SM-MC}

This study compared the fitted models using three criteria: the deviance information criterion (`DIC`) introduced by @Spiegelhalter_et_al_2002, the widely applicable information criterion (`WAIC`) proposed by @Watanabe_2013, and the Pareto Smoothing Importance Sampling criterion (`PSIS`) developed by @Vehtari_et_al_2017. These criteria score models in terms of deviations from *perfect* predictive accuracy, with smaller values indicating less deviation [@McElreath_2020]. Deviations from *perfect* predictive accuracy serve as the closest estimate for the Kullback-Leibler divergence [@Kullback_et_al_1951], which measures the degree to which a probabilistic model accurately represents the *true* distribution of the data. Specifically, `DIC` measures in-sample deviations, while `WAIC` and `PSIS` offer an approximate measure of out-of-sample deviations.

`WAIC` and `PSIS` are regarded as full Bayesian criteria because they encompass all the information contained in the parameter’s posterior distribution, effectively integrating and reporting the inherent uncertainty in predictive accuracy estimates. In addition to predictive accuracy, `PSIS` offers an extra benefit by identifying highly influential data points. To achieve this, the criterion employs a built-in warning system that flags observations that make out-of-sample predictions unreliable. The rationale is that observations that are relatively unlikely, according to the model, exert more influence and render predictions less reliable compared to those that are relatively expected [@McElreath_2020].

However, since researchers are mostly interested in comparing candidate models, it is the distance between the models that is useful, rather than the absolute value of the criteria [see @McElreath_2020, pp. 209, 223-224]. Therefore, this study utilized the differences in `WAIC` and `PSIS` (`dWAIC` and `dPSIS`, respectively) to evaluate how distinct our probabilistic models are from each other, and which one is closer to the *true* distribution of the data. Additionally, while `DIC`, `WAIC` and `PSIS` provide *approximately correct* estimates for the expected accuracy, the criteria are also subject to uncertainty due to the specific sample over which they are computed [see @McElreath_2020, pp. 223]. Thus, this uncertainty should also be taken into account for the criteria and their comparisons. Consequently, this study also presented the associated uncertainty for both criteria calculated as `WAIC` $\pm 1 \cdot$`SE`, `PSIS` $\pm 1 \cdot$`SE`, `dWAIC` $\pm 1 \cdot$`dSE` and `dPSIS` $\pm 1 \cdot$`dSE`. Lastly, this research also reported the models' complexity penalization, as well as their associated weight of evidence. The complexity penalization values `pWAIC` and `pPSIS` are roughly associated with the models' number of parameters, while the `weight` of evidence summarizes the relative support for each model.


## Open Science Statement {#sec-M-SM-OS}

In an effort to improve the transparency and replicability of the analysis, this study provides access to an online walk-through. The digital document contains all the code and materials utilized in the study. Furthermore, the walk-through meticulously follows the When-to-Worry-and-How-to-Avoid-the-Misuse-of-Bayesian-Statistics checklist (WAMBS checklist) developed by @Depaoli_et_al_2017. This checklist outlines the ten crucial points that need careful scrutiny when employing Bayesian inference procedures. The digital walk-through is available at the following URL: [https://jriveraespejo.github.io/paper1_manuscript/](https://jriveraespejo.github.io/paper1_manuscript/)


# Results {#sec-results}

This section presents the results of the Bayesian inference procedures, with particular emphasis on answering the three research questions.

```{r}
#| label: code-models

# load reference models
for(i in 1:12){
  model_nam = paste0( ifelse(i<10, 'model0', 'model'), i)
  model_out = file.path( main_dir, 'real_chain')
  model_fit = file_id( model_out, model_nam )
  assign( model_nam,
          rstan::read_stan_csv( file.path( model_out, model_fit ) ) )
}

```


## Predictive capabilities of the Beta-proportion GLLAMM compared to the Normal LMM (RQ1) {#sec-R-RQ1}

This research question evaluated the effectiveness of the Beta-proportion GLLAMM in handling the features of entropy scores by comparing its predictive accuracy to the Normal LMM. Models $1$, $4$, $7$, and $10$ were specifically chosen for this comparison because their assumptions exclusively addressed the features of the scores, without integrating additional covariate information. As detailed in @tbl-fitted, Model $1$ was a Normal LMM that solely addresses data clustering. Building upon this, Model $4$ introduced a robust feature. Conversely, Model $7$ was a Beta-proportion GLLAMM that deals with boundedness, measurement error and data clustering, and Model $10$ extended this model by incorporating a robust feature. 

The left panel of @fig-rq1-waic-psis displays the models' `DIC`, `WAIC`, and `PSIS` values with their corresponding uncertainty intervals. In contrast, the right panel of the figure shows the models' `dWAIC` and `dPSIS` values with their corresponding uncertainty intervals. @tbl-rq1-waic and [-@tbl-rq1-psis] provide similar information, while also reporting the `pWAIC` and `pPSIS` values and the `weight` of evidence for each model. Overall, all criteria consistently pointed to Model $10$ as the most plausible choice for the data. The model exhibits the lowest values for both `WAIC` and `PSIS`, establishing itself as the model with the least deviation from *perfect* predictive accuracy among those under comparison. Additionally, @fig-rq1-waic-psis visually demonstrates the non-overlapping uncertainty in both `dWAIC` and `dPSIS` values for Models $1$, $4$, and $7$ when compared to Model $10$. This indicates that Model $10$ significantly deviated the least from *perfect* predictive accuracy when compared to the rest of the models. Lastly, the `weight` of evidence in @tbl-rq1-waic and [-@tbl-rq1-psis] underscored that $100\%$ of the evidence aligned with and supported Model $10$.

```{r}
#| output: false
#| label: code-rq1-waic

set.seed(12345)

RQ1_WAIC = compare( func=WAIC,
                    model01, model04,
                    model07, model10 )
RQ1_WAIC = cbind( DIC=with(RQ1_WAIC, WAIC-2*pWAIC),
                  RQ1_WAIC )

```


```{r}
#| output: false
#| label: code-rq1-psis

set.seed(12345)

RQ1_PSIS = compare( func=PSIS,
                    model01, model04,
                    model07, model10 )
RQ1_PSIS = cbind( DIC=with(RQ1_PSIS, PSIS-2*pPSIS),
                  RQ1_PSIS )

```


```{r}
#| label: fig-rq1-waic-psis
#| fig-cap: 'Comparison plot for selected models. **_Note:_** Open, black and blue points describe the posterior means for the criteria. Continuous colored horizontal lines indicate the criteria associated uncertainty.'
#| fig-height: 4
#| fig-width: 10

par(mfrow=c(1,2))
plot_compare( waic_obj=RQ1_WAIC, psis_obj=RQ1_PSIS, ns=1, dM=F )
plot_compare( waic_obj=RQ1_WAIC, psis_obj=RQ1_PSIS, ns=1, dM=T )
par(mfrow=c(1,1))

```


Upon closer examination, the reasons behind the observed disparities in the models become more apparent. Specifically, @fig-rq1-pred-speaker demonstrates that the Normal LMM, as outlined in Model $4$, failed to adequately capture the data's underlying patterns, resulting in predictions that were physically inconsistent. This issue is illustrated by the $95\%$ Highest Probability Density Intervals (HPDI) extending beyond the expected zero to one outcome range. Further insight into this lack of fit is provided by @fig-rq1-pred-speaker_model04. The figure displays score prediction densities for Model $4$ that bore no resemblance to the actual data densities. Furthermore, the top two panels in @fig-rq1-model-outliers reveal that misspecification in the Normal LMM causeed the model to be *more surprised* by extreme entropy scores, leading to their identification as highly unlikely and influential observations. Consequently, the model was rendered unreliable due to the potential biases present in the parameter estimates. In contrast, the Beta-proportion GLLAMM appeared to effectively capture the data patterns, generating predictions within the expected data range. This is evident in @fig-rq1-pred-speaker and complemented by @fig-rq1-pred-speaker_model10 and [-@fig-rq1-model-outliers]. In @fig-rq1-pred-speaker_model10, Model $10$ displayed prediction densities that bore more resemblance to the actual data densities. Furthermore, the bottom two panels in @fig-rq1-model-outliers show the model was *less surprised* by extreme scores, fostering more trust in the model's estimates.

```{r}
#| label: fig-rq1-pred-speaker
#| fig-cap: 'Entropy scores prediction for selected models. **_Note:_** Black points show manifest entropy scores where darker points indicate greater overlap. Orange dots and vertical lines show the posterior mean and 95% HPDI derived from Model 4. Blue dots and vertical lines show similar information from Model 10.'
#| fig-height: 6
#| fig-width: 10

plot_speaker(d=data_H,
             stanfit_obj1=model04,
             stanfit_obj2=model10,
             p=0.95,
             decreasing=F,
             leg=c('model 04','model 10'))

```


## Estimation of speakers’ latent potential intelligibility from manifest entropy scores (RQ2) {#sec-R-RQ2}

The second research question aimed to demonstrate the application of the Beta-proportion GLLAMM in estimating the latent potential intelligibility of speakers. This was achieved by employing the general mathematical formalism outlined in @eq-beta-GLLAMM-structural, along with additional specifications provided in @tbl-fitted. The Bayesian procedure successfully estimated the latent potential intelligibility of speakers under Model $10$ through the following structural equation:

$$
SI_{i} = \alpha + e_{i} + u_{i}
$$ {#eq-beta-GLLAMM-structural-model10}

Moreover, due to its implementation under Bayesian procedures, Model $10$ provided the complete posterior distribution of the speakers' potential intelligibility scores. This provision, in turn, (1) enabled the calculation of summaries, facilitating the ranking of individuals, and (2) supported the assessment of differences among selected speakers. In both cases, the model considered the inherent uncertainty of the estimates resulting from its measurement using multiple entropy scores.

```{r}
#| label: code-rq2-si-model10

SI = pred_SI(d=data_H, stanfit_obj=model10, p=0.95) # <1>
SI = SI[order(SI$mean, decreasing=T), ]
SI = round(SI, 3)
SI$HPDI = paste0('[', SI$HPDI_lower, ', ', SI$HPDI_upper,']')

```

```{r}
#| include: false
#| label: tbl-rq2-si-model10
#| tbl-cap: "Latent potential Intelligibility and additional covariates"

kable( SI[,c(1:5,11)], align='rrrrrr', row.names=F, digits=3,
       col.names=c('Speaker ID','Hearing status',
                   'Chron. Age','Chron. Age (centered)',
                   'Posterior mean','95% HPDI') )

```

```{r}
#| label: fig-rq2-si-model10
#| fig-cap: 'Model 10, latent potential intelligibility of speakers. **_Note:_** Black dots and vertical lines show the posterior means and 95% HPDI intervals.'
#| fig-height: 5
#| fig-width: 10

plot_SI(d=data_H, stanfit_obj=model10,
        p=0.95, decreasing=T)

```


@fig-rq2-si-model10 displays the ranking of speakers in decreasing order based on the posterior means of the latent potential intelligibility. These estimates are accompanied by their associated $95\%$ HPDI. The figure indicates that speaker $6$ stands out as the least intelligible in the sample, followed further behind by speaker $1$, $17$ and $9$. In contrast, the figure highlights speaker $20$ as the most intelligible, closely followed by speakers $23$, $31$ and $3$. Conversely, the full posterior distribution for comparing potential intelligibility between the least and most intelligible speakers against other selected speakers is shown in @fig-si_contr_model10. The figure reveals that only the differences between speakers $6$, $1$, $17$, and $9$, along with the difference between speakers $20$ and $3$ are statistically significant, as their associated $95\%$ HPDI did not overlap with zero (shaded area). The R code to derive these scores and generate the figure is available in the digital walk-through document (see @sec-M-SM-OS Open Science Statement).


```{r}
#| label: code-si-contr-model10

SI = contrast_SI(d=data_H, stanfit_obj=model10, 
                 speakers=c(6,20), p=0.95, raw=T)
SI_contr = SI$SI_contr
SI_contr = round(SI_contr, 3)
SI_contr$names = rownames(SI_contr)
SI_contr$HPDI = paste0('[', SI_contr$HPDI_lower, ', ', SI_contr$HPDI_upper,']')

```


```{r}
#| include: false
#| label: tbl-si-contr-model10
#| tbl-cap: "Latent potential intelligibility contrasts of selected speakers"

idx_comp = c(1,21,13,52,60,6)
kable( SI_contr[idx_comp,c('names','mean','HPDI')], row.names=F, align='ccc',
       col.names=c('Contrasts','Posterior mean','95% HPDI') )

```


```{r}
#| label: fig-si_contr_model10
#| fig-cap: 'Model 10, potential intelligibility comparisons among selected speakers. **_Note:_** Shaded area describes the 95% HPDI.'
#| fig-height: 7
#| fig-width: 10

par(mfrow=c(2,3))

for(i in idx_comp){
  dens( SI$SI_raw[[i]], xlim=c(-2.5,2.5),
        col=rgb(0,0,0,0.7), show.HPDI=0.95,
        xlab='Difference in potential intelligibility')
  abline( v=0, lty=2, col=rgb(0,0,0,0.3))
  mtext( text=names(SI$SI_raw)[i], side=3, adj=0, cex=1.1)
}

par(mfrow=c(1,1))

```


## Testing the influence of speaker-related factors on intelligibility (RQ3) {#sec-R-RQ3}

This research question illustrated how hypotheses on intelligibility can be examined within the model's framework. Specifically, the focus centered on assessing the influence of speaker-related factors on intelligibility, such as chronological age and hearing status. Notably, despite RQ1 indicating the suitability of the Beta-proportion GLLAMM for entropy scores, existing statistical literature suggests that, in certain scenarios, models incorporating covariate adjustment exhibit robustness to misspecification in the functional form of the covariate-outcome relationship [@Tackney_et_al_2023]. Consequently, this study compared all models detailed in Table 2. These models were characterized by different covariate adjustments on entropy scores or the latent potential intelligibility of speakers, namely chronological age and hearing status. Furthermore, some models like the Normal LMMs, potentially exhibited misspecification in the covariate-outcome relationship.

Similar to RQ1, all criteria consistently identified the Beta-proportion GLLAMM outlined in models $12$, $11$ and $10$ as the most plausible models for the data. The models exhibited the lowest values for both `WAIC` and `PSIS`, establishing them as the least deviating models among those under comparison. In addition, @fig-rq3-waic-psis depicts the non-overlapping uncertainty for the models' `dWAIC` and `dPSIS` values with horizontal blue lines. This reveals that, when compared to Model $11$, most models exhibited significantly distinct predictive capabilities. Models $12$ and $10$, however, stood out as exceptions to this pattern. This observation suggests that Models $11$, $12$, and $10$ displayed the least deviation from *perfect* predictive accuracy in contrast to the other models. Lastly, the `weight` of evidence in Tables @tbl-rq3-waic and [-@tbl-rq3-psis], underscored that Model $11$ accumulated the greatest support, followed by Model $12$, and lastly, by Model $10$.

```{r}
#| output: false
#| label: code-rq3-waic
#| fig-cap: ''

set.seed(12345)

RQ3_WAIC = compare( func=WAIC,
                    model01, model02, model03,
                    model04, model05, model06,
                    model07, model08, model09,
                    model10, model11, model12 )
RQ3_WAIC = cbind( DIC=with(RQ3_WAIC, WAIC-2*pWAIC),
                  RQ3_WAIC )

```


```{r}
#| output: false
#| label: code-rq3-psis
#| fig-cap: ''

set.seed(12345)

RQ3_PSIS = compare( func=PSIS,
                    model01, model02, model03,
                    model04, model05, model06,
                    model07, model08, model09,
                    model10, model11, model12 )
RQ3_PSIS = cbind( DIC=with(RQ3_PSIS, PSIS-2*pPSIS),
                  RQ3_PSIS )

```


```{r}
#| label: fig-rq3-waic-psis
#| fig-cap: 'Comparison plot for all models. **_Note:_** Open, black and blue points describe the posterior means for the criteria. Continuous colored horizontal lines indicate the criteria associated uncertainty.'
#| fig-height: 8
#| fig-width: 10

par(mfrow=c(1,2))
plot_compare( waic_obj=RQ3_WAIC, psis_obj=RQ3_PSIS, ns=1, dM=F )
plot_compare( waic_obj=RQ3_WAIC, psis_obj=RQ3_PSIS, ns=1, dM=T )
par(mfrow=c(1,1))

```


A closer examination of two models within this comparison set reveals the reasons behind the largest observed disparities. The Normal LMM, as outlined in Model $6$, continued to face challenges in capturing underlying data patterns, resulting in predictions that are physically inconsistent, falling outside the outcome's range. Additionally, the model persisted in identifying highly unlikely and influential observations, making it inherently unreliable. In contrast, the Beta-proportion GLLAMM described by Model $12$ appeared to be less susceptible to extreme scores, effectively capturing data patterns within the expected data range and thereby instilling greater confidence in the reliability of the model's estimates. This contrast is visually depicted in @fig-rq3-pred-speaker, [-@fig-rq3-pred-speaker_model06], [-@fig-rq3-pred-speaker_model12], and [-@fig-rq3-model-outliers].

Considering the results in @fig-rq3-waic-psis, the model comparisons favored three distinct models: Model $10$, $11$ and $12$. Model $10$, supported by $20.4\%$ of the evidence, estimated a single intercept $\alpha$ and no slope to explain the potential intelligibility of speakers (see @tbl-parameter-model10). In contrast, supported by $45.1\%$ of the evidence, Model $11$ in @tbl-parameter-model11 estimated distinct intercepts for each hearing status group, namely $\alpha_{HS[1]}$ for NH speakers and $\alpha_{HS[2]}$ for the HI/CI counterparts, while maintaining a single slope that gauges the impact of age on potential intelligibility estimates. The $95\%$ HPDI for the comparison of intercepts $\alpha_{HS[2]}-\alpha_{HS[1]}$ revealed significant differences between NH and HI/CI speakers. Lastly, with evidence of $34.1\%$, Model $12$ in @tbl-parameter-model12 estimated different intercepts and slopes per hearing status group, namely $\alpha_{HS[1]}$ and $\beta_{A,HS[1]}$ for the NH speakers, and $\alpha_{HS[2]}$ and $\beta_{A,HS[2]}$ for the HI/CI counterparts. The $95\%$ HPDI for the comparison of intercepts and slopes revealed significant differences solely in the slopes between NH and their HI/CI counterparts ($\beta_{A,HS[2]}-\beta_{A,HS[1]}$).

However, a discerning reader can notice that these models yielded conflicting conclusions regarding the influence of chronological age and hearing status on intelligibility. Model $10$ implied no influence of chronological age and hearing status on the potential intelligibility of speakers. @fig-rq3-intelligibility-model10, however, revealed the reason for the model's low support. Model $10$ failed to capture the prevalent increasing age pattern observed in potential intelligibility estimates. In contrast, Model $11$ identified significant differences in potential intelligibility between NH and HI/CI speakers. The model further suggested that with the progression of chronological age, HI/CI speakers lag behind in intelligibility development, with no opportunity to catch up to their NH counterparts within the analyzed age range, as depicted in @fig-rq3-intelligibility-model11. Finally, Model $12$ indicated no significant differences in intelligibility between NH and HI/CI speakers at $68$ months of age (around $6$ years old). However, the model revealed distinct evolution patterns of intelligibility per unit of chronological age between different hearing status groups, with HI/CI speakers displaying a slower rate of development compared to their NH counterparts within the analyzed age range. The latter is evident in @fig-rq3-intelligibility-model12.

```{r}
#| label: code-parameter-model10
#| fig-cap: ''

par_int = c('a','aHS[1]','aHS[2]','bAm','bAmHS[1]','bAmHS[2]')
model_par = par_recovery(stanfit_obj=model10,
                         p=0.95, est_par=par_int)

vars_int = c('mean','HPDI_lower','HPDI_upper')
model_tab = model_par[,vars_int]
model_tab = round(model_tab, 2)
model_tab$names = c('$\\alpha$')
model_tab$HPDI = paste0('[', model_tab$HPDI_lower, ', ', model_tab$HPDI_upper,']')

```


```{r}
#| label: tbl-parameter-model10
#| tbl-cap: 'Model 10, parameter estimates and 95% HPDI.'

vars_tab = c('names','mean','HPDI')
opts = options(knitr.kable.NA="")
kable( model_tab[,vars_tab], row.names=F, align='ccc', 
       escape=F, digits=3,
       col.names=c('Parameter','Posterior mean','95% HPDI') )

```


```{r}
#| label: fig-rq3-intelligibility-model10
#| fig-cap: 'Model 10, Potential intelligibility per chronological age and hearing status. **_Note:_** Colored dots denote the posterior means, vertical lines describe the 95% HPDI, thick discontinuous line indicate the regression line, thin continuous lines denote regression lines samples from the posterior distribution, and numbers indicate the speaker index.'
#| fig-height: 5
#| fig-width: 10

pred_intel(d=data_H, stanfit_obj=model10,
           p=0.95, ns=500, seed=12345)

```


```{r}
#| label: code-parameter-model11
#| fig-cap: ''

par_int = c('a','aHS[1]','aHS[2]','bAm','bAmHS[1]','bAmHS[2]')
model_par = par_recovery(stanfit_obj=model11,
                         p=0.95, est_par=par_int)
contrs = contrast_intel(stanfit_obj=model11, p=0.95,
                        rope=c(-0.05,0.05))

vars_int = c('mean','HPDI_lower','HPDI_upper')
model_tab = rbind(model_par[,vars_int], rep(NA,3), rep(NA,3), contrs[,vars_int])
model_tab = round(model_tab, 2)
model_tab$names = c('$\\alpha$','$\\alpha_{HS[1]}$','$\\alpha_{HS[2]}$',
                    '$\\beta_{A}$', NA,
                    'Contrasts','$\\alpha_{HS[2]} - \\alpha_{HS[1]}$')
model_tab$HPDI = paste0('[', model_tab$HPDI_lower, ', ', model_tab$HPDI_upper,']')
model_tab[5, 1:5] = NA
model_tab[6, c(1:3,5)] = NA

```


```{r}
#| label: tbl-parameter-model11
#| tbl-cap: 'Model 11, parameter estimates and 95% HPDI.'

vars_tab = c('names','mean','HPDI')
opts = options(knitr.kable.NA = "")
kable( model_tab[,vars_tab], row.names=F, align='ccc', 
       escape=F, digits=3,
       col.names=c('Parameter','Posterior mean','95% HPDI') )

```


```{r}
#| label: fig-rq3-intelligibility-model11
#| fig-cap: 'Model 11, Potential intelligibility per chronological age and hearing status. **_Note:_** Colored dots denote the posterior means, vertical lines describe the 95% HPDI, thick discontinuous line indicate the regression line, thin continuous lines denote regression lines samples from the posterior distribution, and numbers indicate the speaker index.'
#| fig-height: 5
#| fig-width: 10

pred_intel(d=data_H, stanfit_obj=model11,
           p=0.95, ns=500, seed=12345)

```


```{r}
#| label: code-parameter-model12
#| fig-cap: ''

par_int = c('a','aHS[1]','aHS[2]','bAm','bAmHS[1]','bAmHS[2]')
model_par = par_recovery(stanfit_obj=model12,
                         p=0.95, est_par=par_int)
contrs = contrast_intel(stanfit_obj=model12, p=0.95,
                        rope=c(-0.05,0.05))

vars_int = c('mean','HPDI_lower','HPDI_upper')
model_tab = rbind(model_par[,vars_int], rep(NA,3), rep(NA,3), contrs[,vars_int])
model_tab = round(model_tab, 2)
model_tab$names = c('$\\alpha$','$\\alpha_{HS[1]}$','$\\alpha_{HS[2]}$',
                    '$\\beta_{A,HS[1]}$','$\\beta_{A,HS[2]}$', NA,
                    'Contrasts','$\\alpha_{HS[2]} - \\alpha_{HS[1]}$',
                    '$\\beta_{A,HS[2]} - \\beta_{A,HS[1]}$')
model_tab$HPDI = paste0('[', model_tab$HPDI_lower, ', ', model_tab$HPDI_upper,']')
model_tab[6, 1:5] = NA
model_tab[7, c(1:3,5)] = NA

```


```{r}
#| label: tbl-parameter-model12
#| tbl-cap: 'Model 12, parameter estimates and 95% HPDI.'

vars_tab = c('names','mean','HPDI')
opts = options(knitr.kable.NA = "")
kable( model_tab[,vars_tab], row.names=F, align='ccc', 
       escape=F, digits=3,
       col.names=c('Parameter','Posterior mean','95% HPDI') )

```


```{r}
#| label: fig-rq3-intelligibility-model12
#| fig-cap: 'Model 12, Potential intelligibility per chronological age and hearing status. **_Note:_** Colored dots denote the posterior means, vertical lines describe the 95% HPDI, thick discontinuous line indicate the regression line, thin continuous lines denote regression lines samples from the posterior distribution, and numbers indicate the speaker index.'
#| fig-height: 5
#| fig-width: 10

pred_intel(d=data_H, stanfit_obj=model12,
           p=0.95, ns=500, seed=12345)

```


# Discussion {#sec-discussion}

## Findings {#sec-D-F}

This study examined the suitability of the Bayesian Beta-proportion GLLAMM for the quantitative measuring and testing of research hypotheses related to speech intelligibility using entropy scores. The initial findings supported the assertion that Beta-proportion GLLAMMs consistently outperformed Normal LMMs in predicting entropy scores, underscoring its superior predictive performance. The results also emphasized that models neglecting measurement error and boundedness in the outcomes lead to underfitting and misspecification issues, even when robust features are integrated. This was clearly illustrated by the Normal LMMs.

Secondly, the study showcased the Beta-proportion GLLAMM's proficiency in estimating the latent potential intelligibility of speakers based on manifest entropy scores. Implemented under Bayesian procedures, the proposed model offered a valuable advantage over frequentist methods by further providing the full posterior distribution of the speakers' potential intelligibility. This provision facilitated the calculation of summaries, aiding in the construction of individual rankings, and supported the comparisons among selected speakers. In both scenarios, the proposed model accounted for the inherent uncertainty in the intelligibility estimates.

Thirdly, the study illustrated how the proposed model assessed the impact of speaker-related factors on potential intelligibility. The results suggested that multiple models were plausible for the observed entropy scores. This indicated that different speaker-related factor hypotheses were viable for the data, with some presenting contradictory conclusions about the influence of these factors on intelligibility. However, even without unequivocal support for one hypothesis, the divided support among these models informed that certain statistical issues may be hindering the models' ability to distinguish among individuals and, ultimately, among models. These issues may be attributed to factors such as the insufficient sample size of speakers, the inadequate representation of the population of speakers, referred to as selection bias, and the imprecise measurement of the latent variable of interest.

Ultimately, this study introduced researchers to innovative statistical tools that enhanced existing research models. These tools not only assessed the predictability of empirical phenomena but also quantitatively measured the latent trait of interest, namely potential intelligibility, facilitating the comparison of research hypotheses related to this trait. However, the presented tools introduce new challenges for researchers seeking their implementation. These challenges emerge from two distinct aspects: one methodological and the other practical. In the methodological domain, researchers need familiarity with Bayesian methods and the principled formulation of assumptions regarding the data-generating process and research inquiries. This entails understanding and addressing each of the data and research challenges within the context of a statistical (probabilistic) models. Conversely, in the practical domain, researchers need familiarity with probabilistic programming languages (PPLs), which are designed for specifying and obtaining inferences from probabilistic models -the core of Bayesian methods. To ensure the successful utilization of this new statistical tool, this study addressed both challenges by providing comprehensive, step-by-step guidance in the form of a digital walk-through document (see @sec-M-SM-OS Open Science Statement).


## Limitations and further research {#sec-D-LFR}

This study provided valuable insights into the use of a novel approach to simultaneously address the different data features of entropy scores in speech intelligibility research. However, it is important to acknowledge the limitations of this study and explore potential avenues for future research. Firstly, the study interpreted potential intelligibility as an unobserved latent trait of speakers influencing the likelihood of observing a set of entropy scores. These scores, in turn, reflected the transcribers' ability to decode words in sentences produced by the speakers. Despite this practical approach, the construct validity of the latent trait heavily depended on the listeners' appropriate understanding and execution of the transcription task. Construct validity, as defined by @Cronbach_et_al_1955, refers to the extent to which a set of manifest variables accurately represents a concept that cannot be directly measured. Considering the study assumed the transcription task set by @Boonen_et_al_2023 was properly understood and executed, it expected that the potential intelligibility reflected the overall speech intelligibility of speakers. However, the study did not delved into the general epistemological considerations regarding the connection between the latent variable and the concept.

<!-- These differences between transcriptions may be used in a more refined calculation of entropy. In the present study, each deviance of the transcriptions was equally weighed. In other words, each difference equally increased the entropy score. However, some deviances in the transcriptions are fairly small (e.g., kikker ‘frog’ vs. kikkers ‘frogs’), whereas others can really be considered as mismatches (e.g., jongen ‘boy’ vs. hond ‘dog’). Further research is needed for finding fruitful ways to refine the measure by taking into account the (linguistic) distance between different transcriptions. -->

Secondly, the study revealed a notable lack of unequivocal support for one of the models among the compared set. This outcome may be attributed to factors such as the insufficient sample size of speakers, the inadequate representation of the populations of speakers (referred to as selection bias), and the imprecise measurement of the latent variable. Small sample sizes and selection bias yield data with limited outcome and covariates ranges, leading to biased and imprecise parameter estimates [@Everitt_et_al_2010]. Moreover, fueled by the reduced measurement precision, these issues can result in models with diminished statistical power and a higher risk of type I or type II errors [@McElreath_2020]. Consequently, future research should consider extending this study by conducting formal sample size planning. This entails assessing the impact of expanding the speakers' pool on testing research hypotheses or increasing the number of speech samples, transcriptions, and listeners to enhance the precision of the potential intelligibility estimates. With these insights, future investigations could contemplate increasing the speaker sample with a group that adequately represents the population of interest. However, this must be done while mindful of the pragmatic limitations associated with transcription tasks, specifically considering the costs and time-intensiveness of the procedure.

Thirdly, the study presented an illustrative example for the investigation of research hypotheses within the model's framework. However, it did not offer an exhaustive evaluation of all factors influencing intelligibility, which are thoroughly explored in the works of @Niparko_et_al_2010, @Boons_et_al_2012, @Gillis_2018, and @Fagan_et_al_2020. Consequently, the study could not discard the presence of unobservable variables that might bias the parameter estimates, potentially impacting the inferences provided. Hence, future research should consider integrating appropriate causal hypotheses about these factors into the proposed models, as proper covariate adjustment facilitates the production of unbiased and precise parameter estimates [@Cinelli_et_al_2021; @Deffner_et_al_2022].

Lastly, this study proposes two directions for future exploration in speech intelligibility research. Firstly, there is an opportunity to investigate alternative methods for assessing speech intelligibility beyond transcription tasks and entropy scores. The experimental design of transcription tasks imply that the procedure may be time-intensive and costly. Thus, exploring less time-intensive or more cost-effective procedures, that still offer comparable precision in intelligibility estimates, could benefit both researchers and speech therapists alike. One example of such a method is Comparative Judgment (CJ), where judges compare and score the perceived intensity of a trait between two stimuli [@Thurstone_1927]. CJ has gained increasing attention in educational assessment, with several studies demonstrating its validity in assessing various tasks within student work, as shown in @Pollitt_2012a, @Pollitt_2012b, @Lesterhuis_2018, @vanDaal_2020, and @Verhavert_et_al_2019. The work of @Boonen_et_al_2020 illustrates the potential of this methodology to assess intelligibility. In their study, the authors assessed the overall speech quality of hearing-impaired children using pairwise comparisons of uttered speech samples, while scoring the results in a dichotomous manner. Nevertheless, there is significant room for extending their application. For instance, researchers can perform retrospective power analysis to ascertain the power of the study's claims [see @Kruschke_2015, pp. 393-394]. Furthermore, the application can be extended to other unexplored variants of the CJ method, such as Ordered CJ [@Pritikin_2020] or Multidimensional Dichotomous CJ.

Conversely, a second avenue for exploration involves integrating diverse data types and evaluation methods to assess individuals' intelligibility. This can be accomplished by leveraging two features of Bayesian methods: their flexibility and the concept of Bayesian updating. Bayesian methods possess the flexibility to simultaneously handle various data types. Additionally, through Bayesian updating, researchers can integrate information from the posterior distribution of parameters as priors in models for subsequent evaluations. Ultimately, this could enable researchers to assess speakers' intelligibility progress without committing to a specific data type or evaluation method. This advancement could mirror the emergence of second-generation Structural Equation Models proposed by @Muthen_2001, where models facilitate the combined estimation of categorical and continuous latent variables. However, in the context of future research, the proposal would facilitate the estimation of latent variables using a combination of data types and evaluation methods, contingent upon the fulfillment of construct validity by those evaluation methods.


# Conclusion {#sec-conclusion}

This study have highlighted the effectiveness of the Bayesian Beta-proportion GLLAMM to collectively address several key data features when investigating unobservable and complex traits. The study used speech intelligibility and entropy scores as a motivating example. The results have demonstrated that the proposed model consistently outperforms the Normal LMM in predicting the empirical phenomena. Moreover, the model exhibits the ability to quantify the latent potential intelligibility of speakers, allowing for the ranking and comparison of individuals based on the latent trait while accommodating associated uncertainties. Additionally, the proposed model have facilitated the exploration of research hypotheses concerning the influence of speaker-related factors on potential intelligibility, where the integration and comparison of these hypotheses within the model's framework was a straightforward task.

However, the introduction of these innovative statistical tools presents new challenges for researchers seeking implementation. These challenges encompass the principled formulation of assumptions about the data-generating processes and research inquiries, along with the need for familiarity with probabilistic programming languages (PPLs) essential for implementing Bayesian methods. Nevertheless, the study suggests several promising avenues for future research, including causal hypothesis formulation, and the exploration and integration of novel evaluation methods for assessing intelligibility. The insights derived from this study hold implications for both researchers and data analysts interested in quantitatively measuring intricate, unobservable constructs, while predicting accurately the empirical phenomena.



{{< pagebreak >}}

# Declarations {.unnumbered}

**Funding:** The project was founded through the Research Fund of the University of Antwerp (BOF).

**Conflict of interests:** The authors declare no conflict of interest.

**Ethics approval:** This is an observational study. The University of Antwerp Research Ethics Committee has confirmed that no ethical approval is required.

**Consent to participate:** Not applicable

**Consent for publication:** All authors have read and agreed to the published version of the manuscript.

**Availability of data and materials:** The data is delivered upon request. 

**Code availability:** All the code utilized in this research is available in the different notebooks and `CODE LINKS` referenced in the digital document. The digital document is located at: [https://jriveraespejo.github.io/paper1_manuscript/](https://jriveraespejo.github.io/paper1_manuscript/).  

**Authors' contributions:** *Conceptualization:* S.G., S.dM., and J.M.R.E; *Data curation:* J.M.R.E.; *Formal Analysis:* J.M.R.E.; *Funding acquisition:* S.G. and S.dM; *Investigation:* S.G.; *Methodology:* S.G., S.dM., and J.M.R.E; *Project administration:* S.G. and S.dM.; *Resources:* S.G. and S.dM.; *Software: J.M.R.E.*; *Supervision:* S.G. and S.dM.; *Validation:* J.M.R.E.; *Visualization:* J.M.R.E.; *Writing - original draft:* J.M.R.E.; *Writing - review & editing:* S.G. and S.dM.

**Acknowledgements:** We express gratitude to Renaat van Uffelen for providing timely expertise in enhancing the writing of this document.


{{< pagebreak >}}

# Appendix {#sec-appendix}

## Entropy scores calculation {#sec-appA}

This section exemplified the entropy calculation procedure. For that purpose, the words in position two, three, four and five observed in @tbl-alignment were used. These words were assumed present in the first sentence, produced by the first speaker assigned to the first block, and transcribed by five listeners ($w=\{2,3,4,5\}$, $s=1$, $i=1$, $b=1$, $J=5$). For second word, the first four listeners identified the word type *jongen* $(T_{j1})$, while the last identified the word type *hond* $(T_{j2})$. Therefore, two word types were identified ($K=2$), with proportions equal to $\{ p_{1}, p_{2} \}$ $=$ $\{ 4/5, 1/5 \}$ $=$ $\{ 0.8, 0.2 \}$, with an entropy score equal to:

$$ 
H_{2111} = \frac{ -\left[ 0.8 \cdot log_{2}(0.8) + 0.2 \cdot log_{2}(0.2) \right] }{ log_{2}(5)} \approx 0.3109
$$

For the fourth word, two listeners identified the word type *een* $(T_{j1})$, one listener the word type *de* $(T_{j2})$, and another the word *geen* $(T_{j3})$. In addition, a blank space *[B]* is a symbol that defined the absence of a word in a space where a word was expected during the alignment procedure, as compared with other transcriptions. Notice that for calculation purposes, because the blank space was not expected in such position, it was considered as a different word type. Consequently four word types were registered ($K=4$), with proportions equal to $\{ p_{1}, p_{2}, p_{3}, p_{4} \}$ $=$ $\{ 2/5, 1/5, 1/5, 1/5 \}$ $=$ $\{ 0.4, 0.2, 0.2, 0.2 \}$ with an entropy score equal to:

$$ 
H_{4111} = \frac{ -\left[ 0.4 \cdot log_{2}(0.4) + 3 \cdot 0.2 \cdot log_{2}(0.2) \right] }{ log_{2}(5)} \approx 0.8277
$$
For the fifth word, each listener transcribed a different word. It is important to highlight that when a listener did not identify a complete word, or part of it, (s)he was instructed to write *[X]* in that position. However, for the calculation of the entropy score, if more than one listener marked an unidentifiable word with *[X]*, each one of them was considered a different word type. This was done to avoid the artificial reduction of the entropy score, as *[X]* values already indicated the word’s lack of intelligibility. . Consequently, five word types were observed, $T_{j1}=$*kikker*, $T_{j2}=$*[X]*, $T_{j3}=$*kokkin*, $T_{j4}=$*kikkers*, $T_{j5}=$*[X]* ($K=5$), with proportions equal to $\{ p_{1}, p_{2}, p_{3}, p_{4}, p_{5} \}$ $=$ $\{ 1/5, 1/5, 1/5, 1/5, 1/5 \}$ $=$ $\{ 0.2, 0.2, 0.2, 0.2, 0.2 \}$, with an entropy score equal to:

$$ 
H_{5111} = \frac{ -\left[ 5 \cdot 0.2 \cdot log_{2}(0.2) \right] }{ log_{2}(5)} = 1
$$

Lastly, for the third word, the first two listeners identified the word type *ziet* $(T_{j1})$, the next two listeners identified the word type *zag* $(T_{j2})$, while the last one identified the word type *zoekt* $(T_{j3})$. Consequently, three word types were identified ($K=3$), with proportions equal to $\{ p_{1}, p_{2}, p_{3} \}$ $=$ $\{ 2/5, 2/5, 1/5 \}$ $=$ $\{ 0.4, 0.4, 0.2 \}$, with an entropy score equal to:

$$
H_{2111} = \frac{ -\left[ 2 \cdot 0.4 \cdot log_{2}(0.4) + 0.2 \cdot log_{2}(0.2) \right] }{ log_{2}(5)} \approx 0.6555
$$

Importantly, the last example showcased the major difference between entropy and measures of accuracy based on the percentage of (un)intelligible words. Entropy scores employ all word type proportions in their calculations, effectively capturing the agreement and disagreement among listeners’ word transcriptions [@Boonen_et_al_2023]. In contrast, the percentage of (un)intelligible words discards most word type proportions in favor of *simpler* agreement or disagreement percentages. For example, an agreement percentage could be reflected by the proportion of the most frequent word, i.e., $\text{max}\{ 0.4, 0.4, 0.2 \}=0.4$, or by other similar percentages detailed in the works of @Flipsen_2006 and @Lagerberg_et_al_2014.


{{< pagebreak >}}

## Tables {#sec-appB}

```{r}
#| label: tbl-rq1-waic
#| tbl-cap: 'WAIC comparison for selected models. **_Note:_** The table is sorted based on `weight` from most to least plausible model(s) for the data.'

RQ1_WAIC$Model = as.integer( str_sub( rownames(RQ1_WAIC), start=-2) )

kable( round( RQ1_WAIC[,c(8,1:7)], 3), row.names=F, 
       align='crrrrrrr', escape=F, digits=2 )

```


```{r}
#| label: tbl-rq1-psis
#| tbl-cap: 'PSIS comparison for selected models. **_Note:_** The table is sorted based on `weight` from most to least plausible model(s) for the data.'

RQ1_PSIS$Model = as.integer( str_sub( rownames(RQ1_PSIS), start=-2) )

kable( round( RQ1_PSIS[,c(8,1:7)], 3), row.names=F, 
       align='crrrrrrr', escape=F, digits=2 )

```


```{r}
#| label: tbl-rq3-waic
#| tbl-cap: 'WAIC comparison for all models. **_Note:_** The table is sorted based on `weight` from most to least plausible model(s) for the data.'

RQ3_WAIC$Model = as.integer( str_sub( rownames(RQ3_WAIC), start=-2) )

kable( round( RQ3_WAIC[,c(8,1:7)], 3), row.names=F, 
       align='crrrrrrr', escape=F, digits=2 )

```


```{r}
#| label: tbl-rq3-psis
#| tbl-cap: 'PSIS comparison for all models. **_Note:_** The table is sorted based on `weight` from most to least plausible model(s) for the data.'

RQ3_PSIS$Model = as.integer( str_sub( rownames(RQ3_PSIS), start=-2) )

kable( round( RQ3_PSIS[,c(8,1:7)], 3), row.names=F, 
       align='crrrrrrr', escape=F, digits=2 )

```


{{< pagebreak >}}

## Figures {#sec-appC}

```{r}
#| label: fig-rq1-pred-speaker_model04
#| fig-cap: 'Model 4: Entropy scores density for selected speakers. **_Note:_** Black bars denote the true data density, orange bars describe the predicted data density'
#| fig-height: 7
#| fig-width: 10

col_string = rep( rethink_palette[2], 2)
pred_speaker_pairs(speakers=c(20,8,11, 25,30,6),
                   d=data_H, stanfit_obj=model04,
                   p=0.95, nbins=20, col_string=col_string)

```


```{r}
#| label: fig-rq1-pred-speaker_model10
#| fig-cap: 'Model 10: Entropy scores density for selected speakers. **_Note:_** Black bars denote the true data density, blue bars describe the predicted data density'
#| fig-height: 7
#| fig-width: 10

col_string = rep( rethink_palette[1], 2)
pred_speaker_pairs(speakers=c(20,8,11, 25,30,6),
                   d=data_H, stanfit_obj=model10,
                   p=0.95, nbins=20, col_string=col_string)

```


```{r}
#| label: fig-rq1-model-outliers
#| fig-cap: 'Outlier identification and analysis for selected models. **_Note:_** Thin and thick vertical discontinuous line indicate threshold of 0.5 and 0.7, respectively. Number pair texts indicate the observation pair of speaker and sentence index. '
#| fig-height: 7
#| fig-width: 10

par(mfrow=c(2,2))
plot_outlier(d=data_H, stanfit_obj=model01)
plot_outlier(d=data_H, stanfit_obj=model04)
plot_outlier(d=data_H, stanfit_obj=model07)
plot_outlier(d=data_H, stanfit_obj=model10)
par(mfrow=c(1,1))

```


```{r}
#| label: fig-rq3-pred-speaker
#| fig-cap: 'Entropy scores prediction for selected models. **_Note:_** Black points show manifest entropy scores where darker points indicate greater overlap. Orange dots and vertical lines show the posterior mean and 95% HPDI derived from Model 6. Blue dots and vertical lines show similar information from Model 12.'
#| fig-height: 6
#| fig-width: 10

plot_speaker(d=data_H,
             stanfit_obj1=model06,
             stanfit_obj2=model12,
             p=0.95,
             decreasing=F,
             leg=c('model 06','model 12'))

```


```{r}
#| label: fig-rq3-pred-speaker_model06
#| fig-cap: 'Model 6: Entropy scores density for selected speakers. **_Note:_** Black bars denote the true data density, orange bars describe the predicted data density'
#| fig-height: 7
#| fig-width: 10

col_string = rep( rethink_palette[2], 2)
pred_speaker_pairs(speakers=c(20,8,11, 25,30,6),
                   d=data_H, stanfit_obj=model06,
                   p=0.95, nbins=20, col_string=col_string)

```


```{r}
#| label: fig-rq3-pred-speaker_model12
#| fig-cap: 'Model 12: Entropy scores density for selected speakers. **_Note:_** Black bars denote the true data density, blue bars describe the predicted data density'
#| fig-height: 7
#| fig-width: 10

col_string = rep( rethink_palette[1], 2)
pred_speaker_pairs(speakers=c(20,8,11, 25,30,6),
                   d=data_H, stanfit_obj=model12,
                   p=0.95, nbins=20, col_string=col_string)

```


```{r}
#| label: fig-rq3-model-outliers
#| fig-cap: 'Outlier identification and analysis for selected models. **_Note:_** Thin and thick vertical discontinuous line indicate threshold of 0.5 and 0.7, respectively. Number pair texts indicate the observation pair of speaker and sentence index.'
#| fig-height: 7
#| fig-width: 10

par(mfrow=c(2,2))
plot_outlier(d=data_H, stanfit_obj=model05)
plot_outlier(d=data_H, stanfit_obj=model06)
plot_outlier(d=data_H, stanfit_obj=model11)
plot_outlier(d=data_H, stanfit_obj=model12)
par(mfrow=c(1,1))

```


{{< pagebreak >}}

# References {.unnumbered}

:::{#refs}

:::