Bayesian Reference Analysis for the Generalized Normal Linear Regression Model

Abstract

This article proposes the use of the Bayesian reference analysis to estimate the parameters of the generalized normal linear regression model. It is shown that the reference prior led to a proper posterior distribution, while the Jeffreys prior returned an improper one. The inferential purposes were obtained via Markov Chain Monte Carlo (MCMC). Furthermore, diagnostic techniques based on the Kullback–Leibler divergence were used. The proposed method was illustrated using artificial data and real data on the height and diameter of Eucalyptus clones from Brazil.

1. Introduction

Although the normal distribution is used in many fields for symmetrical data modeling, when the data come from a distribution with lighter or heavier tails, the assumption of normality becomes inappropriate. Such circumstances show the need for more flexible models such as the Generalized Normal (GN) distribution [1], which encompasses various distributions such as the Laplace, normal, and uniform distributions.

In the presence of covariates, the normal linear regression model can be used to investigate and model the relationship between variables, assuming that the observations follow a normal distribution. However, it is well known that a normal linear regression model can be influenced by the presence of outliers [2,3]. In these circumstances, as discussed above, it is necessary to use models in which the error distribution presents heavier or lighter tails than normal such as the GN distribution.

Due to its flexibility, the GN distribution is considered a tool to reduce the impact of the outliers and obtain robust estimates [4,5,6]. This distribution has been used in different contexts, with different parameterizations, but the main difficulties to adopting GN distribution modeling have been computational problems since there are no explicit expressions for estimators of the shape parameter. The estimation of the shape parameter should be done through numerical methods. However, in the classical context, there are problems of convergence, as demonstrated in the studies [7,8].

In the Bayesian context, the methods presented in the literature for the estimation of the parameters of the GN distribution are restricted and applied to particular cases. West [9] proved that a scale mixture of the normal distribution could represent the GN distribution. Choy and Smith [10] used the prior GN distribution for the location parameter in the Gaussian model. They obtained the summaries of the posterior distribution, estimated through the Laplace method, and examined their robustness properties. Additionally, the authors used the representation of the GN distribution as a scale mixture of the normal distribution in random effect models and considered the Markov Chain Monte Carlo method (MCMC) to obtain the posterior summaries of the parameter of interest.

Bayesian procedures for regression models with GN errors have been discussed earlier. Box and Tiao [4] used the GN distribution from the Bayesian approach in which they proposed robust regression models as an alternative to the assumptions of normality of errors in regression models. Salazar et al. [11] considered an objective Bayesian analysis for exponential power regression models, i.e., a reparametrized version of the GDN. They derived the Jeffreys prior and showed that such a prior results in an improper posterior. To overcome this limitation, they considered a modified version of the Jeffreys prior under the assumption of independence of the parameters. This assumption does not hold since the scale and the shape parameters correlate with the proposed distribution. Additionally, the use of the Jeffreys prior is not appropriate in many cases and can cause strong inconsistencies and marginalization paradoxes (see Bernardo [12], p. 41).

Reference priors, also called objective priors, can be used to overcome this problem. This method was introduced by Bernardo [13] and enhanced by Berger and Bernardo [14]. For the proposed methodology, the prior information is dominated by the information provided by the data, resulting in a vague influence of the prior distribution. Estimations are made based on priors through the maximization of the expected Kullback–Leibler (KL) divergence between the posterior and prior distributions. The resulting reference prior affords a posterior distribution that has interesting features, such as consistent marginalization, one-to-one invariance, and consistent sampling properties [12]. Some applications of reference priors can be seen for other distributions in [15,16,17,18,19,20].

In this paper, we considered the reference approach for estimating the parameters of the GN linear regression model. We showed that the reference prior led to a proper posterior distribution, whereas the Jeffreys prior brought an improper one and should not be used. The posterior summaries were obtained via Markov Chain Monte Carlo (MCMC) methods. Furthermore, diagnostic techniques based on the Kullback–Leibler divergence were used.

To exemplify the proposed model, we considered the 1309 entries on the height and diameter of Eucalyptus clones (more details are given in Section 8). For these data, a linear regression model using the normal distribution for the residuals was not adequate, and so, we used the GN linear regression approach.

This paper is composed as follows. Section 2 presents the GN distribution with some special cases, and in the following Section 3, we introduce the GN linear regression model. Section 4 shows the reference and Jeffreys priors, respectively, for a GN linear regression model. Then, the model selection criteria, in Section 6, are discussed, as well as their applications. Section 7 shows the proposed method for analyzing compelling cases considering the Bayesian reference prior approach through the use of the Kullback–Leibler divergence. Section 8 presents studies with an artificial and a real application, respectively. Finally, we discuss the conclusions in Section 9.

2. Generalized Normal Distribution

The Generalized Normal distribution (GN distribution) has been referred to in the literature under different names and parametrizations such as the exponential power distribution or the generalized error distribution. The first formulation of this distribution [1] as a generalization of the normal distribution was characterized by the location, scale, and shape parameters.

Here, we considered the form presented by Nadarajah [21]. It is understood that the random variable Y is the GN distribution given its probability density function (pdf) as:

$$f\left(y\right|\mu ,\tau ,s)=\frac{s}{2\tau \Gamma \left(\frac{1}{s}\right)}exp\left\{-{\left(\frac{\left|y-\mu \right|}{\tau }\right)}^s\right\},y,\mu \in \mathfrak{R},\tau ,s>0.$$

(1)

The parameter $\mu$ is the mean; $\tau$ is the scale factor; and s is the shape parameter. In particular, the mean, variance, and coefficient of the kurtosis of Y are given by:

$$E\left(Y\right)=\mu ,Var\left(Y\right)=\frac{{\tau }^2\Gamma \left(\frac{3}{s}\right)}{\Gamma \left(\frac{1}{s}\right)}\mathtt{and}\gamma =\frac{\Gamma \left(\frac{1}{s}\right)\Gamma \left(\frac{5}{s}\right)}{\Gamma {\left(\frac{3}{s}\right)}^2},$$

respectively.

The GN distribution characterizes leptokurtic distributions if $s<2$ and platykurtic distributions if $s>2$. In particular, the GN distribution displays the Laplace distribution when $s=1$ and the normal distribution when $s=2$ and $\tau$ is equal to $\sqrt{2}\sigma$, where $\sigma$ is the standard deviation, and when $s\to \infty$, the pdf converges to a uniform distribution in (1).

This distribution is flexible-symmetrical concerning the average and unimodality. Moreover, it allows a more flexible fit for the kurtosis than a normal distribution. Furthermore, the ability of the GN distribution to provide a precise fit for the data depends on its shape.

Zhu and Zinde-Walsh [22] proposed a reparameterization of the asymmetric exponential power distribution that allows us to observe the effect of the shape parameter on the distribution, which was adapted for the GN distribution. Using a similar reparameterization, $\sigma =\tau \Gamma \left(1+\frac{1}{s}\right)$, the GN distribution in (1) is given by:

$$f\left(y\right|\mu ,\sigma ,s)=2^{-1}{\sigma }^{-1}exp\left\{-{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y-\mu |}{\sigma }\right)}^s\right\}y,\mu \in \mathfrak{R},\sigma >0,s>0.$$

(2)

The reparametrization above is used throughout the paper. Figure 1 shows the density functions of the GN distribution in (2) for various parameter values.

3. Generalized Normal Linear Regression Model

The GN linear regression model is defined as:

$$Y_i={\mathit{x}}_i^{\top }\mathit{\beta }+{ϵ}_i,i=1,\dots ,n,$$

(3)

where $Y_i$ is the vector of the response for the ith case, ${\mathit{x}}_i^{\top }=(x_{i1},\dots ,x_{ip})$ contains the values of the explanatory variables, $\mathit{\beta }={({\beta }_1,\dots ,{\beta }_p)}^{\top }$ is the vector of regression coefficients, and ${ϵ}_i$ is the vector of random errors that follows a GN distribution with mean zero, scale parameter $\sigma$, and shape parameter s.

Therefore, the likelihood function is:

$$L\left(\mathit{y}\right|\mathit{\beta },\sigma ,s)=2^{-n}{\sigma }^{-n}exp\left\{-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^s\right\}.$$

(4)

The log-likelihood function (4) is given by:

$$logL\left(\mathit{y}\right|\mathit{\beta },\sigma ,s)=-nlog2-nlog\sigma -\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^s.$$

(5)

The first-order derivatives of the log-likelihood function in (5) are given by:

$$\begin{array}{cc}\hfill \frac{\partial logL}{\partial \mathit{\beta }}& =\frac{s\Gamma \left(1+\frac{1}{s}\right)}{\sigma }\sum _{i=1}^n{\mathit{x}}_i^{\top }{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^{s-1}sign(y_i-{\mathit{x}}_i^{\top }\mathit{\beta }),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\partial logL}{\partial \sigma }& =-\frac{n}{\sigma }+\frac{s}{\sigma }\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^s,\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{\partial logL}{\partial s}& =-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\beta |}{\sigma }\right)}^s\left[log\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\beta |}{\sigma }\right)-\frac{\Psi \left(1+\frac{1}{s}\right)}{s}\right],\hfill \end{array}$$

(8)

where $\Psi \left(s\right)=\frac{{\Gamma }^{{}^{\prime }}\left(s\right)}{\Gamma \left(s\right)}$ is the digamma function.

The score function obtains the Fisher information matrix. This matrix is helpful to get the reference priors for the model parameters. The following proposition obtains the elements of the Fisher information matrix for the model in (3).

Proposition 1.

Let $I\left(\mathit{\theta }\right)$ be the Fisher information matrix, with $\mathit{\theta }=(\mathit{\beta },\sigma ,s)$. The elements of the Fisher information matrix,

$$I_{ij}\left(\mathit{\theta }\right)=-E\left(\frac{{\partial }^{2}logf\left(y\right|\mathit{\theta })}{\partial {\theta }_i\partial {\theta }_j}\right)=E\left(\frac{\partial logf\left(y\right|\mathit{\theta })}{\partial {\theta }_i}\frac{\partial logf\left(y\right|\mathit{\theta })}{\partial {\theta }_j}\right),i,j=1,2,3,$$

with $I_{ij}\left(\mathit{\theta }\right)=I_{ji}\left(\mathit{\theta }\right)$ and ${\theta }_j$ the jth element of $\mathit{\theta }$, are given by,

$$\begin{array}{cc}\hfill I_{11}\left(\theta \right)& =E\left[\left(\frac{\partial logL}{\partial \beta }\right)\left(\frac{\partial logL}{\partial {\beta }^{\top }}\right)\right]=\frac{\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right)}{{\sigma }^2}\sum _{i=1}^n{\mathit{x}}_i{\mathit{x}}_i^{\top },\hfill \\ \hfill I_{12}\left(\theta \right)& =E\left[\left(\frac{\partial logL}{\partial \beta }\right)\left(\frac{\partial logL}{\partial \sigma }\right)\right]=\mathbf{0},\hfill \\ \hfill I_{13}\left(\theta \right)& =E\left[\left(\frac{\partial logL}{\partial \beta }\right)\left(\frac{\partial logL}{\partial s}\right)\right]=\mathbf{0},\hfill \\ \hfill I_{22}\left(\theta \right)& =E\left[{\left(\frac{\partial logL}{\partial \sigma }\right)}^2\right]=\frac{ns}{{\sigma }^2},\hfill \\ \hfill I_{23}\left(\theta \right)& =E\left[\left(\frac{\partial logL}{\partial \sigma }\right)\left(\frac{\partial logL}{\partial s}\right)\right]=-\frac{n}{\sigma s},\hfill \\ \hfill I_{33}\left(\theta \right)& =E\left[{\left(\frac{\partial logL}{\partial s}\right)}^2\right]=\frac{n}{s^3}\left\{\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right\},\hfill \end{array}$$

where $s>1$ and ${\Psi }^{{}^{\prime }}\left(s\right)=\frac{\partial \Psi \left(s\right)}{\partial s}$ is the trigamma function. The restriction $s>1$ ensures that the elements $I_{ij}\left(\mathit{\theta }\right)$, calculated for $i,j=1,2,3$, are finite and the information matrix $I\left(\mathit{\theta }\right)$ is positive definite.

For further details of this proof, please see Proposition 5 in Zhu and Zinde-Walsh [22]. Then, Fisher’s information matrix is given by:

$$I\left(\theta \right)=\left[\begin{array}{ccc}\frac{\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right)}{{\sigma }^2}\sum _{i=1}^n{\mathit{x}}_i{\mathit{x}}_i^{\top }& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \frac{ns}{{\sigma }^2}& -\frac{n}{\sigma s}\\ \mathbf{0}& -\frac{n}{\sigma s}& \frac{n}{s^3}\left\{\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right\}\end{array}\right].$$

(9)

The corresponding inverse Fisher’s information matrix is given by:

$$B\left(\theta \right)=\left[\begin{array}{ccc}\frac{{\sigma }^2}{\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right){\sum }_{i=1}^n{\mathit{x}}_i{\mathit{x}}_i^{\top }}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \frac{{\sigma }^2}{ns\left[1-\frac{s}{(1+s){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)}\right]}& \frac{\sigma s^2}{n\left[-s+(1+s){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}\\ \mathbf{0}& \frac{\sigma s^2}{n\left[-s+(1+s){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}& \frac{s^4}{n\left[-s+(1+s){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}\end{array}\right].$$

(10)

The matrix in (9) coincides with the Fisher information matrix found by Salazar et al. [11] due to the one-to-one invariance property.

4. Objective Bayesian Analysis

An important class of objective priors was introduced by Bernardo [13] and later developed by Berger and Bernardo [14]. This class of prior is known as reference priors. A vital feature of the method developed by Berger–Bernardo is the specific treatment given to interest and nuisance parameters. The construction of the reference prior, in the case of nuisance parameters, must be made using an ordered parameterization. The parameter of interest is selected, and the procedure is followed below (see Bernardo [12], for a detailed discussion).

Proposition 2.

(Reference priors under asymptotic normality) Let a probability model $f\left(y\right|\theta ,\mathit{\omega }),y\in {\mathbb{R}}^n,\mathit{\omega }=({\omega }_1,\dots ,{\omega }_m),\theta \in \mathrm{\Phi },\mathit{\omega }\in \mathrm{\Omega }={\prod }_{j=1}^m{\mathrm{\Omega }}_j,$ with $m+1$ real-valued parameters, and let θ be the quantity of interest. For instance, the posterior distribution of the parameters $(\theta ,{\omega }_1,\dots ,{\omega }_m)$ is asymptotically normal with covariance matrix $B(\widehat{\theta },{\widehat{\omega }}_1,\dots ,{\widehat{\omega }}_m)$, where $H=B^{-1},B_j$ is the upper left $j\times j$ submatrix of B, $H_j=B_j^{-1}$, and $h_{jj}(\theta ,{\omega }_1,\dots ,{\omega }_m)$ is the lower right element of $H_j$.

Then, it holds that the conditional reference priors can be represented as:

$$\begin{array}{cc}\hfill {\pi }^{\mathsf{R}}\left({\omega }_m\right|\theta ,{\omega }_1,\cdots ,{\omega }_{m-1})& \propto h_{m+1,m+1}^{\frac{1}{2}}(\theta ,{\omega }_1,\dots ,{\omega }_m),\hfill \end{array}$$

and:

$$\begin{array}{cc}\hfill {\pi }^{\mathsf{R}}\left({\omega }_i\right|\theta ,{\omega }_1,\cdots ,{\omega }_{i-1})\propto exp& \{{\int }_{{\mathrm{\Omega }}_{i+1}}\cdots {\int }_{{\mathrm{\Omega }}_m}log{h}_{i+1,i+1}^{\frac{1}{2}}(\theta ,{\omega }_1,\dots ,{\omega }_m)\times \hfill \\ & \left\{\prod _{j=i+1}^m{\pi }^{\mathsf{R}}\left({\omega }_j\right|\theta ,{\omega }_1,\cdots ,{\omega }_{j-1})\right\}d{\omega }_{i+1}\},\hfill \end{array}$$

where $d{\mathit{\omega }}_j=d{\omega }_j\times \cdots \times d{\omega }_m$ and ${\pi }^{\mathsf{R}}\left({\omega }_i\right|\theta ,{\omega }_1,\cdots ,{\omega }_{i-1}),i=1,\cdots ,m$ are all proper. A compact approximation will be only required, for the corresponding integrals, if any of those conditional reference priors is not proper.

Furthermore, the marginal reference prior of θ is:

$${\pi }^{\mathsf{R}}\left(\theta \right)\propto exp\left\{{\int }_{{\mathrm{\Omega }}_1}\cdots {\int }_{{\mathrm{\Omega }}_m}log{B_{11}}^{-\frac{1}{2}}(\theta ,{\omega }_1,\dots ,{\omega }_m)\left\{\prod _{j=1}^m{\pi }^{\mathsf{R}}\left({\omega }_j\right|\theta ,{\omega }_1,\cdots ,{\omega }_{j-1})\right\}d{\mathit{\omega }}_1\right\},$$

where $B_{11}^{-\frac{1}{2}}(\theta ,{\omega }_1,\dots ,{\omega }_m)=h_{11}^{\frac{1}{2}}(\theta ,{\omega }_1,\dots ,{\omega }_m)$.

The reference posterior distribution associated with $\theta$, after $\mathit{y}$, is given by:

$${\pi }^{\mathsf{R}}\left(\theta \right|\mathit{y})\propto {\pi }^{\mathsf{R}}\left(\theta \right){\int }_{{\mathrm{\Omega }}_1}\cdots {\int }_{{\mathrm{\Omega }}_m}f\left(\mathit{y}\right|\theta ,{\omega }_1,\cdots ,{\omega }_m)\left\{\prod _{j=1}^m{\pi }^{\mathsf{R}}\left({\omega }_j\right|{\omega }_1,\cdots ,{\omega }_{j-1})\right\}d{\omega }_1,\cdots ,d{\omega }_m.$$

The proposition is first presented for the presence of one nuisance parameter and further extended to a vector of nuisance parameters (see Bernardo [12] for a detailed discussion and proofs). Here, as our model has a special structure, we also considered an additional result presented in the following corollary that will be used to construct the reference prior to the GN distribution.

Corollary 1.

If the nuisance parameter spaces ${\mathrm{\Omega }}_i$ do not depend on $\{\theta ,{\omega }_1,\dots ,{\omega }_{i-1}\}$ and the functions $h_{11},h_{22},\cdots ,h_{mm},h_{m+1m+1}$ factorize in the form:

$$h_{11}^{\frac{1}{2}}(\theta ,{\omega }_1,\cdots ,{\omega }_m)=f_0\left(\theta \right)g_0({\omega }_1,\cdots ,{\omega }_m)and$$

$$h_{i+1,i+1}^{\frac{1}{2}}(\theta ,{\omega }_1,\cdots ,{\omega }_m)=f_i\left({\omega }_i\right)g_i(\theta ,{\omega }_i,\cdots ,{\omega }_{i-1},{\omega }_{i+1},\cdots ,{\omega }_m),i=1,\cdots ,m,$$

then ${\pi }^{\mathsf{R}}\left(\theta \right)\propto f_0\left(\theta \right)$ and ${\pi }^{\mathsf{R}}\left({\omega }_i\right|\theta ,{\omega }_1,\cdots ,{\omega }_{i-1})\propto f_i\left({\omega }_i\right),i=1,\cdots ,m,$ and there is no need for compact approximations, even if ${\pi }^{\mathsf{R}}\left({\omega }_i\right|\theta ,{\omega }_1,\cdots ,{\omega }_{i-1})$ are not proper.

Under appropriate regularity conditions (see Bernardo [12]), the posterior distribution of $(\theta ,{\omega }_1,\dots ,{\omega }_m)$ is asymptotically normal with mean $(\widehat{\theta },{\widehat{\omega }}_1,\dots ,{\widehat{\omega }}_m)$, the corresponding MLEs and covariance matrix $B(\widehat{\theta },{\widehat{\omega }}_1,\dots ,{\widehat{\omega }}_m)$, where $B(\widehat{\theta },{\widehat{\omega }}_1,\dots ,{\widehat{\omega }}_m)=I^{-1}(\widehat{\theta },{\widehat{\omega }}_1,\dots ,{\widehat{\omega }}_m)$ and $I(\theta ,{\omega }_1,\dots ,{\omega }_m)$ is the corresponding $(m+1)\times (m+1)$ Fisher information matrix; in that case, $H(\theta ,{\omega }_1,\dots ,{\omega }_m)=I(\theta ,{\omega }_1,\dots ,{\omega }_m)$, and the reference prior may be computed from the elements of Fisher matrix $I(\theta ,{\omega }_1,\dots ,{\omega }_m)$.

4.1. Reference Prior

The parameter vector $(\beta ,\sigma ,s)$ is ordered and divided into three distinct groups, according to their inferential importance. We considered here the case in which $\beta$ is the parameter of interest and $\sigma$ and s are the nuisance parameters. To obtain a joint reference prior for the parameters $\beta ,\sigma$, and s, the following ordered parameterization was adopted:

$${\pi }^{\mathrm{R}}(\beta ,\sigma ,s)={\pi }^{\mathrm{R}}\left(s\right|\mathit{\beta },\sigma ){\pi }^{\mathrm{R}}\left(\sigma \right|\mathit{\beta }){\pi }^{\mathrm{R}}\left(\mathit{\beta }\right).$$

Consider the Fisher matrix in (9), the inverse Fisher matrix in (10)$\left(\right)$, and Corollary 1. Let $H\left(\mathit{\theta }\right)=I\left(\mathit{\theta }\right)$; it follows that $h_{33}^{\frac{1}{2}}(\beta ,\sigma ,s)=\sqrt{\frac{n}{s^3}\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)}=f_2\left(s\right)g_2(\mathit{\beta },\sigma ).$ Then,

$${\pi }^{\mathsf{R}}\left(s\right|\mathit{\beta },\sigma )\propto s^{-\frac{3}{2}}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}^{\frac{1}{2}}.$$

Let $H_2\left(\mathit{\theta }\right)=B_2^{-1}\left(\mathit{\theta }\right)$, where $B_2\left(\mathit{\theta }\right)$ is the upper left $2\times 2$ submatrix of $B\left(\mathit{\theta }\right)$; it follows that $h_{22}^{\frac{1}{2}}(\mathit{\beta },\sigma ,s)=\frac{1}{\sigma }\sqrt{ns\left[1-\frac{s}{(1+s){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)}\right]}=f_1\left(\sigma \right)g_1(\mathit{\beta },s).$ Then,

$${\pi }^{\mathsf{R}}\left(\sigma \right|\mathit{\beta })\propto \frac{1}{\sigma }.$$

Finally, let $h_{11}(\mathit{\beta },\sigma ,s)=B_{11}^{-1}(\mathit{\beta },\sigma ,s)$; it follows that $h_{11}^{\frac{1}{2}}(\mathit{\beta },\sigma ,s)=\sqrt{{\sum }_{i=1}^n{\mathit{x}}_i{\mathit{x}}_i^{\top }}\sqrt{\frac{\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right)}{{\sigma }^2}}=f_0\left(\mathit{\beta }\right)g_0(\sigma ,s)$. Then,

$${\pi }^{\mathsf{R}}\left(\mathit{\beta }\right)\propto 1.$$

Therefore, a joint reference prior for the ordered parameter is given by:

$${\pi }^{\mathsf{R}}(\mathit{\beta },\sigma ,s)\propto \frac{1}{\sigma }\times s^{-\frac{3}{2}}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}^{\frac{1}{2}},$$

(11)

where $\mathit{\beta }\in {\mathfrak{R}}^p,\sigma \in {\mathfrak{R}}^{+}$, and $s>1$.

Using the likelihood function (4) and the joint reference prior (11), we obtain the joint posterior distribution for $\mathit{\beta },\sigma$, and s,

$${\pi }^{\mathsf{R}}(\mathit{\beta },\sigma ,s|\mathit{y})\propto {\sigma }^{-(n+1)}{s}^{-\frac{3}{2}}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}^{\frac{1}{2}}exp\left\{-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_{\mathit{i}}^{\top }\mathit{\beta }|}{\sigma }\right)}^s\right\}.$$

