Multidimensional beta-binomial regression model: A joint analysis of patient-reported outcomes

Abstract

Patient-reported outcomes (PROs) are often used as primary outcomes in clinical research studies. PROs are usually measured in ordinal scales and they tend to have excess variability beyond the binomial distribution, a property called overdispersion. Beta-binomial distribution has been previously proposed in this context in order to fit PROs, and beta-binomial regression (BBR) as a good alternative for modelling purposes, including the extension to mixed-effects models in a longitudinal framework. Many PROs have various health dimensions, which are commonly correlated within subjects. However, in clinical analysis, dimensions are separately analysed. In this work, we propose a multidimensional BBR model that incorporates a multidimensional outcome including several PROs in a joint analysis. The proposal has been evaluated and compared to the independent analysis through a simulation study and a real data application with patients with respiratory disease. Results show the advantages that a multidimensional approach offers in terms of parameter significance and interpretation. Additionally, the methods proposed in this work are implemented in the PROreg R-package developed by the authors.

Keywords

Beta-binomial regression patient reported outcomes random-effects multidimensional joint analysis COPD

1 Introduction

Patient-reported outcomes (PROs) are increasingly used as primary outcomes in health services research studies (Weldring and Smith, 2013). For instance, health-related quality of life (HRQoL) measurements are routinely used in both experimental and observational studies when assessing the health status and evolution of patients with chronic diseases (GBD 2016 Disease and Injury Incidence and Prevalence Collaborators, 2018). HRQoL scores are often measured on an ordinal scale. They are built as a sum of ordinal responses to several items collected in a questionnaire form.

The ordinal scaling given by such measures indicates that they tend to generate discrete, bounded and skewed distributed data, often showing U, J or inverse J shapes (Arostegui et al., 2007). In fact, the different items are answered by the same individuals which defines a correlation structure among the ordinal responses that constitute the final score, which increments the variability of the outcome beyond the mean-variance structure of the binomial distribution, a property called overdispersion. In this sense, due to its flexible distributional shape and capacity to account for extra variability, the beta-binomial distribution (BB) has been proposed in the literature to fit PROs (Arostegui et al., 2012), and the beta-binomial regression (BBR) as a good alternative for modelling purposes (Najera-Zuloaga et al., 2018).

Several outcome measures are often collapsed in the same questionnaire. Indeed, most PROs are multidimensional, with the same questionnaire covering different health aspects. For instance, HRQoL is a concept that covers physical, emotional and social components (Testa and Simonson, 1996). In general, HRQoL questionnaires, both generic and disease-specific, are designed to measure different components of health status. Thus, the most common HRQoL questionnaires have various health dimensions or outcomes. The methodology usually used in the design and validation of questionnaires tries to guarantee that each dimension is reflected in a scale that measures a different health concept and that the overlap between these scales is minimal, which is called discriminant validity. Fayers and Machin (2007) state that scale-to-scale correlation in the same questionnaire must be lower than the correlation of each scale with convergent measurements from other questionnaires and, it must be also lower than the reliability measurement of each scale. Nevertheless, in practice, there is always some overlap and consequently, the dimensions of the same questionnaire tend to be correlated (de Vet et al., 2011).

When HRQoL scores have been analyzed as outcome variables in the literature, they have been independently considered as one response variable at each time. Although the interest in a multidimensional approach to analyzing jointly all the dimensions of the same questionnaire has been previously stated (Chang and Reeve, 2005), as far as our knowledge, there is no practical proposal in this context. In our opinion, accounting for the correlation between dimensions, due to within-subject correlation, would better account for the uncertainty in the estimation of the effect of covariates in HRQoL and, consequently, provide a more precise statistical inference.

The extension of BBR in a mixed-effects model framework was proposed recently for analysing longitudinal discrete and bounded outcomes (Najera-Zuloaga et al., 2019). The beta-binomial mixed-effects model (BBMM) is an extension of the BBR model to any hierarchical structure, such as longitudinal or clustered data. However, in practice, the proposal was only limited to longitudinal data, where temporal-specific matrices were defined and applied in both simulation studies and real applications. In this work, we extend the proposed BBMM and define the multidimensional beta-binomial regression model which, based on a mixed effects modelling approach, allows the joint analysis of several outcome variables in a regression setting. Additionally, we updated the PROreg R package to version 1.2, including the implementation of the proposed methodology.

This methodological proposal was originally motivated by a clinical study consisting of 543 patients with Chronic Obstructive Pulmonary Disease (COPD). COPD is considered as a complex, heterogeneous and multisystem disease, hence it is necessary to have tools that enable us to evaluate all those aspects in a broad and integrative approach. HRQoL is considered an important outcome itself or a predictor of other outcomes, such as mortality in COPD patients. In the study, disease-specific HRQoL was assessed using the Spanish validated version of the St. George Respiratory Questionnaire (SGRQ), which provides values bounded in a 0–100 scale where higher scores represent lower health status (Jones et al., 1991). Section S0 in the Supplementary Material gives more detailed information about the patient’s scores recoding process and shows a brief description of the subjects’ scoring patterns, their HRQoL and the most important considered covariates. One of the main aims of the study was to describe the HRQoL of the patients with COPD when they were stable and to test for the impact of several patients’ and disease characteristics on HRQoL, which was assessed by the three dimensions/outcomes that the SGRQ provides, namely activity, impacts and symptoms. Indeed, clinical researchers were interested in a joint multidimensional analysis of the three dimensions, which shows that there is a substantive interest in a joint multidimensional approach in this context, not only from a statistical point of view but also for clinical research.

Accordingly, the goal of this work is two-fold: first, we propose a multidimensional beta-binomial regression model to incorporate several dimensions of the PRO in a joint analysis as multivariate outcomes and, second, we implement the methodology in the PROreg R-package through the BBmm function to provide tools for the clinicians to apply it. The proposal is evaluated in a simulation study including different settings. Moreover, the proposed methodology is applied to patients with COPD to detect and quantify the effect of some factors in disease-specific HRQoL. Results from the multidimensional BBR approach were compared to those obtained from the unidimensional analysis of each dimension at one time, showing a more appropriate interpretation of the fixed effects and a significant improvement in terms of capacity to detect statistically significant covariates.

