BAYESIAN ESTIMATION FOR THE STABLE DISTRIBUTIONS IN THE PRESENCE OF COVARIATES WITH APPLICATIONS IN CLINICAL ISSUES

ABSTRACT

In this paper we explore a Bayesian approach for stable distributions in presence of covariates. This class of distribution has great flexibility for fitting asymmetric and heavy-tailed empirical data. These models are commonly used for data sets in finance and insurance. In this paper we show that these distributions can also be used to fit clinical data. Since there is not an analytical form for the density probability function which implies in serious difficulties to obtain the maximum likelihood estimators for the parameters, we use Bayesian methods with data augmentation techniques to get the inferences of interest. In this study we also discuss the choice of different prior distributions for the parameters considering regression models for the location and scale parameters of the stable distribution. We use MCMC (Markov Chain Monte Carlo) algorithms to generate samples from the posterior distributions in order to evaluate the point and interval estimators. A great simplification is obtained using the OpenBugs software. Two real data examples illustrate the applicability of the proposed modeling approach.

Keywords:
table distributions; Bayesian approach; regression models; prior distributions; MCMC methods

1 INTRODUCTION

In many areas of applications, the usual assumption of normality of the data needed in many statistical models is not satisfied in practice, even using data transformations to improve the symmetry of the data. The assumption of normality is essential for the use of traditional techniques such as ANOVA (analysis of variance) models, linear regression models, hypothesis tests for comparison of means, medians or variances, especially with small sample sizes. In this situation, usually the data analyst uses existing non-parametrical methods, but this approach usually gives poor statistical inference results, and in general it is not possible to get predictions, among many other problems. A recent alternative that has been explored in the literature given the advancement of computational methods and better computers is the use of super models that generalize the normal distribution allowing the adjustment to asymmetric data, multimodal data, among others. In this case we can cite the family of stable distributions that have already been used in areas such as economics and finance. The family of stable distributions arises from a generalized version of the central limit theorem (Gnedenko & Kolmogorov, 1954) where the Gaussian distribution is a particular case of this family of distributions. The stable distribution has great flexibility to fit data with very different shapes, allowing to fit asymmetric and heavy-tailed data, which makes it suitable for modeling in different research areas such as complex systems in physics, engineering, biology, sociology and other fields.Several analytical properties of this class of stable distributions and several practical applications can be seen in Zolotarev (1986).

A major difficulty associated with stable distributions is that, in general, there is no simple and closed formula for the probability density functions (pdf), making it difficult to apply well-known estimation methods, such as parameter estimation using maximum likelihood methods. In theory, the pdf of a stable distribution could be obtained numerically from an integration based on the inversion formula or using a characteristic function-based estimation method, but these techniques are often complex and difficult to be used in practice.

Borak et al. (2005) and Kateregga et al. (2017) reviewed some classical approaches proposed to estimate the parameters of stable distributions and investigated their performance from simulation studies. They observed that the inferences obtained by these approaches can present several computational difficulties and often inaccurate estimates of model parameters.

In the last decades, Bayesian inference methods have been considered as an efficient and robust alternative approach to data analysis, especially using more complex models in general with several parameters and with great flexibility to be fited by the data. The literature presents several studies considering the use of Bayesian estimation methods to get the estimators for the parameters of stable distributions. In this direction, one of the first papers written to explore the use of Bayesian methods for stable distribution was introduced by Buckle (1995) proposing a Bayesian analysis using Markov Chain Monte Carlo (MCMC) methods and data augmentation techniques (introduction of unobserved latent variables) that made it possible to get precise inferences of interest to the stable distribution. Lombardi (2004) introduced a novel approach for Bayesian inference in the setting of stable distributions that resorts to a Fast Fourier Transform (FFT) of the characteristic function in order to approximate the likelihood function; the posterior distributions of the parameters are then produced via a random walk MCMC method. Oral and Erdemir (2012) compared the FFT with their Bayesian approach using Metropolis random walk chain and direct numerical integration.

Achcar et al. (2013a) studied computational aspects for the Bayesian analysis involving stable distributions. Following in the same direction, Achcar et al. (2013b) presented some robustness aspects of linear regression models in the presence of outliers or discordant observations considering the use of stable distributions for the response in place of the usual normality assumptionwhere the main aspects of the methodology, as well as two examples of implementation, were presented: one using data from the industrial area and other from simulated data. The authors elaborated the analysis and discussion including measures of goodness of fit of the models. The influence of outliers in the model fit was also discussed in this publication.

Lemke et al. (2015) considered a fully Bayesian inference for α-stable distributions using a Poisson series representation, and more recently Karling et al. (2021) provided a Bayesian algorithm that makes use only of the power series representation.

The main goal of this paper is to carry out a Bayesian inference with data augmentation techniques and MCMC (Markov Chain Monte Carlo) methods proposed by Buckle (1995) to analyse data under stable distributions in presence of covariates. We provide major applications of the stable distributions in the medical field involving covariates in order to determine important prognostic factors.

We also discuss the choice of appropriate prior distributions for the parameters of the model and some computational issues to obtain accurate Bayesian inferences. The popular existing free Openbugs software (Spiegelhalter et al., 2003) is used to obtain the posterior summaries of interest.

The paper is organized as follows: in Section 2, we review the basic concepts of the stable distribution family and some basic properties of the model. Section 3 presents a Bayesian analysis for the stable distributions not considering the presence of covariates. Section 4 presents a Bayesian analysis of the model assuming uniform prior distributions for the parameters of the model and also considering regression models for the location and scale parameters. Section 5 presents two applications with real medical data sets. Finally, section 6 provides some concluding remarks on the obtained results.