(12)

The posterior conditional probability densities are given by,

$${\pi }^{\mathsf{R}}\left(\mathit{\beta }\right|\sigma ,s,\mathit{y})\propto exp\left\{-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^s\right\},$$

(13)

$${\pi }^{\mathsf{R}}\left(\sigma \right|\mathit{\beta },s,\mathit{y})\propto {\sigma }^{-(n+1)}exp\left\{-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^s\right\},$$

(14)

$${\pi }^{\mathsf{R}}\left(s\right|\mathit{\beta },\sigma ,\mathit{y})\propto s^{-\frac{3}{2}}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}^{\frac{1}{2}}exp\left\{-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{s}\right)|y_i-{\mathit{x}}_i^{\top }\mathit{\beta }|}{\sigma }\right)}^s\right\}.$$

(15)

The densities in (13) and (15) do not belong to any known parametric family, and the densities in (14) can be easily reduced to an inverse-gamma distribution form by the transformation $\lambda ={\sigma }^s$. The parameters of interest were obtained by Monte Carlo methods in a Markov Chain (MCMC). Thus, the posterior densities were evaluated by applying the Metropolis–Hastings algorithm; see, e.g., Chen et al. [23].

4.2. A Problem with the Jeffreys Prior

The use of the Jeffreys prior in the multiparametric case is often controversial. Bernardo ([12], p. 41) argued that the use of the Jeffreys prior is not appropriate in many cases and can cause strong inconsistencies and marginalization paradoxes. This prior is obtained from the square root of the determinant of the Fisher information matrix of (9),

$${\pi }^{\mathrm{J}}(\mathit{\beta },\sigma ,s)\propto \sqrt{det\left(I\right(\theta \left)\right)}=\sqrt{det\left(I_{11}\right)[I_{22}{I}_{33}-I_{23}^2]},$$

where:

$$det\left(I_{11}\right)={\left[\frac{\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right)}{{\sigma }^2}\right]}^{p}det\left(\sum _{i=1}^n{\mathit{x}}_{\mathit{i}}{\mathit{x}}_{\mathit{i}}^{\top }\right),$$

and:

$$[I_{22}{I}_{33}-I_{23}^2]=\frac{n^2}{{\sigma }^2{s}^2}\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)-1\right],$$

is given by:

$${\pi }^{\mathrm{J}}(\mathit{\beta },\sigma ,s)\propto {\sigma }^{-(p+1)}{\left[\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right)\right]}^{\frac{p}{2}}{s}^{-1}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)-1\right]}^{\frac{1}{2}}.$$

(16)

Such a prior was also presented in Salazar et al. [11]. Both priors are from the family of prior distributions represented as:

$$\pi (\mathit{\beta },\sigma ,s)\propto \frac{\pi \left(s\right)}{{\sigma }^c},a\in \mathfrak{R};$$

(17)

they are usually improper, and c is a hyperparameter, while $\pi \left(s\right)$ is the prior related to the shape parameter.

The Jeffreys prior and reference prior are of the form (17) with, respectively,

$$c=1\mathrm{and}{\pi }^{\mathsf{R}}\left(s\right)\propto s^{-\frac{3}{2}}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)\right]}^{\frac{1}{2}},$$

(18)

$$c=p+1\mathrm{and}{\pi }^{\mathrm{J}}\left(s\right)\propto s^{-1}{\left[\left(1+\frac{1}{s}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{s}\right)-1\right]}^{\frac{1}{2}}{\left[\Gamma \left(\frac{1}{s}\right)\Gamma \left(2-\frac{1}{s}\right)\right]}^{\frac{p}{2}}.$$

(19)

The posterior distribution associated with the prior in (17) is proper if:

$${\int }_1^{\infty }L\left(s\right|\mathit{y})\pi \left(s\right)ds<\infty ,$$

(20)

where $L\left(s\right|\mathit{y})$ is the integrated likelihood for s, given by:

$${\int }_{{\mathfrak{R}}^p}{\int }_0^{\infty }L(\mathit{\beta },\sigma ,s|\mathit{y}){\sigma }^{-c}d\sigma d\mathit{\beta }.$$

(21)

Corollary 2.

The marginal prior for s given in Equations (18) and (19) is a continuous function in $[1,\infty ]$, and as $s\to \infty$, we have ${\pi }^{\mathsf{R}}\left(s\right)=O\left(s^{-3/2}\right)$ and ${\pi }^{\mathrm{J}}\left(s\right)=O\left(s^{\frac{p-2}{2}}\right)$.

Proof. 

See Appendix A. □

Corollary 3.

For $n>p+1-c$, the likelihood function for the parameter s, $L\left(s\right|\mathit{y})$, under the class of priors (17), is a continuous function bounded in $[1,\infty ]$ and with complexity $O\left(1\right)$ when $s\to \infty$.

Proof. 

See Appendix A. □

Proposition 3.

The reference prior given in (11) yields a proper posterior distribution, and the Jeffreys prior given in (16) leads to an improper posterior distribution.

Proof. 

See Appendix A. □

Therefore, the Jeffreys prior leads to an improper posterior distribution and cannot be used in a Bayesian analysis. Another objective prior known as the maximum data information prior could be considered [16,24,25,26]; however, such a prior is not invariant under one-to-one transformations, which limits its use. Additionally, the main aim was to consider objective priors. We avoided the use of normal or gamma priors due to the lack of invariance in the parameters. Moreover, the prior may depend on hyperparameters that are not easy to elicit for the GN distributions, and the posterior estimates may vary depending on the included information. Finally, Bernardo [12] pointed out that the use of our “flat” priors to assign “non-informative” priors should be strongly discouraged because they often result in the suppression of important inappropriate and unjustified assumptions that can easily have a strong influence on the analysis, or even make it invalid.

5. Metropolis–Hastings Algorithm

Here, the Metropolis–Hastings algorithm is considered to sample from $\mu ,\sigma$, and s. In this case, the following conditional distributions are considered: ${\pi }^{\mathsf{R}}\left(\mu \right|\sigma ,s,\mathbf{y})$, ${\pi }^{\mathsf{R}}\left(\sigma \right|\mu ,s,\mathbf{y})$, and ${\pi }^{\mathsf{R}}\left(s\right|\mu ,\sigma ,\mathbf{y})$, respectively. On the other hand, $\mu \in \mathfrak{R},\sigma >0,s>1$, we considered the following change of variables $(\mu ,\sigma ,s)$ for $\mathit{\omega }=({\omega }_1,{\omega }_2,{\omega }_3)=(\mu ,log\left(\sigma \right),log(s-1))$. This modification leads the parametric space to ${\mathfrak{R}}^3$, which allows us to sample from the posterior distribution in a more efficient way. Considering the Jacobian transformation, the posterior distribution is given by:

$$\begin{array}{cc}\hfill {\pi }^{\mathsf{R}}\left(\mathit{\omega }\right|\mathbf{y})& \propto exp{\left({\omega }_2\right)}^{-n}exp\left({\omega }_3\right){\left[\left(1+\frac{1}{(1+exp\left({\omega }_3\right))}\right){\Psi }^{{}^{\prime }}\left(1+\frac{1}{(1+exp\left({\omega }_3\right))}\right)\right]}^{\frac{1}{2}}\hfill \\ \hfill & \times {(1+exp\left({\omega }_3\right))}^{-\frac{3}{2}}exp\left\{-\sum _{i=1}^n{\left(\frac{\Gamma \left(1+\frac{1}{(1+exp\left({\omega }_3\right))}\right)|y_i-{\omega }_1|}{exp\left({\omega }_2\right)}\right)}^{(1+exp\left({\omega }_3\right))}\right\}.\hfill \end{array}$$

The construction of the Metropolis–Hastings algorithm is done using a random walk for the parameter ${\omega }_1$, for instance the transition is obtained using $q({\omega }_1,{\omega }_1^{*})$, where ${\omega }_1^{*}$ is generated from ${\omega }_1^{*}={\omega }_1+kz$ where $z\sim N(0,{\eta }^2)$ and k is a constant that controls the acceptance rate. We have that ${\eta }^2$ is the diagonal element of the covariance matrix from the joint posterior distribution of $\mathit{\omega }$, obtained assuming the maximum a posterior estimate of ${\pi }^{\mathsf{R}}\left(\mathit{\omega }\right|\mathbf{y})$.

The computational stability was improved considering the logarithm scale. The steps to sample from the posterior distribution are:

• Set the values ${\mathit{\omega }}^{\left(0\right)}=\left({\omega }_1^{\left(0\right)},{\omega }_2^{\left(0\right)},{\omega }_3^{\left(0\right)}\right)=\left({\mu }^{\left(0\right)},log\left({\sigma }^{\left(0\right)}\right),log(s^{\left(0\right)}-1)\right)$.
• Generate ${\omega }_1^{*}$ from the proposal distribution $q({\omega }_1^{\left(0\right)},{\omega }_1^{*})$.
• Sample u from a uniform distribution $U(0,1)$.
• If $u\le min\left\{1,exp\left[log{\pi }^{\mathsf{R}}\left({\omega }_1^{*}|{\omega }_2^{\left(0\right)},{\omega }_3^{\left(0\right)},\mathbf{y}\right)-log{\pi }^{\mathsf{R}}\left({\omega }_1^{\left(0\right)}|{\omega }_2^{\left(0\right)},{\omega }_3^{\left(0\right)},\mathbf{y}\right)\right]\right\},$ then update ${\omega }_1^{\left(1\right)}$ from ${\omega }_1^{*}$; otherwise, use the value of ${\omega }_1^{\left(0\right)}$, i.e., ${\omega }_1^{\left(1\right)}={\omega }_1^{\left(0\right)}$.
• Repeat the same steps above for ${\omega }_2^{\left(1\right)}$ and ${\omega }_3^{\left(1\right)}$.
• Repeat Steps 2–5 until we obtain the target sample size.