In the following section, we introduce the details of the multidimensional beta-binomial regression model, including the shared random-effects approach for the joint estimation. We continue with a simulation study where we analyse the performance of our proposal under different scenarios and remark on the main differences compared to the unidimensional approach. Next, we show the analysis of the COPD study, where clinically significant results are presented and compared to the unidimensional models results. Finally, we discuss the results and provide some conclusions and final remarks.

2 A beta-binomial model for multivariate outcomes

In this section, we define the multidimensional regression approach for the joint analysis of different health dimensions provided by the same subject and drawn from beta-binomial distributions. The multivariate outcome is considered through the inclusion of a common random effect among the individuals, the so-called shared random effect approach (McCulloch, 2008).

For an univariate analysis of the SGRQ dimensions with covariates, a beta-binomial regression was first proposed in Williams (1982) applying a logit link function to the probability parameter of the beta-binomial distribution which is connected to the set of covariates of interest.

Details about the definition of the beta-binomial regression are shown in Section S0 of the Supplementary Material. Najera-Zuloaga et al. (2019) developed a beta-binomial mixed-effects model (BBMM) and an estimation procedure extending the proposal by Williams (1982), including normally distributed random-effects in the linear predictor of the beta-binomial regression model.

Let $Y = (Y_{1}, \dots, Y_{n})^{'}$ be the vector of the outcome variables, where conditional on the random-effects $u$ the elements of $Y$ are independent and drawn from a beta-binomial distribution. Furthermore, assume that the random-effects $u$ follow a multivariate normal distribution with zero mean and variance–covariance matrix $D$ , which depends on a variance parameter vector $λ$ and it is non-singular. Namely, the model assumptions are defined as

Y_{i} | u \sim BB (m, p_{i}, ϕ) independentand u \sim N (0, D (λ)), i = 1, \dots, n .

where $p_{i}$ is the probability of obtaining one more point in the score of the outcome, $ϕ$ is the dispersion parameter and $m$ is the maximum of the score.

With the purpose of simplyfing the notation we will define as $θ = (ϕ, λ)^{'}$ the vector of variance components. If we apply a logit link function to the conditional expectation of the model and link it with the given $X_{1}, \dots, X_{t}$ covariates, we have that

η_{i} = \log (\frac{p_{i}}{1 - p_{i}}) = x'_{i} β + z'_{i} u, i = 1, \dots, n,

(2.1)

where $β = (β_{0}, β_{1}, \dots, β_{t})^{'}$ is the vector of the coefficients of the fixed effects, $x'_{i}$ is the $i$ th row of a full rank $n \times (t + 1)$ matrix $X$ composed by the given $t$ covariates, $u$ is a $q$ length vector composed by the random-effects and $z'_{i}$ is the $i$ th row of a $n \times q$ model matrix $Z$ defined by the random structure of the model.

In order to obtain the estimates of the fixed and the random-effects, a first-order Laplace approximation, which is equivalent to integrating the random effects out (Lee and Nelder, 2001), is applied to the marginal log-likelihood function of the model, decomposing it into the joint log-likelihood and an adjusted term. For the estimation of the dispersion parameters, as the maximum likelihood estimation might be substantially biased due to the previous estimation of the fixed-effects, an adjusted profile h-likelihood approach is used (Lee and Nelder, 1996). The general estimation procedure of the BBMM iterates between estimates obtained from the previously mentioned adjusted log-likelihoods until convergence. As it was shown in Najera-Zuloaga et al. (2019), this approach outperforms other beta-binomial implementations such as the one implemented in the gamlss R-package (Rigby and Stasinopoulos, 2005) with the beta-binomial distribution family (i.e., family=BB) due to the fact that in this approach the dispersion parameter is penalized.

More details on the estimation and inference procedure are given in Najera-Zuloaga et al. (2019).

2.1 A shared random-effects approach

Following the proposal by McCulloch (2008), we consider a shared random-effects approach in the BBMM context to account for a multidimensional outcome. The inclusion of a common random effect in the linear predictor of each dimension allows the construction of a correlation structure among them. The basic idea is to use a random effect to account for the correlation within subjects between the dimensions, assuming that conditional on the random effects $u$ the dimensions are independent. The conditionally independent assumption reflects the belief that a common set $u$ of underlying characteristics of the individual governs the outcome process.

For instance, consider that we have two vectors of beta-binomial random variables $Y^{1}$ and $Y^{2}$ associated with the measurements of the same $n$ subjects for two different dimensions, and assume that we are interested in modelling them as a function of a set of covariates $X$ . For the sake of clarity, we only consider two dimensions and one covariate although an extension to more dimensions and covariates is immediate. If we analysed the dimensions separately through BBR models, we would not accommodate the correlation between $Y_{i}^{1}$ and $Y_{i}^{2}$ , $i = 1, \dots, n$ . Therefore, we would have some loss of information in the analysis procedure as covariate mean effects could be masked by subject-specific effects. To overcome this limitation, we propose to introduce a random effect per individual that will be shared by both dimensions. Namely, we propose the following model,

\begin{matrix} Y_{i}^{1} | u_{i} \sim BB (m_{1}, p_{i}^{1}, ϕ_{1}) indp . i = 1, \dots, n \\ logit (p_{i}^{1}) = γ_{0}^{1} + γ_{1}^{1} x_{i} + u_{i}, \\ Y_{i}^{2} | u_{i} \sim BB (m_{2}, p_{i}^{2}, ϕ_{2}) indp . i = 1, \dots, n \\ logit (p_{i}^{2}) = γ_{0}^{2} + γ_{1}^{2} x_{i} + u_{i}, \\ u_{i} \sim N (0, σ_{u}^{2}) . \end{matrix}

(2.2)

The inclusion of a shared random-effect per individual sets up a specific correlation structure in the data. It is well-known that, in non-linear regression models, the specification of the linear predictor does not allow the inclusion of any error term, which does not allow the specification of a variance–covariance structure of the residuals, and hence, of the dimensions involved in the analysis. In this context, a latent variable approach results of main interest as it allows estimating the variability or correlation that exists between the errors of underlying variables associated with each observation. Namely, we assume that there is an unobservable continuous latent random variable $Y_{i}^{l *}$ defined in $ℝ$ for each binary variable $Y_{i j}^{l}$ that constructs the final observed beta-binomial score $Y_{i}^{l}$ , $i = 1, \dots, n, j = 1, \dots, m, l = 1,2$ , such that the observed binary variable $Y_{i j}^{l}$ takes the value 1 if and only if $Y_{i}^{l *}$ exceeds a certain threshold $ζ$ , which without loss of generality can be assumed as 0. Hence,

p_{i}^{l} = P (Y_{i j}^{l} = 1) = P (Y_{i}^{l *} > ζ) = P (Y_{i}^{l *} > 0), i = 1, \dots, n, j = 1, \dots, m, l = 1,2.

The fact that an outcome variable is modelled by a set of covariates implies that the latent variable of the outcome can also be modelled by the same covariates. Therefore, in order to model the dependence of the correlated data, we define the following linear models for the latent variables,

\begin{matrix} Y_{i}^{1 *} = γ_{0}^{1} + γ_{1}^{1} x_{i} + u_{i} + ϵ_{i}^{1}, \\ Y_{i}^{2 *} = γ_{0}^{2} + γ_{1}^{2} x_{i} + u_{i} + ϵ_{i}^{2}, \end{matrix}

(2.3)

where $ϵ_{i}^{l}$ are independent error terms which are assumed to have a distribution with cumulative density function $F^{l} (ϵ_{i}^{l})$ .

Under this model, the probability of observing a positive outcome is given by,

p_{i}^{l} = P (Y_{i}^{l *} > 0) = P (ϵ_{i}^{l} > - η_{i}^{l}) = 1 - F^{l} (- η_{i}^{l})

where $η_{i}^{l} = γ_{0}^{l} + γ_{1}^{l} x_{i} + u_{i}$ is the linear predictor in Equation (2.2) and Equation (2.3), $l = 1,2$ . If the distribution of the errors is symmetric around zero, we have that,

p_{i}^{l} = 1 - F^{l} (- η_{i}^{l}) = F^{l} (η_{i}^{l}) .

Therefore, the distribution of the error terms in the latent variable model is completely related with the link function of the original model in Equation (2.2).

The multidimensional BBR model in Equation (2.2) is based on a logit link function which defines error terms following the standard logistic distribution. This conclusion allows a complete correlation study between latent variables associated which each dimension for the $i$ th individual as,

ℂ orr [Y_{i}^{1 *}, Y_{i}^{2 *}] = \frac{ℂ ov [Y_{i}^{1 *}, Y_{i}^{2 *}]}{\sqrt{V ar [Y_{i}^{1 *}]} \sqrt{V ar [Y_{i}^{2 *}]}} = \frac{σ_{u}^{2}}{σ_{u}^{2} + π^{2} / 3}, i = 1, \dots, n,

(2.4)

where $π^{2} / 3$ is the variance of the standard logistic distribution and $σ_{u}^{2}$ is the variance of the random effects. Therefore, the variance of the random-effects determines a specific correlation structure between the observations of the same individual in different dimensions or outcomes.

Comparison of the results provided by the multidimensional and the unidimensional approaches are not straightforward. It is known that differences in the behaviour of regression coefficients in the marginal and conditional models are based on the fact that it is not possible to compare them like with like. In fact, it was shown by Lee and Nelder (2004) that these differences are mainly caused by the choice of unidentifiable constraints on the random-effects.

For the marginal expectation of each individual outcome, based on the relationship defined in the multidimensional model given by Equation (2.2), we have that for $l = 1,2$ ,

E [Y_{i}^{l}] = m_{l} E [\frac{\exp (γ_{0}^{l} + γ_{1}^{l} x_{i} + u_{i})}{1 + \exp (γ_{0}^{l} + γ_{1}^{l} x_{i} + u_{i})}] \approx m_{l} Φ (\frac{c}{\sqrt{1 + c^{2} σ_{u}^{2}}} (γ_{0}^{l} + γ_{1}^{l} x_{i})),

(2.5)

where $Φ$ is the cumulative function of the standardized normal distribution and $c = (16 \sqrt{3}) / (15 π)$ . Derivation of the formula can be found in Section S3 of the Supplementary Material.

Alternatively, if we calculate the marginal expectation of the outcome based on the separate analysis where the $l^{t h}$ dimension is modelled as,

\begin{matrix} Y_{i}^{l} \sim BB (m_{l}, p_{i}^{l}, ϕ_{l}), iid i = 1, \dots, n \\ logit (p_{i}^{l}) = β_{0}^{l} + β_{1}^{l} x_{i}, \end{matrix}

(2.6)

we have that,

E [Y_{i}^{l}] = m_{l} E [\frac{\exp (β_{0}^{l} + β_{1}^{l} x_{i})}{1 + \exp (β_{0}^{l} + β_{1}^{l} x_{i})}] \approx m_{l} Φ (c (β_{0}^{l} + β_{1}^{l} x_{i})) .

(2.7)

Therefore, equalizing Equation (2.5) and Equation (2.7) we derive an expression to compare results provided by the two models,

β_{r}^{l} = \frac{γ_{r}^{l}}{\sqrt{1 + c^{2} σ_{u}^{2}}}, r = 0,1, and l = 1,2.

(2.8)

As aforementioned, model formulation for $l$ any number of dimensions ( $l > 2$ ) and for $t$ covariates ( $t > 1$ ) is straightforward from the previously described particular methodology for $l = 2$ and $t = 1$ .

Finally, it is worth mentioning that the computational cost of the estimation in mixed-effects models, specially with non-exponential family distributions, is very high. In fact, one of the main characteristics of the shared random-effects approach is that the number of the dimensions involving the analysis does not alter the dimensionality of the joint density as increasing the number of dimensions does not alter the number of random-effects in the model,

f (y) = \int_{- \infty}^{\infty} f (y^{1}, \dots, y^{L} | u) d u = \int_{- \infty}^{\infty} \prod_{l = 1}^{L} f (y^{l} | u) f (u) d u,

where $L$ is the number of dimensions.

3 Comparing unidimensional and multidimensional beta-binomial regression approaches

In this section, we perform a simulation study to illustrate and compare the main differences when fitting to correlated beta-binomial outcomes univariate beta-binomial regressions and our multivariate beta-binomial proposal with the shared random-effects aproach. The goal of the simulation study is two-fold: a) analyse the performance of the multidimensional BBR model when estimating parameters in terms of bias, dispersion and coverage probabilities; b) illustrate the misleading conclusions in terms of interpretation and covariate significance when univariate models are applied to correlated outcomes.

3.1 Simulation scenarios

We have generated 100 random realizations of $n = 100$ observations of a vector of outcomes of dimension $L = 3$ , where each component $Y^{l}$ conditioned on some simulated random-effects follows a beta-binomial distribution with a fixed maximum score number, $m_{l}$ , a probability parameter $p^{l}$ and a dispersion parameter $ϕ_{l}$ , $l = 1,2,3$ , namely,

\begin{matrix} Y_{i}^{l} | u_{i} \sim BB (m_{l}, p_{i}^{l}, ϕ_{l}) indp . i = 1, \dots, n \\ logit (p_{i}^{l}) = γ_{0}^{l} + γ_{1}^{l} x_{i} + u_{i}, a n d u_{i} \sim N (0, σ_{u}^{2}) . \end{matrix}

For the sake of clarity, we only consider a single covariate in the linear predictor ( $t = 1$ ), which follows an Uniform distribution in the $(0,1)$ interval. Thus, only two parameters, $γ_{0}^{l}$ and $γ_{1}^{l}$ , were included for each dimension, $l = 1,2,3$ , in the fixed effects. Realizations of the random-effects were simulated assuming a normal distribution of mean 0 and standard deviation $σ_{u}$ .

The simulation study has been divided into several scenarios that depend on the variability of the random-effects, $σ_{u}$ , and the parameters of the beta-binomial distribution, $m_{l}$ , $ϕ_{l}$ and $p^{l}$ , for each dimension, where $p^{l}$ is fixed by the vector of coefficients of the linear predictor, $γ^{l} = (γ_{0}^{l}, γ_{1}^{l})^{'}$ , $l = 1,2,3$ . The values of the parameters in the simulation scenarios were selected to reflect different forms for the outcome, a variety of relationships between the predictor and the outcome, and different values for the variance of the random-effect, which dominates the correlation structure in the multidimensional outcome. The combination of values that we selected for the definition of the scenarios are shown in Table 1.

Table 1

Simulated scenarios based on the values of the parameters: $γ^{'} = {(γ_{0}^{'}, γ_{1}^{'})}^{'}, ϕ_{I}$ and m_l, with σ_u = (0,0.5,1,1.5)

	$l = 1$	$l = 2$	$l = 3$
$γ_{0}^{l}$	1	−2	0.5
$γ_{1}^{l}$	−1.5	2.5	0.5
$ϕ_{l}$	0.5	1	1.5
$m_{l}$	20	10	4

Simulated scenarios based on the values of the parameters: $γ^{l} = (γ_{0}^{l}, γ_{1}^{l})^{'}$ , $ϕ_{l}$ and $m_{l}$ , with $σ_{u} = (0,0.5,1,1.5)$

In particular, the case of $σ_{u} = 0$ stands for null random-effects in the model, therefore, the absence of within-subjects correlation between the outcomes. Additionally, we considered values of $σ_{u}$ equal to $0.5,1$ and $1.5$ to analyse scenarios with correlated structures between different dimensions. It is worth mentioning that when beta-binomial mixed-effects simulation studies have been performed in the literature, the variability of the random effects, or in other words the correlation structure, has been defined in a more cautious way. For instance, Wu et al. (2017) defined a standard deviation of the random effects equal to $0.5$ when simulated the performance of the developed Generalized Estimating Equation (GEE) procedure to perform estimations in beta-binomial mixed-effects models. Therefore, we know that some of the proposed scenarios are extreme cases, particularly because correlation structures between dimensions are lower in real practise. Nevertheless, we decided to include larger values of the standard deviation of the random effects in the simulation study to analyse the performance of our proposal in those extreme situations.

The function that implements the multidimensional BBR model is called BBmm and it is available in the PROreg R-package version 1.2.

3.2 Simulation results

Several conclusions can be obtained from this simulation study. A detailed summary of the results is reported in Table S2 of the Supplementary Material. On the one hand, we can appreciate that in the first two scenarios, no-correlation and a cautious value of $σ_{u}$ , the mean estimates of the fixed parameters show a small bias and small mean squared error (MSE). Regarding the dispersion of the estimates, similar values of the empirical standard deviation (ESD), the average estimated standard deviation (ASD) provides evidence about the assumed empirical distribution of the estimators. In terms of the percentage the real value of the fixed effects is included in the estimated 95% confidence interval (PCI), two of the dimensions offer adequate values in both scenarios. However, although the remaining one gets similar results in the rest of the previously analysed parameters, it shows PCI values lower than the empirical value of the 95%, specially in the intercept term.

On the other hand, as the defined standard deviation of the random effects increases and we consider more extreme scenarios ( $σ_{u} = 1$ and $σ_{u} = 1.5$ ), we can appreciate that although estimated mean values, and bias, are still quite near from the real values, the difference between the ESD and ASD begin to increase, leading to PCIs that distance from the expected 95%. However, it is worth mentioning that, although in $σ_{u} = 1$ scenario the PCI are lower than in previous, three out of six are located around the 95% threshold assumed in the model. Additionally, even in the largest variability scenario, the mean estimates are not too far from the true value and PCIs are located around 88%.

3.3 Interpreting the estimated effects

In this section we show the different interpretation between the estimates of the fixed effects when unidimensional and multidimensional approaches are applied to several correlated outcomes. The simulation study shows that misleading interpretations may occur if independent unidimensional regression models are applied to outcomes that are correlated.

Figure 1 shows the box plots of the estimated coefficients for the slope in both unidimensional and multidimensional approaches for each dimension in each defined scenario. One could interpret the deviation showed in Figure 1 as bias, however, as it was mentioned before, the fact that the approaches are measuring the effect of the covariates in a different way, marginally and conditionally, makes the direct comparison inappropriate. In fact, the results are along with the relationship in Equation (2.8) developed to compare regression coefficients in marginal and conditional approaches. It can be appreciated that, results are very similar in the null random-effects scenario ( $σ_{u} = 0$ ), no correlation. On the contrary, following Equation (2.8), it can be appreciated in Figure 1 that the magnitude of the estimates from the unidimensional analysis are smaller as the variance of the random effects, that is, the correlation between the dimensions, increases.

Figure 1:

Note: The red line represents the true value, $R = 100$ realizations.

In PROs analysis, as the same individuals respond to the items that construct the outcome variables, individual-specific characteristics may mask the relationship between the outcomes and covariates, leading to erroneous marginal effect interpretation. The good point of the conditional approach is that it fixes the effects of the covariates by individual-specific characteristics, which concludes with more realistic effects of the covariates in the outcome. If we used a marginal approach to make the interpretations, the effect will always show a more moderate relationship (in absolute value) than the conditional as has been shown in Figure 1. Therefore, this approach could lead to erroneous interpretations. Moreover, following Equation (2.8), we could transform the conditional estimates to interpret them in a marginal approach if there is some interest in that. However, the contrary cannot be done, as with the marginal approach, we do not estimate any variance of the random effects that allows transforming the estimates to a conditional approach.

Additionally, the limitation of the unidimensional approach to accommodate the correlation existing in the data tends to overestimate the dispersion parameter of the assumed beta-binomial distribution. Figure 2 shows the estimations of the dispersion parameter of the beta-binomial distribution, $ϕ^{l}$ , in the $L = 3$ dimensions for each simulation scenario. It can be appreciated that, while in the non-correlation scenario both approaches perform similarly, as the correlation increases the estimation of $ϕ$ in the unidimensional approach tends to depart from the real value.

Figure 2:

Note: Red line stands for the true value, $R = 100$ realizations.

The simulation study shows that the unidimensional approach obtains more cautious (smaller in absolute value) estimates of the regression coefficient as the variance of the random effects increases. Moreover, Figure 2 shows that the unidimensional approach tends to overestimate the dispersion parameter of the beta-binomial distribution, increasing the variance of the estimations of the regression parameters. Therefore, when the variance of the random effects, that is, the correlation between the responses provided by the same individual, is big enough, the unidimensional approach can lead to covariate effects that are not statistically significant, as the estimations of the regression coefficients are cautious and have high variance.

4 Analysis of COPD data with beta-binomial approaches

In this section, we describe the application of the multidimensional beta-binomial regression model defined in Section 2 to the patients of the COPD study. The objective is to test for the impact of several patients’ and disease characteristics on the activity, impacts and symptoms dimensions provided by the SGRQ.

Previously, Esteban et al. (2016) identified four subtypes of patients with COPD in the presented study, three of them had marked respiratory profiles with a continuum in the severity of clinical variables, while the fourth subtype had a more systemic profile, with intermediate severity as regards to respiratory disease. Although this subtype classification accounted for 73% of the variability in the original clinical variables, researchers believed that the remaining variability could still be relevant to the HRQoL of the patients. Therefore, we decided to include in the analysis not only the subtype but also other clinical variables collected in the study and candidates to significantly influence HRQoL beyond the subtype.

The beta-binomial fitting requires a previous recoding of the original scores to a binomial form as described in Arostegui et al. (2013). The recoding of the SGRQ scores to a binomial form has been performed based on a minimum clinically important difference. Jones (2002) showed that 4 points of change in the SGRQ could be considered as a clinically significant change and, hence, we have divided the 0–100 scale into 4 subintervals of equal length. As a result of the recoding process, the recoded SGRQ scores are presented as ordinal scores with a range from 0 to 24. Figure 3 shows the fitted beta-binomial distribution of the SGRQ scores.

Figure 3:

Histograms of the recoded values of the scores for the three dimensions of the SGRQ along with the estimated density function for the fitted beta-binomial distribution: $B B (24,0.4845,0.3066)$ for the activity dimension, $B B (24,0.3134,0.2271)$ for the impacts dimension and $B B (24,0.4435,0.2228)$ for the symptoms dimension.

Health dimensions were correlated as expected. Indeed, Pearson’s correlation coefficients between the three dimensions of the SGRQ were: 0.51 (symptoms - activity); 0.62 (symptoms - impacts) and 0.78 (activity - impacts), being very similar to Spearman’s correlation coefficients. Therefore, independent analysis of each dimension could be inefficient or confusing, and hence, a joint analysis of the three HRQoL dimensions was preferred. Thus, the three dimensions of the SGRQ were included in the model as a multivariate outcome. Nevertheless, an independent analysis of the three HRQoL dimensions was also performed for comparison purposes.

In terms of software, functions that implement both the multidimensional and unidimensional beta-binomial regression approaches are included in the PROreg R-package version 1.2 (https://CRAN.R-project.org/package=PROreg). Specifically, the unidimensional approach is implemented by BBreg function, while the function that implements the multidimensional approach is named BBmm. More information about the functions and implementing examples are provided in the package’s documentation.

4.1 Results

Table 2 shows the results of the multidimensional beta-binomial regression model fitted to COPD data along with the results of the unidimensional analysis performed for the activity dimension.

Table 2

Results of the multidimensional and the unidimensional beta-binomial regression models for the activity dimension of the SGRQ.

SGRQ		Multidimensional			Unidimensional
Dimension	Covariate	$\hat{β}$	$S D (\hat{β})$	$p$	$\hat{β}$	$S D (\hat{β})$	$p$
Activity	Mild COPD
	Moderate COPD	0.488	0.109	$< 0.001$	0.472	0.154	0.002
	Severe COPD	1.007	0.135	$< 0.001$	0.923	0.190	$< 0.001$
	Systemic COPD	0.280	0.121	0.020	0.283	0.170	0.096
	Age $\leq 65$
	Age $(65,75]$	−0.376	0.075	$< 0.001$	−0.360	0.105	$< 0.001$
	Age $> 75$	−0.605	0.093	$< 0.001$	−0.558	0.132	$< 0.001$
	Treatment 1
	Treatment 2	0.286	0.107	0.007	0.273	0.147	0.064
	Treatment 3	0.322	0.095	$< 0.001$	0.289	0.131	0.028
	6mWT (25m)	−0.052	0.011	$< 0.001$	−0.043	0.015	0.005
	FEV1 (100ml)	−0.042	0.008	$< 0.001$	−0.041	0.011	$< 0.001$
	PA: $< 2$ h/w
	PA: 2–4 h/w	−0.643	0.130	$< 0.001$	−0.575	0.181	0.002
	PA: $> 4$ h/w	−0.590	0.132	$< 0.001$	−0.541	0.184	0.003
	PA: Intense	−0.779	0.151	$< 0.001$	−0.699	0.211	$< 0.001$
	CCI $< 2$
	CCI $\geq 2$	0.164	0.073	0.025	0.130	0.103	0.210
	PrevHospi 0
	PrevHospi 1–2	0.058	0.080	0.472	0.069	0.113	0.544
	PrevHospi $> 2$	0.046	0.181	0.798	0.022	0.243	0.928

Note: SGRQ: St George Respiratory Questionnaire; COPD: Chronic obstructive pulmonary disease; FEV1: forced expiratory volume in the first second; 6mWT: 6-minute walking test; CCI: Charlson comorbidity index; PA: Physical activity; PrevHospi: Number of hospitalizations in the previous year.

Similarly, Tables S3–S4 in the Supplementary Material show the same results for the symptoms and impacts dimensions. When we compare results from the unidimensional analysis with the results from the multidimensional analysis, we do not observe important differences in the estimation of the beta coefficients, as the estimated random effects’ standard deviation was not very large. Nevertheless, we observe remarkable differences in their variances. The variance of the estimated beta coefficients is consistently larger in the unidimensional approach than in the multidimensional approach, as shown in the simulation study due to the overestimation of the beta-binomial dispersion parameter. In consequence, the significance of the coefficients is more likely in the latter than in the former approach. For instance, from the multidimensional approach, we conclude that subjects with moderate, severe or systemic COPD had significantly poorer HRQoL on activity and impacts than subjects with mild COPD. However, neither the difference between systemic and mild in both dimensions nor the difference between moderate and mild in activity was significant in the unidimensional approach. Many other differences between both approaches could be observed in the three tables, all of them in the same direction, showing that the multidimensional approach has more capacity to detect significant covariates than the unidimensional approach.

Table 3 shows the results of the dispersion parameters obtained in the multidimensional and unidimensional beta-binomial regression models. Values of the dispersion parameter of the beta-binomial distribution for the three dimensions of the SGRQ questionnaire are shown on a logarithmic scale and they were consistently smaller and had more variability for the multidimensional approach than for the unidimensional approach.

Table 3

Estimation of the dispersion parameters in the beta-binomial regression models.

		Multidimensional		Unidimensional
Parameter	SGRQ Dimension	Estimate	SD	Estimate	SD
$\ln (\hat{ϕ})$
	Activity	−2.830	0.116	−1.491	0.073
	Impacts	−3.937	0.219	−1.746	0.075
	Symptoms	−2.441	0.090	−1.672	0.072
$\hat{σ}$		0.827	0.026	–	–

Note: SGRQ: St George Respiratory Questionnaire; SD: standard deviation.

As previously mentioned, the $ϕ$ coefficient captures the correlation between items regarding the same HRQoL dimension. It was also shown in the simulation study (Figure 2) that the unidimensional approach tends to overestimate the parameter $ϕ$ as it does not take into account the within-subject variability in the dimensions. Based on the result provided by the unidimensional approach, we would overestimate the correlation between items in an HRQoL dimension, leading to misleading conclusions about the way the questionnaire is constructed. The estimated variance of the random effect that accounts for the variability within individuals was 0.684. Consequently, we can estimate the correlation between the three HRQoL dimensions within individuals accounted by the random effects from Equation (2.4), which takes the value 0.173.

5 Discussion

A proposal for the joint analysis of various dimensions from the same questionnaire as a multivariate outcome has been presented, specifically designed for the beta-binomial regression framework, which has been proposed in the literature for the analysis of overdispersed binomial data, such as PROs. The method was developed in the context of mixed-effects models based on a shared random effect for all the dimensions to be jointly analyzed and estimation was performed using a penalized joint log-likelihood approach. The model can be implemented in the PROreg R-package version 1.2.

Wu et al. (2017) proposed a Generalized Estimating Equation (GEE) procedure to overcome the estimation of the parameters in a beta-binomial mixed effects model, where the application was limited to longitudinal beta-binomial data. However, we have based the inferential process of our proposal on the conditional joint log-likelihood approach developed by Najera-Zuloaga et al. (2019). We decided to extend the joint log-likelihood approach mainly for two reasons: (i) the conditional approach offers the possibility to estimate subject-specific effects that, in the PROs framework, are of main importance to evaluating the psychological state of each patient, and; (ii) GEE approach is more restrictive regarding missing data (Gibbons et al., 2010), which could be quite common in this type of clinical studies. Finally, it is also worth mentioning that the logit link function is the most widely used link function when analysing PROs, however, the logistic regression requires more computing effort than the probit regression with GEE method (Wu et al., 2017).

Generally speaking, when several health dimensions were measured within the same subjects (i.e., there exists a correlation between dimensions), the simulation study showed that results from the individual analysis and multidimensional approach could lead to misleading interpretations of the effect of the covariates.

Indeed, the differences between the unidimensional and the multidimensional results increase considerably as the correlation between responses increases. If we pay attention to the coefficients of the fixed effects, specifically the slope, which quantifies the effect of the covariate $x$ in the outcome, we observe that the unidimensional approach provides lower estimates in magnitude, that is, tends to provide more conservative estimates of the effects of the covariates, which can also be concluded with Equation (2.8). In practice, conservative estimates together with an overestimation of the dispersion parameter would imply that the effect of some significant variables will not be detected in the independent analysis of the outcomes, which in fact may also be observed in the analysis of the real data.

The main differences between the estimated regression coefficients in the unidimensional and the proposed multidimensional approach are derived by the fact of comparing directly estimates obtained in conditional and marginal models. Once the estimates in the conditional approach have been transformed to marginal using Equation (2.8), both models end with very similar estimates. We believe that the interpretation of the results should be done in a conditional way, where the effect of the covariates is fixed by the individual specific characteristics. However, the marginal model estimates cannot be transformed to conditional as the model does not provide any variance of the random effects, which is essential to apply Equation (2.8). Beyond the discussion of the interpretation of the regression coefficients in conditional or marginal approaches, it is worth mentioning that as it was shown in the COPD data analysis since the multivariate approach accommodates the individual correlation, the variance of the regression estimates reduces considerably, increasing the capacity to detect statistically significant covariates. Results of the estimation of the dispersion parameter shown in Figure 2 are along with the previous statement where the overestimation of the dispersion parameter increases the standard error of the estimates of the regression coefficients in the individual regressions.

According to the rules for discriminant validity in the design of PROs, scale-to-scale correlation in a questionnaire must be lower than the reliability index (Fayers and Machin, 2007). Reliability is generally measured by Cronbach’s alpha coefficient, which is considered acceptable for values higher than 0.7 (Nunnally and Bernstein, 1994). Thus, the correlation between dimensions of the same questionnaire, although expected to be low, could be as high as 0.7. For instance, the correlation between the three dimensions of the SGRQ in the validation study of the Spanish version, measured by the Pearson correlation coefficient, was 0.57 (symptoms - activity), 0.65 (symptoms - impacts) and 0.79 (activity - impacts). These values are very similar to the values we found in our study. The maximum value that we have considered for $σ_{u}$ in the simulated scenarios is 1.5, which leads to a correlation of 0.406. Nevertheless, we should mention that application to real data shows that the random effect does not account for all the existing correlations between the dimensions. The estimated correlation between the dimensions was 0.17, whereas the empirical correlation ranged from 0.51 to 0.78. However, part of the correlation between dimensions in real cases could be due to the effect that some common covariates have on the outcomes. When this effect is accounted for with the inclusion of the fixed effects, the residual correlation is what is estimated through the random effect. This last result strengthens the extreme definition of the simulation study when $σ_{u}$ was fixed and equalized to $1.5$ .

Although we have applied the model assuming that all the dimensions come from the same questionnaire, several PROs could be measured for the same subjects, showing information on different health aspects. For instance, the subjects of the COPD study were asked for generic HRQoL, disease-specific HRQoL and anxiety and depression symptoms, each of the outcomes measured by a different questionnaire with 8, 3 and 2 dimensions, respectively. The three PROs, and hence the 13 health dimensions, could be jointly analysed using the multidimensional BBR model. In addition, the number of dimensions involved in the analysis does not affect the computational cost, as increasing the number of dimensions does not alter the number of random effects in the model.

The main limitation of the shared random-effects approach presented here is the rigidity of a common unique random effect for capturing the correlations between all the dimensions, which could be quite unrealistic in some real practices. A very conservative recommendation could be to limit the proposed method to questionnaires with only two dimensions, with a clear common correlation structure (Martin et al., 2019). However, our study showed that including several dimensions in a multidimensional approach is more appropriate than performing independent analyses for each dimension, especially when the correlation between dimensions is high.

Theoretically, the extension of the shared random effects to longitudinal data is straightforward. However, when analysing repeated measurements with this approach, the correlation between different dimensions is determined by the dimension’s intrinsic correlation given by the repeated measurements. Therefore, the correlation structure within each dimension dictates the association between the different dimensions, which could be too rigid and unrealistic in real situations. The restriction could be relaxed by allowing different but correlated random effects for the various dimensions. Nevertheless, it would increase the dimensionality of the joint marginal density, which combined with the complexity of the beta-binomial distribution itself, could lead to inappropriate approximations of Laplace’s approximation method (Shun and McCullagh, 1995).

Finally, when the proposed methodology is applied to patients with COPD to detect and quantify the effect of some factors in disease-specific HRQoL, results are highly valuable for clinical researchers. For instance, as regards comorbidities, a cut-off of three comorbidities has been mentioned in the literature to establish significant differences in HRQoL (Jones et al., 2011). In our study, this result was not obtained through independent analysis, whereas it was detected through multidimensional analysis. In fact, Charlson’s comorbidity index greater or equal to 3 was a significant variable, showing a minimal clinically important difference in activity and symptoms for patients with the same other characteristics. In addition, the model estimated that the subtype previously defined as systemic COPD, with a high presence of comorbidities (Esteban et al. reported a mean Charlson’s comorbidity index of 4.8), had a significantly lower HRQoL than the subtype with a mild COPD, both in the dimension of activity and impact (Esteban et al., 2016).

In summary, we propose a multidimensional approach to analyze several PROs in a joint multivariate beta-binomial regression framework that offers more reliable results than independent beta-binomial regression analyses of the outcomes. Therefore, we recommend using the multidimensional BBR approach to analyze the PROs when more than one dimension is included in the questionnaire.

Footnotes

Acknowledgements

We gratefully acknowledge the clinical researchers for giving us the opportunity to participate in the data analysis of the COPD study that motivated this work and also all the patients who participated in the study.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

Financial support for this study was provided, in part, by the Department of Education, Language Policy and Culture from the Basque Government (IT1456-22) and BERC 2022-2025 program), the Spanish Ministry of Economy and Competitiveness MINECO and FEDER (PID2020-115882RB-I00) funded by (AEI/FEDER, UE) and acronym ‘S3M1P4R’, and the Ministry of Science and Innovation: BCAM Severo Ochoa accreditation CEX2021-001142-S/MICIN/AEI/10.13039/501100011033.

References

Arostegui

, Nu´n˜ez-Anto´n

and Quintana

(2007) Analysis of the Short Form-36 (SF-36): The beta-binomial distribution approach. Statistics in Medicine , 26, 1318–1342.

Arostegui

, Nu´n˜ez-Anto´n

and Quintana

(2012) Statistical approaches to analyse patient-reported outcomes as response variables: An application to health related quality of life. Statistical Methods in Medical Research , 21, 189–214.

Arostegui

, Nu´n˜ez-Anto´n

and Quintana

(2013) On the recoding of continuous and bounded indexes to a binomial form: An ap- plications to quality-of-life scores. Journal of Applied Statistics , 40, 563–582.

Chang

C-H

and Reeve

(2005) Item response theory and its applications to patient- reported outcomes measurement. Evalua- tion & the Health Professions , 28, 264–282.

de Vet

HCW

, Terwee

, Mokkink

and Knol

(2011) Measurement in Medicine . Cam- bridge, UK: Cambridge University Press.

Esteban

, Arostegui

, Aburto

, Moraza

, Quintana

, Garc´ıa-Loizaga

, Basualdo

, Aramburu

, Aizpiri

, Uranga

and Capelastegui

(2016) Chronic Obstructive Pulmonary Disease Subtypes. Transitions over time. PLOS ONE , 11, e0161710. doi: 10.1371/journal.pone.0161710.

Etxeberria

, Urdaneta

and Galdona

(2019) Factors associated with health related qual- ity of life (HRQoL): Differential patterns de- pending on age. Quality of Life Research , 28, 2221–2231.

Fayers

and Machin

(2007) Quality of life. The Assessment, Analysis and Interpretation of Patient-Reported Outcomes (2nd ed). Lon- don: John Wiley and Sons.

Forcina

and Franconi

(1988) Regression anal- ysis with the beta-binomial distribution. Riv- ista di Statistica Applicata , 21.

10.

Gardiner

, Luo

and L-A

Roman

(2019) Fixed effects, random effects and GEE: What are the differences? Statistics in Medicine , 2, 221–239. URL http://doi.wiley.com/10.1002/sim.3478.

11.

GBD 2016 Disease and Injury Incidence and Prevalence Collaborators. (2018). Global, re- gional, and national incidence, prevalence, and years with disability for 354 diseases and injuries for 195 countries and territo- ries, 1990–2017: A systematic analysis for the Global Burden of Disease Study. The Lancet , 392, 1789–1858.

12.

Gibbons

, Hedeker

and DuToit

(2010) Analysing cognitive test data: Distributions and non-parametric random effects. Annual Review of Clinical Psychology , 6, 79–107.

13.

Johnson

, Kemp

and Kotz

(2005). Univariate Discrete Distributions (3rd ed.). Wiley Series in Probability and Statistics.

14.

Jones

, Quirk

and Baveystock

(1991) The St George’s Respiratory Questionnaire. Respiratory Medicine , 85, 25–31.

15.

Jones

(2002) Interpreting thresholds for a clin- ically significant change in health status in asthma and COPD. European Respiratory Journal , 19, 398–404.

16.

Jones

, Brusselle

, Dal Negro

, Ferrer

, Kardos

, Levy

, Perez

, Soler-Catalun˜a

, van der Molen

, Adamek

and Banik

(2011) Health-related quality of life in pa- tients by COPD severity within primary care in Europe. Respiratory Medicine , 105, 57–66.

17.

Lee

and Nelder

(2001) Hierarchical general- ized linear models: A synthesis of generalized linear models, random-effects and structured dispersions. Biometrika , 88, 987–1006.

18.

Lee

and Nelder

(1996) Hierarchical general- ized linear models. Journal of the Royal Sta- tistical Society , 58, 619–678.

19.

Lee

and Nelder

(2004) Conditional and marginal models: Another view. Statistical Science , 19, 219–238.

20.

Martin

, Arostegui

, Loron˜o

, Padierna

, Najera-Zuloaga

and Quintana

(2019) Anxiety and depressive symptoms are re- lated to core symptoms, general health out- come and medical comorbidities in eating disorders. European Eating Disorders Review , 27, 600–613. doi: 10.1002/erv.2677.

21.

McCullagh

and Nelder

(1989) Generalized Linear Models (2nd edition). London: Chap- man & Hall.

22.

McCulloch

(2008) Joint modelling of mixed outcome types using latent variables. Statis- tical Methods in Medical Research , 17, 53– 73.

23.

Molenberghs

, Verbeke

, Iddi

and Deme´trio

(2012) A combined beta and normal random-effects model for repeated, overdis- persed binary and binomial data. Journal of Multivariate Analysis , 111, 94–109.

24.

Najera-Zuloaga

, D-J

Lee

and Arostegui

(2018) A comparison of beta-binomial regression model approaches to analyze health-related quality of life data. Statistical Methods in Medical Research , 27, 2989–3009. doi: 10.1177/0962280217690413.

25.

Najera-Zuloaga

, D-J

Lee

and Arostegui

(2019) Beta-binomial mixed-effects model for analyzing longitudinal binomial data with overdispersion. Biometrical Journal , 61, 600–615. doi: 10.1002/bimj.201700251.

26.

Nunnally

and Bernstein

(1994) Psychometric Theory (3rd ed.). New York: McGraw-Hill.

27.

Pawitan

(2001) In All Likelihood: Statistical Modelling and Inference Using Likelihood . Oxford University Press.

28.

Reinsel

(1984) Estimation and prediction in a multivariate random-effects generalized lin- ear model. Journal of the American Statistical Association , 79, 406–414.

29.

Rigby

and Stasinopoulos

(2005) General- ized Additive Models for Location, Scale and Shape. Applied Statistics , 54, 507–554.

30.

Shun

and McCullagh

(1995) Laplace approx- imation of high dimensional integrals. Series B: Journal of the Royal Statistical Society, 57, 749–760.

31.

Testa

and Simonson

(1996) Assessment of quality of life outcomes. New England Jour- nal of Medicine , 334, 835–840.

32.

Weldring

and Smith

(2013) Patient-Reported Outcomes (PROs) and Patient-Reported Outcome Measures (PROMs). Health Services Insights , 6, 61–68.

33.

Williams

(1982) Extra-binomial variation in lo- gistic linear models. Series C: Journal of the Royal Statistical Society, 31, 144–148.

34.

, Zhang

and Long

(2017) Longi- tudinal beta-binomial modeling using GEE for overdispersed binomial data. Statistics in Medicine , 36, 1029–1040.

Multidimensional beta-binomial regression model: A joint analysis of patient-reported outcomes

Abstract

Keywords

1 Introduction

2 A beta-binomial model for multivariate outcomes

3.1 Simulation scenarios

Table 1

Simulated scenarios based on the values of the parameters: γ ′ = γ 0 ′ , γ 1 ′ ′ , ϕ I and ml, with σu = (0,0.5,1,1.5)

3.3 Interpreting the estimated effects

Figure 3:

Table 2

Results of the multidimensional and the unidimensional beta-binomial regression models for the activity dimension of the SGRQ.

Estimation of the dispersion parameters in the beta-binomial regression models.

Footnotes

Acknowledgements

Declaration of Conflicting Interests

Funding

References

Simulated scenarios based on the values of the parameters: $γ^{'} = {(γ_{0}^{'}, γ_{1}^{'})}^{'}, ϕ_{I}$ and m_l, with σ_u = (0,0.5,1,1.5)