2 DEFINITION AND BASIC PROPERTIES

The stable distribution does not have an analytic closed form for the probability density and distribution functions, but can be expressed by its characteristic function given by

$\varphi \left(t;\alpha ,\beta ,\mu ,\sigma \right)=E\left[e^{itX}\right]=\left[\begin{array}{c}exp\left\{i\mu t-{\left|\sigma t\right|}^{\alpha }(1-i\beta \left(signt\right)\tan \left(\frac{r\alpha }{2}\right)\right\},\alpha \ne 1\\ exp\left\{i\mu t-\alpha \left|t\right|\left(1+i\beta \frac{2}{\pi }\left(signt\right)\ln \left|t\right|\right)\right\},\alpha =1\end{array}\right].$ (1)

where signt = 1 if t > 0; signt = 0 if t = 0 and signt = -1, if t < 0.

The stable distribution, denoted by S _α (β, µ, σ), requires four parameters α, β , µ and σ to complete descrition. The parameter $\alpha \in (0,2]$ defines the index of stability or the shape parameter (heaviness of the tails), and when α = 2 this class reduces to Gaussian distributions. The parameter $\beta \in \left[-1,1\right]$ is the skewness parameter, where for β = 0 we have symmetric distributions. The location and scale parameters are, respectively, $\mu \in \left(-\infty ,\infty \right)$ and $\sigma \in \left(0,\infty \right)$. An important property of the stable distribution is that if a random variable has a stable distribution, that is, $X^{~}{S}_{\alpha }\left(\beta ,\mu ,\sigma \right)$ then $Z=\left(X-\mu \right)/{\sigma }^{~}{S}_{\alpha }\left(\beta ,0,1\right)$.

A great difficulty associated to stable distributions S _α (β, µ, σ) is that in general there is no simple closed form for their probability density functions except for the case of a Gaussian (α = 2), Cauchy (α = 1, β = 0) and Lévy (α = 0.5, β = ±1) distributions. However, it is known that the probability density functions of stable distributions are continuous and unimodal. The support of all stable distributions is given in (-∞, ∞), except for α < 1 and |β| = 1 when the support is (-∞, 0) for β = 1 and (0, ∞) for β = -1. It is important to point out that, if α < 1, the variance is infinite and the mean of the stable distribution does not exist. When α > 1, the mean of the distribution exists and is equal to µ. When the skewness parameter β is positive, the distribution is skewed to the right; when it is negative it is skewed to the left. When β = 0, the distribution is symmetric about µ. As α approaches 2, β loses its effect and the distribution approaches the Gaussian distribution regardless the values of β . If the scale parameter σ = 1 and the location parameter µ = 0, the distribution is called standard stable.

Figure 1 shows the density plots for different parameter values of the stable distribution.

Figure 1
Densities of the stable distribution S _α (β, µ, σ) with parameters µ = 0, σ = 1. Right panel: closed form for the known densities Gauss, Cauchy and Levy. Left panel: stable density functions for α = 1.2 and β = 0, 0.5, 0.8 and 1.

Some properties of the α-stable distribution are:

• The tail of the density function decays like a power function, that is, $P\left(\left|X\right|>x\right)\propto C{x}^{-\alpha }$ when x → ∞ for some constant C.

• $E\left({\left|X\right|}^p\right)<\infty ,0<p<\alpha ;E\left({\left|X\right|}^p\right)=\infty ,p\ge \alpha .$

• $E\left(\left|X\right|\right)=\mu ,if\alpha >1;E\left(\left|X\right|\right)=\infty ,if\alpha \le 1.$

• The stability property is preserved under linear transformation.

Although this class of distributions is a good alternative for data modeling in different areas, we usually have great difficulties to obtain estimates under classical approaches due to the lack of closed form expressions for their probability density functions. One possibility in applications is to compute the probability density function from the inversion formula,

$f\left(x\right)=1/2\pi {\int }_{-\infty }^{\infty }{e}^{-itx}\Phi \left(t\right)dt$ (2)

where Φ(t) is the characteristic function. In applications, we must use numerical methods to solve the integral in (2), usually taking a great computational time.

The literature presents different studies leading to different techniques to compute the density functions of the family of stable distributions S _α (β, µ, σ) (Oral & Erdemir, 2012).

3 BAYESIAN INFERENCES FOR THE PARAMETERS OF THE STABLE DISTRIBUTION

An alternative estimation method to obtain inferences for stable distributions is the use of Bayesian inference (Gelman et al., 2013; O’Hagan & Forster, 2004; Bernardo & Smith, 2007).

In Bayesian inference, a prior distribution is specified for the parameters of the distribution, which reflects existing knowledge (or lack of knowledge) regarding the probability of the possible values. The ability to specify a prior distribution is seen as an advantage of the Bayesian approach because the uncertainty surrounding each parameter is taken into consideration.

It is important to note that, in the Bayesian approach, the inference is based on the posterior density of the parameters, which is proportional to the product of likelihood and prior distribution, that is $p\left(\theta |\mathit{x}\right)\propto f\left(\mathit{x}|\theta \right){\pi }_0\left(\theta \right)$. However, for stable distributions the likelihood function based on the the data is not possible to derive directly from the density (2) since this one is not given in a closed form. Therefore, an alternative to estimate the parameters of S _α (β, µ, σ) is to use latent or artificial variables.

If it is not possible to obtain the probability density function of x in closed form whereas the joint density function of x and y exists, then the posterior density is found by taking integration,

$\pi \left(\theta |x\right)\propto \int f\left(x,y|\theta \right){\pi }_0\left(\theta \right)dy$ (3)

where π ₀ is the prior distribution for the parameter.

Considering a latent or an auxiliary variable, Buckle (1995) proved a theorem that is very useful to simulate samples of the joint posterior distribution for the parameters α, β, γ and δ. This theorem establishes that a stable distribution for a random variable Z defined in (-∞, ∞) is obtained as the marginal of a bivariate distribution for the random variable Z itself and an auxiliary or latent random variable Y. This non-observed random variable Y is defined in the interval (-0.5, a _α,β ), when $Z\in \left(-\infty ,0\right)$, and in (a _α,β , 0.5), when $Z\in \left(0,\infty \right)$.

The quantity a _α,β is given by

$a_{\alpha ,\beta }=-b_{\alpha ,\beta }\alpha \pi$ (4)

where $b_{\alpha ,\beta }=\left(\beta \pi /2\right)min\left\{\alpha ,2-\alpha \right\}$.

The joint probability density function for the random variables Z and Y is given by

$f\left(z,y|\alpha ,\beta \right)=\left\{\alpha /\left|\alpha -1\right|\right\}\left\{1/\left|z\right|\right\}\left\{{\left|z/t_{\alpha ,\beta }\left(y\right)\right|}^{\theta }\right\}exp\left\{-{\left|z/t_{\alpha ,\beta }\left(y\right)\right|}^{\theta }\right\}/\sigma$ (5)

where $\theta =\alpha /\left(\alpha -1\right)$,

$t_{\alpha ,\beta }\left(y\right)=\left\{\sin \left(\pi \alpha y+b_{\alpha ,\beta }\right)/\cos \left(\pi y\right)\right\}{\left\{\cos \left(\pi y\right)/\left[\cos \left(\pi \left(\alpha -1\right)y\right)+b_{\alpha ,\beta }\right]\right\}}^{1/\theta }$ (6)

and $Z=\left(X-\mu \right)/\sigma \text{for}\sigma \ne 0$.

Buckle (1995) showed that the marginal distribution for the random variable Z of the bivariate density (5), is a stable S _α (β, 0, 1) distribution.

Let us assume that $\mathit{x}=\left(x_1,x_2,\dots ,x_n\right)$ is a random sample of size n, where $X~S_{\alpha }\left(\beta ,\mu ,\sigma \right)$ and π ₀(α, β, µ, σ) is a joint prior distribution for the parameters α, β, µ and σ. Then, from the bivariate density (5), the joint posterior distribution for α, β, µ and σ is given by:

$\pi \left(\alpha ,\beta ,\mu ,\sigma |\mathit{x}\right)\propto \in {\left(\frac{\alpha }{\left|\alpha -1\right|\sigma }\right)}^{n}exp\left\{-\sum _{i=1}^n{\left|\frac{z_i}{t_{\alpha ,\beta }\left(y\right)}\right|}^{\theta }\right\}\prod _{i=1}^n{\left|\frac{z_i}{t_{\alpha ,\beta }\left(y\right)}\right|}^{\theta }\frac{1}{\left|z_i\right|}{\pi }_0\left(\alpha ,\beta ,\mu ,\sigma \right)dy.$ (7)

where $\theta =\alpha /\left(\alpha -1\right),Z_i=\left(X_i-\mu \right)/\sigma$, for $i=1,\dots ,n,\alpha \in (0,2],\beta \in \left[-1,1\right],\mu \in \left(-\infty ,+\infty \right),\sigma \in \left(0,+\infty \right)$.

We assume a uniform U(-0.5, 0.5) prior distribution h(y_i) for the latent random quantities Y_i for i=1, 2, . . ., n. Then the joint posterior probability distribution for α, β, µ, σ and y is given by

$\pi \left(\alpha ,\beta ,\mu ,\sigma ,y|x\right)\propto {\left(\frac{\alpha }{\left|\alpha -1\right|\sigma }\right)}^{n}exp\left\{-\sum _{i=1}^n{\left|\frac{z_i}{t_{\alpha ,\beta }\left(y_i\right)}\right|}^{\theta }\right\}\prod _{i=1}^n{\left|\frac{z_i}{t_{\alpha ,\beta }\left(y_i\right)}\right|}^{\theta }\frac{1}{\left|z_i\right|}\prod _{i=1}^{n}h\left(y_i\right){\pi }_0\left(\alpha ,\beta ,\mu ,\sigma \right).$ (8)

The Gibbs sampler can be implemented through the Openbugs software by the following procedure. Firstly, start with the initial values α ⁰, β ⁰, µ ⁰, σ ⁰. For each observation x _i , then y _i is generated from $\pi \left(y_i|\alpha ,\beta ,\mu ,\sigma ,x_i\right)$, for i=1, 2, . . . , n. After generating the vector $\mathit{y}=\left(y_1,y_2,\dots ,y_n\right)$ then the parameters α, β, µ and σ are generated from the conditional distributions $\pi \left(\alpha /\beta ,\mu ,\sigma ,\mathit{x},\mathit{y}\right),\pi \left(\beta /\alpha ,\mu ,\sigma ,\mathit{x},\mathit{y}\right),\pi \left(\mu /\alpha ,\beta ,\sigma ,\mathit{x},\mathit{y}\right)$ and $\pi \left(\sigma /\alpha ,\beta ,\mu ,\mathit{x},\mathit{y}\right)$, respectively.

In this study we assume uniform prior distributions on the domain of each parameter, that is, $\alpha ~U\left(0,2\right),\beta ~U\left(-1,1\right),\sigma ~U\left(0,a\right),\mu ~U\left(-b,b\right)$, where the symbol ~ denotes that the parameter has the distribution, a and b are known hyperparameters, usually very large to have approximately non-informative prior distributions. Other possibility in the data analysis, is to assume a gamma prior distribution for the scale parameter σ and a normal prior distribution for the location parameter µ, but it was observed in the examples considered in this study that, assuming uniform priors for all parameters, the obtained posterior inferences were better, especially in terms of convergence of the simulation algorithms.

With this choice of priors, the simulated Gibbs sampling for the joint posterior distribution is carried out to obtain the posterior summaries of interest.

4 A BAYESIAN ANALYSIS ASSUMING REGRESSION MODELS FOR STABLE DISTRIBUTIONS

In this section, let us assume the presence of a vector $\mathit{v}=\left(v_1,v_2,\dots ,v_k\right)$ of covariates associated to the variable X_i, i = 1,. . . , n, under a stable distribution ${X_i}^{~}{S}_{\alpha }\left(\beta ,{\mu }_i,{\sigma }_i\right)$, which implies that $Z_i=\left(X_i-{\mu }_i\right)/{{\sigma }_i}^{~}{S}_{\alpha }\left(\beta ,0,1\right)$, where the location and the scale parameters µ _i and σ _i are assumed to be dependent of the covariates as follows:

$$\begin{gathered} {\mu }_i=d_0+d_1{v}_{1i}{d}_2{v}_2+\dots +d_k{v}_{ki} \\ \text{and} \\ {\sigma }_i=exp\left(e_0+e_1{v}_{1i}+e_2{v}_{2i}+\dots +e_k{v}_{ki}\right), \end{gathered}$$ (9)

i=1, 2, . . ., n, with $\mathit{d}=\left(d_0,d_1,\dots ,d_k\right)$ and $\mathit{e}=\left(e_0,e_1,\dots ,e_k\right)$ as the regression coefficients.

Assuming a joint prior distribution π ₀ (α, β, d, e ) for α, β , d and e , Buckle (1995) shows that the joint posterior distribution for the parameters α, β , d and e, is given by

$\pi \left(\alpha ,\beta ,\mathit{d},\mathit{e}|\mathit{x}\right)\propto \int {\left(\frac{\alpha }{\left|\alpha -1\right|{\sigma }_i}\right)}^{n}exp\left\{\sum _{i=1}^n{\left|\frac{z_i}{t_{\alpha ,\beta }\left(Y_i\right)}\right|}^{\theta }\right\}\prod _{i=1}^n\left|\frac{z_i}{t_{\alpha ,\beta }\left(y_i\right)}\right|\frac{1}{\left|z_i\right|}{\pi }_0\left(\alpha ,\beta ,\mathit{d},\mathit{e}\right)dy,$ (10)

where $\theta =\alpha /\left(\alpha -1\right),Z_i=\left(X_i-{\mu }_i\right)/{\sigma }_i$, for $i=1,\dots ,n,\alpha \in (0,2],\beta \in \left[-1,1\right],\mathit{x}=\left(x_1,x_2,\dots ,x_n\right)$ and $\mathit{y}=\left(y_1,y_2,\dots ,y_n\right)$ are, respectively, the observed and non-observed data vectors. Note that when α = 2, we have θ = 2 and b_α,β = 0. In this case, we have a Gaussian distribution with mean equal to µ and variance equal to 2σ ².

For a Bayesian analysis of the proposed model, we assume uniform U(a,b) prior distributions for the parameters α, β , d ₀ , d ₁ , d ₂ , . . ., d _k , e ₀ , e ₁ , e ₂ , . . ., e _k where the hyperparameters a and b are assumed to be known in each application following the restrictions $\alpha \in (0,2]$ and $\beta \in \left[-1,1\right]$. We further assume prior independence among all parameters.

In the simulation algorithm to obtain a Gibbs sample for the parameters α, β , d and e, having the joint posterior distribution (10), it is also assumed uniform distributions for the latent variable Y defined in the interval (-0.5, 0.5). With this choice of prior distributions, it is used the standard free Openbugs software to get the posterior summaries of interest, which gives great simplification to obtain the simulated Gibbs samples for the joint posterior distribution of interest.

From the posterior distribution (10), the joint posterior probability distribution for α, β , d, e and $\mathit{y}=\left(y_1,y_2,\dots ,y_n\right)$ is given by,

$\pi \left(\alpha ,\beta ,\mathit{d},\mathit{e},\mathit{y}|\mathit{x}\right)\propto {\left(\frac{\alpha }{\left|\alpha -1\right|{\sigma }_i}\right)}^{n}exp\left\{\sum _{i=1}^n\left|\frac{z_i}{t_{\alpha ,\beta }\left(y_i\right)}\right|\right\}\prod _{i=1}^n{\left|\frac{z_i}{t_{\alpha ,\beta }\left(y_i\right)}\right|}^{\theta }\frac{1}{\left|z_i\right|}·\prod _{i=1}^{n}h\left(y_i\right){\pi }_0\left(\alpha ,\beta ,\mathit{d},\mathit{e}\right),$

where θ and t_α,β (y) are respectively defined in (7) and (8) and the variable Y is defined in the interval (-0.5, 0.5).

Since we are using the Openbugs software to simulate samples for the joint posterior distribution we do not present all full conditional distributions needed for the Gibbs sampling algorithm. This software only requires the data distribution and prior distributions of the interested random quantities. This gives great computational simplification to get the posterior summaries of interest as shown in the applications as follows.

5 APPLICATIONS WITH REAL DATA

In this section two examples with real data sets showing the usefulness of the proposed methodology are introduced.

5.1 Effects of smoking on health

The literature presents many studies on the effects of smoking on health. In this direction, Tager et.al (1983) reported analyses of a study aimed at assessing children’s breathing in the absence or presence of smoking cigarettes, as well as exposure to passive smoke from at least one parent (n = 654 individuals). These studies present some of the earliest attempts to document the obvious signs of reduced pulmonary function from smoking and from exposure to second-hand smoke. The data in this investigation comes from an observational study where the subjects self-select which group they think they belong, to the smoking or non-smoking group. Subjects also selfreport smoking status. Using a spirometer, FEV (Forced Expiratory Volume) is recorded for each subject.

The following variables are reported associated to each individual: age (age in years), FEV (forced expiratory volume in litres), ht (height of subject in cm), gender (female (0), male (1)), smoker (non-smoker (0), smoker (1)). Figure 2 shows the histograms and normal plots of the response variable FEV in the original and logarithmic scales. From these plots, it is observed that the normality is improved in the logarithmic scale. Not considering the logarithm transformation, the distribution is approximately skewed (Figure 2). In this way, we assume a stable distribution to have better fit for the data in the original scale.

Figure 2
Histograms (original and transformed scales).

As a first statistical analysis, it is assumed a linear regression model with normal errors under a classical approach (least squares estimators-LSE) relating the response FEV in the logarithm scale with the covariates age, height, gender and smoker. The fitted regression model is given by

$\log \left(FEV\right)=-1.9414+0.016819\text{height}+0.02363\text{age}+0.0288\text{gender}-0.0471\text{smoker}.$ (12)

Table 1 shows the summaries for LSE analysis (S = 0.145882; R-sq = 80.96%). From the results presented in Table 1, it is observed that assuming a linear multiple regression model with normal errors (response in logarithmic scale) all covariates, height, age, gender and smoking, show significative effects on the response FEV since the p-values are smaller than 0.05 (assumed significance level).

Table 1
Linear normal regression model (FEV data).

	LSE	SE	T-Value	P-Value
Constant	-1.9414	0.0790	-24.58	< 0.001
Age	0.0236	0.0033	7.04	< 0.001
gendercoded	0.0288	0.0118	2.45	0.015
smokercoded	-0.0471	0.0210	-2.25	0.025
Height	0.0168	0.0006	25.34	< 0.001

Older people have larger FEV; males have larger FEV when compared to females; smokers have smaller FEV when compared to non-smokers and larger height implies in larger FEV. Figure 3 presents the residual analysis, from where it is observed that the needed assumptions (normality and constant variance of the residuals) assuming the logarithmic scale for the responses (a lognormal distribution for FEV) are only approximately verified.

Figure 3
Residual plots (normal linear regression model-FEV data).

5.1.1 A Bayesian approach using a stable distribution not considering the presence of covariates

Now, let us assume a stable distribution for the dataset FEV in the original scale not considering the presence of covariates, assuming the following prior distributions for the parameters α, β , µ and σ of the stable distribution: ${\alpha }^{~}U\left(1,2\right),{\beta }^{~}U\left(-1,1\right),{\mu }^{~}U\left(-2,2\right)$ and ${\sigma }^{~}U\left(0,2\right)$, where U(a, b) denotes an uniform distribution in the interval (a, b). Using informative priors for the parameters of the model, some information on the parameters of the stable distribution could be obtained from the histograms presented in Figure 2, as close symmetry leading to prior distributions that should be concentrated in α = 2. In this way, it is considered a uniform prior for the parameter α concentrated in the interval (1, 2). Table 2 shows the posterior summaries of interest (burn-in sample = 100,000 and other 400,000 Gibbs samples taking every 100^th sample) obtained using the Openbugs software. Convergence of the MCMC algorithm was verified from standard traceplots of the simulated samples. The Openbugs code is presented in ^{Appendix 1}.

From the results of Table 2, it is possible to see that the posterior mean for the parameter α has a Monte Carlo estimator based on the simulated Gibbs samples close to the value 2 (Monte Carlo estimator of the posterior mean equal to 1.912) indicating that the fitted stable distribution is close to the usual normal distribution. The convergence of the simulated Gibbs samples was verified from trace plots of the generated samples for each parameter.

Table 2
Bayesian analysis FEV data no covariates.

	mean	s.d.	Lower 95%	Upper 95%
α	1.912	0.0467	1.83	1.993
β	-0.0755	0.4961	-0.9024	0.9151
µ	1.997	0.0024	1.991	2.0
σ	0.7389	0.0261	0.6908	0.7926

5.1.2 A Bayesian approach assuming a stable distribution in presence of covariates

For a second statistical analysis now considering the presence of a vector of covariates, it is assumed the regression model for the location and scale parameters of the stable distribution defined by (9) under a Bayesian approach (use of MCMC methods) relating the location and scale parameters µ and σ of the stable distribution for the response FEV with the covariates age, height, gender and smoke, that is,

$$\begin{gathered} {\mu }_i=d_0+d_{1}heigh{t}_i+d_{2}ag{e}_i+d_{3}gende{r}_i+d_{4}smoke{r}_i \\ {\sigma }_i=exp\left(e_0+e_{1}heigh{t}_i+e_{2}ag{e}_i+e_{3}gende{r}_i+e_{4}smoke{r}_i\right). \end{gathered}$$ (13)

For a Bayesian analysis, we assume an uniform U(1.5, 2) prior distribution for α, an uniform U(-0.5, 0.5) prior distribution for β and uniform prior U(-1, 1) distributions for the regression parameters d₀, d₁, d₂, d₃, d₄, e₀, e₁, e₂, e₃ , e₄. Table 3 shows the posterior summaries of interest.

Table 3
Bayesian analysis - (FEV data presence of covariates).

	Mean	s.d.	Lower 95%	Upper 95%
α	1.975	0.0215	1.923	1.999
β	-0.0138	0.2831	-0.4783	0.4713
d0	2.015	0.0157	2.0	2.061
d1	-0.0119	3,11E-01	-0.0124	-0.0112
d2	0.2394	0.0049	0.2292	0.2479
d3	0.2355	0.0426	0.1491	0.3178
d4	0.0499	0.1060	-0.1615	0.2431
e0	-0.8975	0.1053	-0.9978	-0.613
e1	-0.0066	0.0013	-0.0095	-0.0044
e2	0.0854	0.0156	0.0549	0.1167
e3	0.3526	0.0581	0.2427	0.4694
e4	0.1683	0.1060	-0.0281	0.3937

From the results presented in Table 3, it is observed that assuming the stable distribution in presence of the covariates, the covariate smoker do not show significative effect on the location and scale parameters of the stable distribution assumed for the response FEV (Forced Expiratory Volume) since zero is inside the 95% credible intervals for the corresponding regression parameters d₄ and e₄. All other covariates (age, height and gender) show significative effects on the location and scale parameters of the stable distribution since the zero value is not included in the respective 95% credible intervals for the regression parameters d₁, d₂, d₃, e₁, e₂ and e₃. The convergence of the simulated Gibbs samples was verified from trace plots of the generated samples for each parameter.

Figure 4 shows the scatter plots of the response FEV versus each covariate, from where it is possible to see the great linear dependence between the response FEV and the covariates age and height.

Figure 4
Scatter plots of the response FEV versus each covariate.

5.2 Association between hemoglobin (HbA1c) levels and some covariates for type 2 diabetes mellitus (T2DM)

As another illustrative example, we assume in this example, a data set related to the association between some covariates for type 2 diabetes mellitus (T2DM) patients with different glycated hemoglobin (HbA1c) levels (data set introduced by Shu et al., 2017). The sample size consists of n = 154 patients. In this application, we consider some covariates related to the response glycated hemoglobin (HbA1c) for each patient: age, education in years,duration diagnosis in months, medication compliance, barrier score, dietary knowledge score, BMI (body mass score), gender and total HEI (healthy eating index) score. Figure 5 shows the scatter plots of the response HbA1c average versus each covariate.

Figure 5
Scatter plots of the response HbA1c versus each covariate.

A first statistical analysis is considered assuming a classical multiple linear regression model in presence of the covariates (age, education in years, duration diagnosis in months, medication compliance, barrier score, dietary knowledge score, BMI (body mass score), gender and total HEI (healthy eating index) score. Table 4 shows the summaries for the LSE analysis. The fitted regression model (use of Minitab^® software) is given by

$HBA1c.average=12.37+0.0188age+0.0010educ.years+0.00221duration.diagnosis-0.0263medication.compliance-0.0100barrier.score-0.1267dietary.knowledge.score+0.0581BMI+0.057gender-0.0221total.HEI.score$ (14)

From the results of Table 4, it is observed that assuming a multiple linear regression model with normal errors, the covariates medication compliance, dietary knowledge skill score and BMI show significative effects on the response FEV since the respective p-values are smaller than 0.05 (assumed significance level). Figure 6 shows the residual analysis, from where it is observed that the needed assumptions (normality and constant variance of the residuals) assuming the logarithmic scale for the responses are not well verified, especially the normality of the residuals.

Table 4
Linear normal regression model (HbA1c data).

	Coef	SE	T-value	P-Value
Constant	12.37	2.39	5.18	< 0.001
Age	0.0188	0.0207	0.91	0.366
Educ years	0.0010	0.0459	0.02	0.983
Durationdiagnosis in months	0.0022	0.0014	1.52	0.131
Medication compliance	-0.0263	0.0127	-2.07	0.040
Barrier score	-0.0100	0.0166	-0.60	0.549
Dietary knowledge skill score	-0.1267	0.0446	-2.84	0.005
BMI	0.0581	0.0289	2.01	0.046
Gender	0.057	0.387	0.15	0.883
total HEI score	-0.0221	0.0179	-1.23	0.220

Figure 6
Residual plots (normal linear regression modelHbA1c data).

5.2.1 A Bayesian approach using a stable distribution not considering the presence of covariates

In this subsection, let us assume a stable distribution for the dataset, first not considering the presence of covariates assuming the same uniform prior distributions for the parameters α, β , µ and σ assumed in section 5.1.1. Table 5 presents the posterior summaries of interest (burn-in sample = 300,000 and other 511,000 Gibbs samples taking every 100th sample) obtained using the Openbugs software. Convergence of the MCMC algorithm was verified from standard traceplots of the simulated samples. The convergence of the simulated Gibbs samples was verified from trace plots of the generated samples for each parameter.

Table 5
Bayesian analysis stable distribution not considering the presence of covariates (HbA1c data).

	mean	s.d.	Lower 95%	Upper 95%
α	1.65	0.1424	1.389	1.909
β	-0.2276	0.7033	-0.9878	0.975
µ	8.791	0.1734	8.484	9.153
σ	1.393	0.103	1.208	1.616

5.2.2 A Bayesian approach assuming a stable distribution in presence of covariates

Now, considering the presence of a vector of covariates (age, education in years, duration diagnosis in months, medication compliance, barrier score, dietary knowledge score, BMI, gender and total HEI, it is assumed the regression models for the location and scale parameters of the stable distribution defined by (9) under a Bayesian approach. That is,

$$\begin{gathered} {\mu }_i=d_0+d_{1}ag{e}_i+d_{2}educ.year{s}_i+d_{3}duration.diagnosi{s}_i+d_{4}medication.complianc{e}_i+d_{5}barrier.scor{e}_i+d_{6}dietary.knowlwdge.scor{e}_i+d_{7}BM{I}_i+d_{8}total.HEI.scor{e}_i+d_{9}gende{r}_i \\ {\sigma }_i=exp\left(e_0+e_{1}ag{e}_i+e_{2}educ.year{s}_i+e_{3}duration.diagnosi{s}_i+e_{4}medication.complianc{e}_i+e_{5}barrier.scor{e}_i+e_{6}dietary.knowledge.scor{e}_i+e_{7}BM{I}_i+e_{8}total.HEI.scor{e}_i+e_{9}gende{r}_i\right) \end{gathered}$$ (15)

For a Bayesian analysis, it is assumed an uniform U(0.5,2) prior for α , an uniform U(-1, 1) prior for β , an uniform U(-10,10) prior for the regression parameter d ₀, an uniform prior U(-1, 1) for the regression parameters d₉ and e ₀ and uniform prior distributions U(-0.5, 0.5) for the other regression parameters. Table 6 presents the posterior summaries of interest.

From the results presented in Table 6, it is observed that assuming the stable distribution in presence of the covariates, the covariate duration of diagnosis, dietary knowledge score and BMI show significative effect on the location parameter of the stable distribution assumed for the response HbA1c data since zero is not inside the 95% credible intervals for the corresponding regression parameters d₃, d₆ and d₇. Also the covariate age shows significative effect in the scale parameter since zero is not inside the 95% credible interval for the corresponding regression parameter e₁. The convergence of the simulated Gibbs samples was also verified from trace plots of the generated samples for each parameter.

Table 6
Bayesian analysis - stable regression model (HbA1c data).

	mean	s.d.	Lower 95%	Upper 95%
α	1.82	0.1479	1.47	1.994
β	0.0196	0.6218	-0.9489	0.9678
d0	8.92	0.2138	8.506	9.3380
d1	0.0267	0.0186	-0.0098	0.0629
d2	-4,75E-01	0.0413	-0.0804	0.0879
d3	0.0021	0.0014	4.09E-4	0.0059
d4	-0.0164	0.0129	-0.0404	0.0121
d5	-0.0428	0.0228	-0.0897	0.0014
d6	-0.1468	0.0404	-0.2229	-0.0646
d7	0.0734	0.0309	0.0100	0.1313
d8	-0.0251	0.0179	-0.0603	0.0093
d9	0.2887	0.3430	-0.4444	0.91
e0	0.192	0.0995	0.0026	0.3804
e1	0.0157	0.0080	0.0031	0.0316
e2	-0.0213	0.0193	-0.0592	0.0163
e3	-6,85E-01	0.0010	-0.0027	0.0012
e4	0.0014	0.0058	-0.0113	0.0117
e5	-0.0116	0.0067	-0.0246	0.0013
e6	-0.0384	0.0206	-0.0813	0.0020
e7	0.0155	0.0126	-0.0105	0.0393
e8	-0.0079	0.0081	-0.0233	0.0098
e9	0.2315	0.1399	-0.0732	0.4706

6 CONCLUDING REMARKS

The use of stable distributions can be a good alternative to analyze data from different application areas when the data indicate in some way (descriptive analyses, theoretical interpretations) that the normality assumption usually necessary for several traditional statistical techniques as ANOVA models, multiple linear regression models, among many others, is not satisfied. The stable distribution presents great fit flexibility for asymmetric or heavy tail data sets and good robustness to the presence of outliers, as pointed out by Achcar et al. (2013a,b). For the situation of non-normality of the data, usually data analysts (especially non statisticians) use non-parametric techniques that often do not respond adequately to the objectives of a research, especially in the medical area. The literature presents several alternatives to analyze data when there is no indication of normality, such as the use of appropriate distributions for heavy tails or finite mixtures of parametric distributions, among several other possibilities, but the use of stable distributions could be a good alternative in the data analysis. The two examples with medical data presented in the study show that stable distributions can be very useful in clinical studies, in addition to what is already considered in other areas of application, such as in particular, the area of finance where the use of stable distributions is already very popular.

The results of this study show, especially for regression models, that the assumption of errors with stable distribution could be considered in medical studies, where the computational costs for the determination of Bayesian estimators are not very high using MCMC simulation methods and data augmentation techniques, especially using existing Bayesian software to simulate samples of the joint posterior distribution of interest like the Openbugs software used in this study.

We also point out that the use of the proposed methodology could be useful to check the normality of the errors in multiple regression parameters since the parameter α of the stable distribution indicates a Gaussian distribution if α = 2, that is, if the value 2 is inside a 95% credible interval for α, there is an indication of normality for the residuals.

Although the stable distribution probability density function does not have a closed form, we obtained satisfactory results for parameter estimation in the presence of a vector of covariates by applying the Bayesian inference method based on data augmentation techniques.

In both examples presented in Section 5, it was shown that using joint regression models for the location and scale parameters of the stable distribution under a Bayesian approach and using existing MCMC methods assuming uniform priors for all parameters, good inferences were obtained for the model parameters. It is important to note that using a simple linear regression model for the location parameter of the stable distribution, assuming normal distributions as prior distributions for the regression parameters, it was not possible to obtain good convergence results from the MCMC simulation algorithm even after a large number of iterations.

We also showed that the MCMC methods for a Bayesian analysis of these models using the Openbugs software give a great simplification to obtain very accurate posterior summaries of interest, not requiring large computational time to perform, even when the simulation of a large number of Gibbs samples is necessary for the convergence of the algorithm.

Other important issue to use Bayesian methods is the choice of prior distributions for the parameters of the proposed model. In the two applications with medical data, we have used noninformative prior distributions since we did not have prior opinion of medical experts for the elicitation of priors. The data sets were obtained from the literature to illustrate the proposed methodology. In practical work, usually we can use informative priors based on the information of medical experts leading to more accurate inference results. This is a great advantage of the use of Bayesian methods in applications in medical studies and applied work in general.

References

• ACHCAR JA, LOPES SRC, MAZUCHELI J & LINHARES R. 2013a. A Bayesian Approach for Stable Distributions: Some Computational Aspects. Open Journal of Statistics, 3: 268-277.
• ACHCAR JA, ACHCAR A & MARTINEZ EZ. 2013b. Robust Linear Regression Models: Use of a Stable Distribution for the Response Data. Open Journal of Statistics, 3(6): 409-416.
• BERNARDO JM & SMITH AFM. 2007. Bayesian Theory, 2nd edition. Wiley Series in Probability and Statistics. Wiley, New York.
• BORAK S, HARDLE WK & WERON R. 2005. Stable Distributions. Statistical Tools for Finance and Insurance. 10.1007/3-540-27395-61.
  » https://doi.org/10.1007/3-540-27395-61
• BUCKLE DJ. 1995. Bayesian inference for stable distributions. Journal of the American Statistical Association, 90(430): 605-613.
• GELMAN A, CARLIN JB, STERN HS, DUNSON DB, VEHTARI A & RUBIN DB. 2013. Bayesian data analysis. London, United Kingdom: Chapman and Hall/CRC.
• GNEDENKO BV & KOLMOGOROV AN. 1954. Limit Distributions for Sums of Independent Random Variables. Cambridge, MA: Addison-Wesley.
• KARLING MJ, LOPES SRC & SOUZA RM. 2021. A Bayesian approach for estimating the parameters of an α-stable distribution. Journal of Statistical Computation and Simulation, 91(9):1713-1748.
• KATEREGGA M. 2017. Parameter estimation for stable distributions with application to commodity futures log-returns. Cogent Economics & Finance, 5(1):1318813.
• LEMKE T, RIABIZ M & GODSILL SJ. 2015. Fully Bayesian inference for α-stable distributions using a Poisson series representation, Digital Signal Processing, 47(C), 96-115.
• LOMBARDI MJ. 2004. Bayesian inference for alpha-stable distributions: a random walk MCMC approach, Econometrics Working Papers Archive wp2004 11, Universita’ degliStudi di Firenze, Dipartimento di Statistica, Informatica, Applicazioni “G. Parenti”. https://ideas.repec.org/s/fir/econom.html
  » https://ideas.repec.org/s/fir/econom.html
• O’HAGAN A & FORSTER JJ. 2004. Bayesian Inference, 2nd edition, v. 2B of Kendall’s Advanced Theory of Statistics. Arnold, London.
• ORAL E & ERDEMIR C. 2012. A Bayesian Estimation of Stable Distributions, Journal of Statistical and Econometric Methods, 1(3):39-52.
• SHU PS, CHAN YM & HUANG SL. 2017. Higher body mass index and lower intake of dairy products predict poor glycaemic control among type2 diabetes patients in Malaysia. PLoS One, 12:e0172231.
• SPIEGELHALTER DJ, THOMAS A, BEST NG & LUNN D. 2003. WinBUGS User’s Manual, MRC Biostatistics Unit, Cambridge.
• TAGER I, WEISS SAM, ROSNER B & SPEIZER F. 1983. Longitudinal study of the effects of maternal smoking on pulmonary function. New England Journal of Medicine, 309(12):699-703.
• ZOLOTAREV VM. 1986. One-Dimensional Stable Distributions, Vol 65 of Translations of Mathematical Monographs, American Mathematical Society, Providence, Rhode Island.

APPENDIX 1 (OPENBUGS-APPLICATION FEV DATA OF APPLICATION 4.1)

• Stable distribution for the data (no presence of covariates)

• Stable distribution (regression in the location parameter and scale parameters)

Publication Dates

• Publication in this collection
  10 June 2022
• Date of issue
  2022

History

• Received
  22 July 2021
• Accepted
  14 Mar 2022