After we generated the values of $\mathit{\omega }$, we computed $\left(\mu ,\sigma =exp\left({\omega }_2\right),s=1+exp\left({\omega }_3\right)\right)$. It was assumed that $\eta$ has the same values for all steps. The value of k is defined as aiming for the acceptance rate to be between $20\%$ and $80\%$ [27]. To confirm the convergence of the chains, the Geweke diagnostic [28] was used, as well as graphical analysis.

6. Selection Criteria for Models

In the Bayesian context, there are a variety of criteria that can be adopted to select the best fit between a collection of models. This paper considered the following criteria: the Deviation Information Criterion (DIC) defined by Spiegelhalter et al. [29] and the Expected Bayesian Information Criterion (EBIC) proposed by Brooks [30]. These criteria are based on the posterior mean deviation, $E\left(D\right(\omega \left)\right)$, which is the deviation evaluated at the posterior mean of $\omega$; thus, it is estimated with $\widehat{D}=\frac{1}{Q}{\sum }_{q=1}^{Q}D\left({\omega }^{\left(q\right)}\right)$, where the index q indicates the $q-th$ realization of a total of Q realizations and $D\left(\omega \right)=-2{\sum }_{i=1}^{n}log\left(f\left(y_i\right|\omega )\right)$, where $f(.)$ is the probability density of the GN distribution. The criteria DIC and EBIC are given by:

$$DIC=\overline{D}+p_D=2\overline{D}-\widehat{D}andEBIC=\overline{D}+2b,$$

where b is the number of parameters in the model and $p_D$ is the effective number of parameters, defined as $E\left[D\right(\omega \left)\right]-D\left[E\right(\omega \left)\right]$, where $D\left[E\right(\omega \left)\right]$ is the mean posterior deviation, which can be estimated as $\widehat{D}=D\left(\frac{1}{Q}{\sum }_{q=1}^{Q}D\left({\omega }^{\left(q\right)}\right)\right)$. Smaller values for the DIC and EBIC indicate the preferred model.

Another widely accepted criterion for model selection is the Conditional Predictive Ordinate (CPO). A detailed description of this selection criterion and the CPO statistic, as well as the applicability of the method for selecting models can be found in Gelfand et al. [31,32]. Let $\mathcal{D}$ denote the complete data set and ${\mathcal{D}}_{(-i)}$ denote the data with the i-th observation excluded. Consider $\pi \left(\omega \right|{\mathcal{D}}_{(-i)})$ for $i=1,\dots ,n$, the posterior density of $\omega$ given ${\mathcal{D}}_{(-i)}$. Thus, we can define the i-th observation of the CPO${}_i$ by:

$$\mathrm{C}\mathrm{P}{\mathrm{O}}_i={\int }_{\omega }f\left(y_i\right|\omega )\pi \left(\omega \right|{\mathcal{D}}_{(-i)})d\omega ={\left\{{\int }_{\omega }\frac{\pi \left(\omega \right|{\mathcal{D}}_{(-i)})}{f\left(y_i\right|\omega )}d\omega \right\}}^{-1},i=1,\dots ,n,$$

where $f\left(y_i\right|\omega )$ is the probability density function. High values of CPO${}_i$ indicate the best model. An estimate for CPO${}_i$ can be obtained using an MCMC sample of the posterior distribution of $\omega$ given $\mathcal{D}$, $\pi \left(\omega \right|\mathcal{D})$. For this, let ${\omega }_1,\dots ,{\omega }_Q$ be a sample of the distribution $\pi \left(\omega \right|\mathcal{D})$, of size Q. A Monte Carlo approximation of the CPO${}_i$ [23] is given by:

$$\widehat{\mathrm{C}\mathrm{P}{\mathrm{O}}_i}={\left\{\frac{1}{Q}\sum _{q=1}^Q\frac{1}{f\left(y_i\right|{\theta }_q)}\right\}}^{-1}.$$

A summary statistic of the CPO${}_i$’s is ${\sum }_{i=1}^{n}log\left({\widehat{\mathrm{CPO}}}_i\right)$, wherein the higher the value of B, the better the fit of the model. To illustrate the proposed methodology, a comparison between the normal and GN models is presented in Section 8.

7. Bayesian Case Influence Diagnostics

One way to observe the influence of observations on a fit of the model is via the global diagnosis; for instance, one can remove some cases from the analysis and analyze the effect of the removal [33]. The diagnostics of the case influence in a Bayesian perspective are based on the Kullback–Leibler divergence (K-L). Let $K(\pi ,{\pi }_{(-i)})$ denote the K-L divergence between $\pi$ and ${\pi }_{(-i)}$, where $\pi$ denotes the posterior distribution of $\omega$ for all data and ${\pi }_{(-i)}$ denotes the posterior distribution of $\omega$ without the i-th case. Specifically,

$$K(\pi ,{\pi }_{(-i)})=\int \pi \left(\omega \right|D)log\left\{\frac{\pi \left(\omega \right|D)}{\pi \left(\omega \right|D_{(-i)})}\right\}d\omega ,$$

(22)

and therefore, $K(\pi ,{\pi }_{(-i)})$ measures the effect of deleting the i-th case from the full one on the joint posterior distribution of $\omega$. The calibration can be obtained by solving the equation:

$$K(\pi ,{\pi }_{(-i)})=K(B\left(0.5\right),B\left(p_i\right))=\frac{-log\left\{4p_i(1-p_i)\right\}}{2},$$

where the Bernoulli distribution $B\left(p\right)$ is expressed by the parameter p with success probability[34] and $p_i$ is a calibration measure of the K-L divergence. After some algebraic manipulation, the obtained expression is:

$$p_i=0.5\left(1+\sqrt{1-exp\{-2K(\pi ,{\pi }_{(-i)})\}}\right),$$

bounded by $0.5\le p_i\le 1$. Therefore, if $p_i$ is significantly higher than 0.5, then the i-th case is influential.

The posterior expectation (22) can also be written in the form:

$$\begin{array}{cc}\hfill K(\pi ,{\pi }_{(-i)})& =log{E}_{\omega |D}\left\{{\left[f\left(y_i\right|\omega )\right]}^{-1}\right\}+E_{\omega |D}\{log\left[f\left(y_i\right|\omega )\right]\}\hfill \\ \hfill & =-log\left(\mathrm{C}\mathrm{P}{\mathrm{O}}_i\right)+E_{\omega |D}\{log\left[f\left(y_i\right|\omega )\right]\},\hfill \end{array}$$

(23)

where $E_{\omega |D}(.)$ denotes the expectation from the posterior $\pi \left(\omega \right|D)$. Thus, (23) can be estimated by using the MCMC methods to achieve the sample from the posterior $\pi \left(\omega \right|D)$. Therefore, if ${\omega }_1,\dots ,{\omega }_Q$ is a sample of size Q of $\pi \left(\omega \right|D)$, then:

$$\widehat{K(\pi ,{\pi }_{(-i)})}=-log\left(\widehat{\mathrm{C}\mathrm{P}{\mathrm{O}}_i}\right)+\frac{1}{Q}\sum _{q=1}^{Q}log\left[f\left(y_i\right|{\omega }^{\left(q\right)})\right].$$

8. Applications

In this section, the proposed method is illustrated using artificial and real data.

8.1. Artificial Data

An artificial sample of size $n=500$ was generated in accordance with (3) with $p=2,$${{\mathit{x}}_i}^{\top }=(1,x_{i1}),x_{i1}\sim N(2.5,1),\mathit{\beta }={(2,-1.5)}^{\top },\sigma =1$, and $s=2.5$. The posterior samples were generated by the Metropolis–Hastings technique through the MCMC implemented in the R software [35]. A single chain of dimensions 300,000 was considered for each parameter, and also, we discarded the first 150,000 iterations (burn-in), aiming to reduce correlation problems. A space with a size of 15 was used, resulting in a final sample size of 10,000. The convergence of the chain was verified by the criterion proposed by Geweke [28].

Table 1. Artificial data. Posterior mean, median, standard deviation (SD), and $95\%$ HDP intervals for the parameters of the model.

Parameters	Mean	Median	SD	95% HDP
${\mathit{\beta }}_1$	$1.995$	$1.996$	$0.086$	$(1.826;2.164)$
${\mathit{\beta }}_2$	$-1.510$	$-1.510$	$0.032$	$(-1.572;-1.447)$
$\sigma$	$1.027$	$1.028$	$0.055$	$(0.922;1.135)$
s	$2.657$	$2.631$	$0.320$	$(2.042;3.293)$

shows the posterior summaries for the parameters of the GN linear regression model. It can be seen that the estimates were close to the true values, and the $95\%$ HPD credibility intervals covered the true values of the parameters.

We used the same sample previously simulated to investigate the K-L divergence measure in detecting the GN linear regression model’s influential observations. We selected Cases 50 and 250 for perturbation. For each of these two cases and also considering both cases simultaneously, the response variable was disturbed as follows: ${\tilde{y}}_i=y_i+5S_y$, where $S_y$ is the standard deviation of $y_i$. The MCMC estimates were done similarly to those in the last section. Note that, due to the invariance property, $\mu$ can be computed for the standard GN distribution using the Bayes estimates of ${\mathit{\beta }}_1$ and ${\mathit{\beta }}_2$, that is $\mu ={\mathit{x}}_i^{\top }\mathit{\beta }$.

To reveal the impact of the influential observations in the estimates of ${\mathit{\beta }}_1,{\mathit{\beta }}_2,\sigma$, and s, we calculated the measure of Relative Variation (RV), which was obtained by $RV=\left|{\displaystyle \frac{\widehat{\theta }-{\widehat{\theta }}_0}{\widehat{\theta }}}\right|\times 100\%$, where $\widehat{\theta }$ and ${\widehat{\theta }}_0$ are the posterior averages of the model parameters considering the original data and perturbed data, respectively.

Table 2. Artificial data. Posterior mean, $RV(\%)$, and $95\%$ HDP intervals for the parameters of the model.

Data Names	Perturbed Case	Parameters	Mean	RV (%)	95% HDP
a	None	${\beta }_1$	$1.995$	−	$(1.826;2.164)$
		${\beta }_2$	$-1.510$	−	$(-1.572;-1.447)$
		$\sigma$	$1.027$	−	$(0.922;1.135)$
		s	$2.657$	−	$(2.041;3.293)$
b	50	${\beta }_1$	$2.026$	$1.554$	$(1.841;2.211)$
		${\beta }_2$	$-1.517$	$0.464$	$(-1.585;-1.449)$
		$\sigma$	$0.970$	$5.550$	$(0.866;1.073)$
		s	$2.188$	$17.651$	$(1.763;2.613)$
c	250	${\beta }_1$	$1.966$	$1.454$	$(1.775;2.154)$
		${\beta }_2$	$-1.495$	$0.993$	$(-1.564;-1.423)$
		$\sigma$	$0.916$	$10.808$	$(0.816;1.016)$
		s	$1.867$	$29.733$	$(1.565;2.165)$
d	$\{50,250\}$	${\beta }_1$	$2.002$	$1.332$	$(1.807;2.198)$
		${\beta }_2$	$-1.506$	$0.265$	$(-1.578;-1.436)$
		$\sigma$	$0.902$	$12.171$	$(0.801;1.003)$
		s	$1.773$	$33.271$	$(1.496;2.052)$

shows the posterior estimates for the artificial data and the RVs of the estimates of the parameters concerning the original simulated data. The data set denoted by (a) consisted of the original simulated data set without perturbation, and the data sets denoted by (b) to (d) consisted of data sets resulting from perturbations in the original simulated data set. Higher values of the relative variations for the parameters $\sigma$ and s showed the presence of influential points in the data set. However, the estimate of s did not differ so much from the perturbed Cases (c) and (d).

Considering the samples generated from the posterior distribution of the GN linear regression model parameters, we estimated the measures of the K-L divergence and their respective calibrations for each of the cases considered (a–d), as described in Section 7. The results in

Table 3. Diagnostic measures for the artificial data.

Data Names	Case Number	$\mathit{K}(\mathit{\pi },{\mathit{\pi }}_{(-\mathit{i})})$	Calibration
a	50	0.0121	0.5774
	250	0.0014	0.5262
b	50	16.1593	1.0000
c	250	19.2236	1.0000
d	50	2.8796	0.9992
	250	18.2292	1.0000

show that for the data without perturbation (a), the selected cases were not influential because they had small values for $K(\pi ,{\pi }_{(-i)})$ and calibration close to $0.577$. However, when the data were perturbed (b,d), the values of $K(\pi ,{\pi }_{(-i)})$ were more extensive, and their calibrations were close or equal to one, indicating that these data were influential.

8.2. Real Data Set

In order to illustrate the proposed methodology, recall the Brazilian Eucalyptus clones data set on the height (in meters) and the diameter (in centimeters) of Eucalyptus clones. The data belong to a large pulp and paper company from Brazil. As a strategy for the rising rentability of the forestry enterprise and keeping pulp and paper production under control, the company needs to keep an intensive Eucalyptus clone silviculture. The height of the trees is an important measure for selecting different clone species. Moreover, it is desirable to have trees of similar heights, possibly with a slight variation, and consequently with a distribution function with lighter tails.

The objective is to relate the tree’s diameter (explanatory variable) with its height (response variable). The GN and normal linear regression models were fit to the data via the Bayesian reference process. The posterior samples were generated by the Metropolis–Hastings technique, similar to the simulation study, in which we considered a single chain of dimensions 300,000 for each parameter, with a burn-in of 150,000 iterations; additionally, jumps with a size of 15 were used, resulting in a final sample size of 10,000. The Geweke criteria verified the convergence of the chain.

Table 4. Real data. Posterior mean and $95\%$ HPD intervals for the parameters of the model and Bayesian comparison criteria.

Model	Parameters	Mean	95% HDP	DIC	EBIC	B
	${\beta }_1$	7.066	(6.398; 7.724)
Generalized	${\beta }_2$	0.948	(0.907; 0.990)	5829.180	5847.96	−2914.595
Normal	$\sigma$	3.237	(3.034; 3.440)
	s	2.865	(2.436; 3.305)
	${\beta }_1$	7.002	(6.331; 7.671)
Normal	${\beta }_2$	0.949	(0.907; 0.991)	6248.652	6262.98	−3126.614
	$\sigma$	1.993	(1.916; 2.068)

shows the posterior summaries for the parameters of both distributions and the model selection criteria. The GN linear regression model was the most suitable to represent the data as it performed better than the normal linear regression model for all the criteria used. For the fit regression model chosen, note that ${\beta }_2$ was significant, and then, for every one-unit increase in the diameter, the average height of Eucalyptus 0.95 meters increased. In particular, the analysis of the shape parameter ($s>2$) provided strong evidence of a platykurtic distribution for the errors, and this favored the GN linear regression model. This was further confirmed by graphical analyses of the quantile residuals of the GN linear regression model presented in Figure 2d.

Figure 2a and Figure 3a show the scatter plot of the data and the adjusted normal and GN linear regression models. It was observed that, on average, the estimated heights were close to those observed, indicating that the models considered had a good fit. The residuals graph by the fit values and the residuals graph by the observations were also quite similar for both models. The presence of heteroscedasticity (see Figure 2b and Figure 3b) was noted, as well as the quadratic trend (see Figure 2c and Figure 3c) for the height-to-diameter ratio of the Eucalyptus. The graphs of the quantiles of the GN distribution and the normal distribution for the residuals of the models are presented in Figure 2d and Figure 3d, respectively. It can be seen that in the setting of the normal linear regression model, many points were far from their tails, indicating an inadequate specification of the error distribution for the model. On the other hand, using our proposed approach in the GN distribution, we observed a good fit as points were following the line, indicating that the theoretical residuals were close to the observed residuals. Therefore, there was evidence that the model chosen outperformed the normal linear regression model in fitting the data.

To investigate the influence of height and diameter data from Eucalyptus on the fit of the generalized normal linear regression model chosen, we calculated the K-L divergence measures and their respective calibrations. Figure 4 shows the K-L measurements for each observation. Note that Cases 335, 479, and 803 exhibited higher values of the K-L divergence when compared with other observations. The K-L divergences and calibrations concerning three observations that showed the highest calibration values are presented in

Table 5. Diagnostic measures for the height and diameter data of Eucalyptus.

Case Number	$\mathit{K}(\mathit{\pi },{\mathit{\pi }}_{(-\mathit{i})})$	Calibration
335	0.1683	0.7673
479	0.1734	0.7707
803	0.4936	0.8960

. It can be seen that Observation 803 was possibly an influential case, which was also shown to be an outlier from the visual analysis. To assess whether this observation altered the parameter estimates of the GN linear regression model, we carried out a sensitivity analysis.

Table 6. Posterior estimates and $RV(\%)$ for Eucalyptus height and diameter data following the removal of the influential case.

Cases Deletions	Parameter	Mean	$\mathit{RV}(\%)$	95% HDP
803	${\beta }_1$	7.067	0.014	(6.375; 7.753)
	${\beta }_2$	0.948	-	(0.905; 0.991)
	$\sigma$	3.252	0.463	(3.054; 3.454)
	s	2.936	2.478	(2.486; 3.390)

shows the new estimates of the model parameters after excluding the case with the greatest calibration value and the relative variations $\left(RV\right)$ for these estimates regarding the Eucalyptus data. Here, the relative variations were obtained by $RV=\left|{\displaystyle \frac{\widehat{\theta }-{\widehat{\theta }}_0}{\widehat{\theta }}}\right|\times 100\%$, where $\widehat{\theta }$ and ${\widehat{\theta }}_0$ are the posterior averages of the model parameters obtained from the original data and from the data without influential observation, respectively. We noted a slight change in the RV of the s parameter when we excluded influential observations. However, such a change was insignificant. This indicated that the GN linear regression model was not affected by the compelling cases.

Overall, regardless of the unit measurement/scale in both cases (synthetic and real data sets), the visual representation corroborated the explainability/adjustment of the GN model, by the quantile-quantile plot, residuals versus fits plot, and standardized errors showing no pattern left to explain equal error variances, as well as no outliers.

9. Discussion

In this paper, we presented the generalized normal linear regression model from objective Bayesian analysis. The Jeffreys prior and reference prior for the generalized normal model were discussed in detail. We proved that the Jeffreys prior leads to an improper posterior distribution and cannot be used in a Bayesian analysis. On the other hand, the reference prior leads to a proper posterior distribution.

The parameter estimates were based on a Bayesian reference analysis procedure via MCMC. Diagnostic techniques based on the Kullback–Leibler divergence were built for the generalized typical linear regression model. Studies with artificial and real data were performed to verify the adequacy of the proposed inferential method.

The result of the application to a set of actual data showed that the generalized normal linear regression model outperformed the normal linear regression model, regardless of the model selection criteria. Furthermore, through a study of artificial data and real data, the Kullback–Leibler divergence effectively detected the points that were influential in the fit of the generalized normal linear regression model. The withdrawal of such influential points from the set of real data showed that the generalized normal model was not affected by influential observations. This result was corroborated by the fact that the generalized normal distribution was considered a tool for reducing outliers and achieving robust estimates. The proposed methodology showed consistent marginalization and sampling properties and thus eliminated the problem of estimating the parameters of this important regression model. Moreover, adopting the reference (objective) prior, we obtained one-to-one consistent results, under the Bayesian paradigm, enabling the GN distribution in a practicalform.

Further works can explore a great number of extensions using this study. For instance, the method developed in this article may be applied to other regression models such as the Student t regression model and the Birnbaum–Saunders regression model, among others. Additionally, other generalizations of the normal distribution should be considered [36,37,38].