title: Improved estimators in beta prime regression models
authors: Francisco M. C. Medeiros Mariana C. Araújo Marcelo Bourguignon
text: 
  
  
  Concentration for random product formulas
      Joel A. Tropp
      December 30, 2023
  =========================================
  
  
  
  
  In this paper, we consider the beta prime regression model recently proposed by <cit.>, which is tailored
  to situations where the response is continuous and restricted to the positive real line with skewed and long tails
  and the regression structure involves regressors and unknown parameters.
  We consider two different strategies of bias correction of the maximum-likelihood estimators
  for the parameters that index the model. In particular, we discuss bias-corrected estimators for the mean and the dispersion
  parameters of the model. Furthermore, as an alternative to the two analytically bias-corrected estimators discussed,
  we consider a bias correction mechanism based on the parametric bootstrap.
  The numerical results show that the bias correction
  scheme yields nearly unbiased estimates.
  An example with real data is presented and discussed.
  
  
  
  
    
  Keywords: Beta prime distribution; Bias correction; Bootstrap; Dispersion covariates; Maximum-likelihood.
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  The beta prime (BP) distribution  (known as inverted beta distribution or beta distribution of the second kind as well) is a two-parameter distribution on the positive real line, which can be interpreted as the distribution of the odds ratio of a variable distributed according to the beta distribution, i.e., if X has a beta distribution with parameters α and β, then Y = X/(1-X) has a BP distribution with α > 0 and β > 0 both shape parameters. We are adopting the parameterization for the BP distribution in terms of the mean and precision parameters which was proposed by <cit.>. An advantage of using this parameterization is that we can introduce regression structures for each mean and precision parameters and the interpretation of the regression coefficients is straightforward in terms of them as in generalized linear models. Thus, the BP random variable Y (<cit.>) is defined as follows: Let Y be a random variable with probability density function (pdf) given by
  
      f(y; μ, ϕ) = y^μ(ϕ + 1)-1(1 + y)^-[μ(ϕ + 1) + ϕ + 2]/B(μ(1+ ϕ), ϕ + 2),    y > 0,
  
  where μ > 0 and ϕ > 0 are mean and precision parameters, respectively, B(μ(1+ ϕ), ϕ + 2) = Γ(μ(1+ ϕ))Γ(ϕ + 2)/Γ(μ(1+ ϕ) + ϕ + 2) is the beta function and Γ(μ(1+ ϕ)) = ∫_0^∞ω^μ(1+ ϕ)-1e^-ωdω is the gamma function. From now on, we use the notation Y ∼BP(μ, ϕ) to indicate that Y is a random variable following a BP distribution. The mean and variance of Y are
  
      E[Y] = μ  and  Var[Y] = μ(1+μ)/ϕ.
  
  
  Some features of the BP model are <cit.>: first, the variance function of the BP model assumes a quadratic form similar to the gamma distribution. However, the variance function of the proposed model is larger than the variance function of gamma distribution, which may be more appropriate in certain practical situations; second, the BP hazard rate function can have an upside-down bathtub or increasing depending on the parameter values. The most classical two-parameter distributions such as Weibull and gamma distributions have monotone hazard rate functions; third, the skewness and kurtosis of the BP distribution can be much larger than those of the gamma and inverse gaussian distributions; fourth, there are some stochastic representation of the BP random variable.
  
  
  
  
  
  In the literature there are only a few works dealing with the BP distribution.
  <cit.> discussed its properties and obtained the maximum likelihood (ML) estimates of the model parameters.
  Bias-corrected versions of the MLEs of the parameters that index the BP distribution were obtained by <cit.>. It is worth mention that all the works related above have considered the usual parameterization of the BP distribution. Considering the parameterization we adopted, <cit.> used the the ML method for estimating the parameters that index the BP regression model. However, as can be seen in Table 2 in <cit.>,  in small-sized samples, the ML estimators of these parameters (especially for precision structure) may be extremely biased.
  So, it is important consider alternative estimators with smaller biases when the number of observations are small.
  
  Investigates how the maximum likelihood estimator behaves in small-sized sample, in particular bias analysis, is an important research area. In regular parametric statistical models the maximum likelihood estimator bias is generally of the order 𝒪(n^-1) for large sample size n and are, in practice, usually ignored since that the asymptotic standard error is of order O(n^-1/2). When dealing with small-sized sample, however, bias can be a problematic issue, thus it can not be neglected. So, it is important to obtain bias correction in these cases. Bias reduction was studied by several authors. In uniparametric models, <cit.> obtained an expression for the 𝒪(n^-1) bias from the maximum likelihood estimator. Assuming independent, but not necessarily identically distributed observations, <cit.> obtained a general expression for the 𝒪(n^-1) bias of the maximum likelihood estimator in multiparametric models. This result has become widely used in the literature to obtain general expressions for the 𝒪(n^-1) bias and to propose bias-corrected estimators in various parametric models. For instance, <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. 
  
  Usually, the approach to obtain bias-corrected versions of the MLEs uses the second order bias. In this procedure the adjustment is made after the MLEs were computed. Additionally, an alternative approach was proposed by <cit.>, who suggested that a bias reduction method by modifying the score function previous to obtain the parameter estimates. This method is called the preventive method and has been studied in parametric models where maximum likelihood estimates can be unstable (infinite or belonging to the parametric space boundary) such as <cit.>, <cit.>, <cit.>, <cit.> and <cit.>. Another possible way to perform bias correction is through bootstrap resampling, which requires no explicit derivation of the bias function.
  In this context, the main goal of this paper is to derive a closed-form expression for the second order biases of the ML estimators in the BP regression model which can be used to define bias corrected ML estimators to order 𝒪(n^-1).
  
  This paper is organized as follows: this introductory
  section. In Section <ref>, the BP regression model is introduced and some of
  its basic properties are outlined.
  In Section <ref>, we obtain the second order biases of the MLEs of the means of the responses and precision parameters of the model.
  Section <ref> discusses the numerical results.
  In Section <ref>, we consider an empirical example. Finally, Section <ref> concludes the paper.
  
  
  
  § BETA PRIME REGRESSION MODEL
  
  
  
  Consider n independent random variables Y_1, …, Y_n where each Y_i, i = 1, …, n has BP distribution with pdf given by (<ref>) with mean μ_i and precision parameter ϕ_i. <cit.> proposed the BP regression model which is defined by (<ref>) and by two functional relations
  
      g_1(μ_i) = η_1i = 𝐱^⊤_iβ  and   g_2(ϕ_i) = η_2i = 𝐳^⊤_iν,
  
  where g_1: ℝ→ℝ^+ and g_2: ℝ→ℝ^+ are strictly monotone, positive and at least twice differentiable link functions, η_1i and η_2i are the linear predictors, β = (β_1, …, β_p)^⊤ (β∈ℝ^p) and ν = (ν_1, …, ν_q)^⊤ (ν∈ℝ^q,   q< n-p) are unknown parameter vectors to be estimated, and 𝐱_i = (x_i1, …, x_ip)^⊤ and 𝐳_i = (z_i1, …, z_iq)^⊤ are observations on p and q known regressors, for i = 1, …, n. Additionally, we assume that the covariate matrices 𝐗 = (𝐱_1, …, 𝐱_n)^⊤ and 𝐙 = (𝐳_1, …, 𝐳_n)^⊤ have rank p and q, respectively. Besides the interpretation of the regression coefficients being in terms of the mean and precision parameters, another advantage of the model proposed by (<ref>) and (<ref>) is that it is suitable for modeling asymmetric data, being an alternative to the generalized linear models when dealing with asymmetric dataset.
  
  
  The log-likelihood function for (β, ν) given the observed values y_1, …, y_n is
  
      ℓ(β, ν) = ∑_i=1^nℓ(μ_i, ϕ_i),
  
  being
  
      ℓ(μ_i, ϕ_i)   =    [μ_i(1+ϕ_i) - 1]log(y_i) - [μ_i(1+ ϕ_i) + ϕ_i +2]log(1+y_i) 
         -   log[ Γ(μ_i(1+ϕ_i))]- log[ Γ(ϕ_i+2)] + log[ Γ(μ_i(1+ϕ_i) + ϕ_i +2)];
  
  μ_i=g_1^-1(η_1i) and ϕ_i=g_2^-1(η_2i) are functions of β and ν, respectively, as defined in (<ref>). A method for obtaining the parameters estimates of the BP model defined by (<ref>) and (<ref>) is described in details in <cit.>. They consider the gamlss function for this purpose.
  
  We assume that the log-likelihood function (<ref>) satisfies the usual regularity conditions of large sample likelihood theory <cit.>. Thus, when n is large and under some regular conditions we have that
  
      ([ β; ν ])
      ∼a_p + q(
      ([ β; ν ])
      , 𝐊(β,ν)^-1),
  
  
  where ∼a means “approximately distributed” and 𝐊(β,ν)^-1 is the inverse of Fisher's information matrix evaluated at β and ν, which can be approximated by J(β,ν)^-1, where -J denotes the (p+q) × (p+q) Hessian matrix evaluated at (β^⊤,ν^⊤)^⊤. Fisher's information matrix 𝐊(β,ν) is presented in Appendix A.
  
  
  
  
  
  
  
  
  
  
    
  
  
  
  
  
  
  
  
  
  
  
  
  
    
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
       
       
  
  
  
  
  
  
  
       
       
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § BIAS CORRECTION OF THE MLES
  
  
  
  Let θ=(β^⊤, ν^⊤)^⊤ be the unknown parameter vector of the BP regression model. We now obtain an expression for the second order biases of the MLEs of the components of θ using Cox and Snell's (<cit.>) general formula. In order to obtain this expression, we first introduce some notation. The lower subscripts r, s, t, u, … and the upper subscripts R, S, T, U, … denote, respectively, the components of β and ν vectors. Therefore, the partional derivatives of the log-likelihood (<ref>) with respect to the components of β and ν are presented as  U_r= ∂ℓ/∂β_r, U_rS= ∂^2 ℓ/ ∂β_r ∂ν_S, U_rST=∂^3 ℓ/ ∂β_r ∂ν_S ∂ν_T, etc. The moments of the log-likelihood derivatives are represented by κ_rs=E(U_rs), κ_r,s=E(U_r U_s), κ_r, ST=E(U_r U_ST), etc, where all κ's regard to a total covering the whole sample and are, in general, of order 𝒪(n^-1). The moments derivatives are defined by κ_rs^(t)= ∂κ_rs/∂β_t, κ_rs^(T)= ∂κ_rs/∂ν_T, etc. Finally, we denote the elements of the inverse of Fisher's information matrix K(β, ν)^-1=K(θ)^-1, which are 𝒪(n^-1), as κ^r,s = -κ^rs, κ^r,S = -κ^rS, κ^R,s = -κ^Rs and κ^R,S = -κ^RS.
  
  From the general Cox and Snell's () formula we can obtain the 𝒪(n^-1) bias of the MLE for the ath component of the parameter vector θ=(θ_1, …, θ_p,
  θ_p+1, …, θ_p+q)^⊤ = (β^⊤, ν^⊤)^⊤ as:
  
  
      B_θ(θ_a)    =   ∑_r,s,uκ^arκ^su{κ_rs^(u) - 1/2κ_rsu} + ∑_R, s, uκ^aRκ^su{κ_Rs^(u) - 1/2κ_Rsu}
         +   ∑_r, S, uκ^arκ^Su{κ_rS^(u) - 1/2κ_rSu} + ∑_r,s,Uκ^arκ^sU{κ_rs^(U) - 1/2κ_rsU}
         +   ∑_R, S, uκ^aRκ^Su{κ_RS^(u) - 1/2κ_RSu} + ∑_R, s, Uκ^aRκ^sU{κ_Rs^(U) - 1/2κ_RsU}
         +   ∑_r,S,Uκ^arκ^SU{κ_rS^(U) - 1/2κ_rSU} + ∑_R, S, Uκ^aRκ^SU{κ_RS^(U) - 1/2κ_RSU}.
   
  
  From (<ref>), we can observe that β and ν are not orthogonal, hence all terms in (<ref>) must be considered. In order to save space, all cumulants needed to obtain (<ref>) are given in Appendix B. After a long algebra presented in details in Appendix B we achieve to the expressions for the second order biases of β and ν, given in matrix form respectively by
  
      B_β(β)     =     K^ββ X^⊤[ M_1 P_ββ + (M_2 + M_3) P_βν + M_5 P_νν] 
          +     K^βν Z^⊤[ M_2 P_ββ + (M_4 + M_5)P_βν + M_6 P_νν]
  
  and
  
      B_ν(ν)     =     K^βν X^⊤[M_1P_ββ + (M_2 + M_3)P_βν + M_5 P_νν]
          +     K^νν Z^⊤[ M_2P_ββ + (M_4 + M_5)P_βν + M_6 P_νν],
  
  where K^ββ, K^βν and K^νν represent matrices which components are respectively the (r,s)th, (r,S)th and (R,S)th elements of the inverse of Fisher's information matrix, M_1 to M_6 are presented in Appendix B, P_ββ, P_βν and P_νν are vectors with the same n × 1 dimension and which elements are the diagonal elements of X K^ββ X^⊤, X K^βν Z^⊤ and Z K^νν Z^⊤, respectively.
  
  We now assume the 2n × 1 vector δ_1 defined as
  
      δ_1 = ([ M_1P_ββ + (M_2 + M_3)P_βν + M_5 P_νν; M_2P_ββ + (M_4 + M_5)P_βν + M_6 P_νν ]),
   and consider K^β *=(K^ββ  K^βν) and K^ν *=(K^νβ   K^νν), the p × (p+q) upper and q × (p+q) lower blocks of the matrix K(β, ν)^-1, respectively. Thus, we can express the second-order biases of β and ν as
  
      B_β(β) = K^β *𝕏^⊤δ_1       B_ν(ν) = K^ν *𝕏^⊤δ_1.
  
  From the expressions above, we can obtain in matrix form the second order bias of the MLE of the joint vector θ=(β^⊤, ν^⊤)^⊤ expressed as
  
      B_θ(θ)=(𝕏^⊤K𝕏)^-1𝕏^⊤δ_1.
  
  
  
  Now, we define the bias-corrected estimator as
  
      θ= θ - B_θ(θ),
  
  
  
  
  
  where B_θ(θ) is bias of the θ with the unknown parameters replaced by their MLEs. Considering the assumptions assumed in Section <ref>, we have that the asymptotic distribution of θ is _p + q(θ, 𝐉(θ)^-1 ), where 𝐉(θ) = 𝐉(β,ν)^-1
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  A second approach to correct the second order bias of the MLE of θ=(β^⊤, ν^⊤)^⊤ is considering the “preventive” method proposed by <cit.>. This method basically consists of modify the original score function in order to remove the 𝒪(n^-1) bias. The modified score function is given by
  
      U^*(θ) = U(θ) -
      K(θ)B_θ(θ),
  
  being K(θ) the information matrix and B_θ(θ) the 𝒪(n^-1) bias. Considering the BP regression model and replacing the expression obtained for B_θ(θ), the modified score function has the following form:
      U^*(θ) = U(θ) -  𝕏^⊤δ_1.
   The second order bias corrected MLEs θ̌ is the solution of U^*(θ)=0. Also, θ̌ is asymptotically normal distributed as _p + q(θ, 𝐉(θ)^-1 ), with 𝐉(θ) as given previously.
  
  Another way to bias-correcting the MLEs of the regression parameters is by the bootstrap technique <cit.>. In this paper, in order to reduce the computational burden, we shall adopt the warp-speed bootstrap method of <cit.> for evaluating the proposed resampling scheme. The warp-speed bootstrap method follows the steps described below. Instead of computing the MLEs for each Monte Carlo sample r=1,2,…,m (with m being the total number of Monte Carlo replications) on the basis of B bootstrap samples, just one resample (i.e. B=1) is generated from the assumed model with the parameters replaced by estimates of maximum likelihood computed using the original sample for each Monte Carlo sample and, hence, estimates of maximum likelihood, say θ^*, is computed for that sample. Therefore, the bootstrap bias estimates θ is
  
      B_θ(θ^*) = θ^* - θ.
  
  By using the bootstrap bias estimate presented above, we arrive at the following bias-corrected, to order O(n^-1), estimator:
  
      θ̃^b = 2θ - θ^*.
  
  
  For a good discussion to the bootstrap method, see <cit.>. Finally, it is worth mentioning that the idea behind the warp-speed bootstrap method is that taking just one bootstrap draw for each simulated sample is sufficient to provide a useful approximation to the bias of estimator. Applying this insight to Monte Carlo evaluation of bootstrap-based bias yields evaluation methods that work with M=1 <cit.>. Due to the resulting dramatic computational savings, <cit.> called their method as “Warp-Speed” Monte Carlo method. Therefore, the bootstrap-based bias on the basis of warp-speed bootstrap method become a viable alternative to inferential improvements in small samples when there are impeditive or too costly analytical difficulties.
  
  
  
  
  
  § NUMERICAL RESULTS
  
  
  
  We now present a Monte Carlo simulation study to investigate and compare the performance of the MLEs along with their corrected versions proposed in this article in small and moderate-sized samples. We use a BP regression models with dispersion covariates and a log link. We consider the model
  
      log(μ_i) = β_0 + ∑_ℓ=1^p β_ℓx_iℓ  and  log(ϕ_i) =ν_0 + ∑_ℓ=1^q ν_ℓx_iℓ,    i=1,2,…,n,
  
  where the true values of the parameters were taken as 1. The covariates values are taken as random draws from the U(0,1) distribution and their values were held constant throughout the simulations. We consider different values for the number of regression parameters (p and q) and the sample size (n=30, 40 and 60). The number of Monte Carlo replicates was 10.000 and all the simulations were performed using the R language <cit.>. In each Monte Carlo replica, we computed the MLEs of the parameters, their corrected versions from the corrective method <cit.>, preventive method <cit.>, and the parametric version of the bootstrap method <cit.>. In order to analyze the results, we computed, for each sample size and for each estimator, the mean of estimates, bias,
  variance and mean square error (MSE). The results are presented in Tables <ref>, <ref> and <ref> for p=q=1, p=q=2, and p=q=3, respectively.
  
  Tables <ref>-<ref> summarize the simulation results for the β's and ν's varying the sample size n and the number of regression parameters (p and q). As can be seen in Tables <ref>-<ref>, for most part of the parameters the estimated biases, in absolute value, of the original MLEs were larger than the others. In general, for the β's, the biases of preventive estimators were smaller than those of the corrective estimators and bootstrap estimators. For the ν's, the biases of the bootstrap estimators were, in general, smaller than those of the corrective estimators and the preventive estimators. These performances are independent of the number of parameters to be estimated. For instance, in absolute value, when p=q=1 and n=30 the bias of the parameter β_0 were 0.0048 (MLE), 0.0009 (Cox-Snell), 0.0000 (Firth) and 0.0013 (p-boot) and the bias of the parameter ν_0 were 0.1044 (MLE), 0.0148 (Cox-Snell), 0.0083 (Firth) and 0.0064 (p-boot); see Table <ref>. However, for all parameters, in most cases the MSE of the corrective estimators were the smallest and the MSE of the bootstrap estimators were the largest, followed by the MSE of the preventive estimators. When we increase the sample size for the β's, the bootstrap estimators tends to shows the smallest bias, although they have the largest MSE. For instance, for p = q = 3 and n = 30  we have that the bias, in absolute value, for β_1 were 0.0026 (MLE), 0.0024 (Cox-Snell), 0.0016 (Firth) and 0.0031 (p-boot), while when n=60 they were 0.0018 (MLE), 0.0016 (Cox-Snell), 0.0015 (Firth) and 0.0014 (p-boot); see Table <ref>. Comparing the results presented in each Table, we can observe that as the sample size increases, in general, the bias of the estimators reduces, as expected.
  
  The previous findings are confirmed by the box plots shown in Fig. 1, which were obtained for sample size n = 30.  In summary, the bias of the MLEs, especially for the ν's parameters, are larger than the bias of the corrected estimators. Box plots for different values of n, p, and q (not shown) exhibited a similar pattern. Therefore, we recommend the use of method (Cox-Snell, Firth or parametric bootstrap) to reduce bias in small and moderate sample size.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § AN APPLICATION
  
  
  
  In order to illustrate the proposed methodology, in this section, we apply the estimation methods considered in the previous section to a real situation. We consider the real dataset used in <cit.>. The main purpose is to assess sediment quality using the freshwater clam Corbiculafluminea to determine its adequacy as a biomonitoring tool in relation to theoretical risk indexes and regulatory thresholds. The study contains 27 observations (small-sized sample), which measured, among other characteristics, the dry weight tissue of the clams (, in g), wet weight tissue (, in g), and the concentrations of caesium () in the soft tissue. Such minerals were considered in 100 micrograms per liter (100μ g /L).
  
  
  We adopted a BP model to fit the dry weight tissue of the clams; that is,
  we consider that _i ∼ BP(μ_i,ϕ_i) with systematic components given by
  
      log(μ_i)     = β_0 + β_1_i + β_2_i 
      log(ϕ_i)    = ν_0 + ν_1_i,    i=1, …, 27.
  
  
  An R implementation for obtaining MLEs along with their corrected versions proposed in this article related to the data used is available at GitHub BiasBPR[<https://github.com/sesiommedeiros/BiasBPR>] repository.
  
  Table <ref> presents the maximum likelihood estimates along with their corrected versions and the corresponding estimates of asymptotic standard errors in parentheses. Note that for the parameters that model the precision, ν_0 and ν_1, the maximum likelihood estimates are smaller than the bias corrected ones, while for the parameters that model the mean, β_0, β_1 and β_2, the estimates are quite close.
  
  
  
  
  
  
  
  
  In the Table <ref>, we present the relative changes (RCs). The RCs are
  calculated from RC(θ) = |(θ- θ_o)/θ_o| × 100%,
  where θ denotes the MLE of θ and θ_o denotes the bias-corrected MLE of θ.
  From Table <ref>, the bias-corrected MLEs for ν_0 and ν_1 present similar results.
  In contrast to the Cox-Snell and p-boot bias-corrected estimators, the preventive method (Firth) gives estimates that dramatically change for β_1 and β_2. For example, the second-order bias is 36.585% of the total amount of the MLE of β_2. Thus, this real example illustrates that bias corrections can have a great effect on the conclusions.
  
  
  
  
  
  
  
  § CONCLUDING REMARKS
  
  
  
  In this paper, we have examined a wide range of estimators for the unknown parameter vector of the BP regression model.
  In particular,  we have derived a closed-form expression, in matrix form, for the second order biases of the ML estimators of the parameters that index the BP regression model proposed by <cit.>. For this, we use the expressions obtained through Cox and Snell's <cit.> formulae and Firth's <cit.> estimating equation. We also considered a bias correction based on parametric bootstrap. The numerical evidence here presented shows that our proposed estimators has good finite-sample behavior, even when the sample size is small. For the mean structure, we observe that the MLE presents a very small bias  (even in small samples). In this case, it is not necessary to use the bias-corrected estimators. However, for the precision structure, we observe
  that the MLE can become considerably biased and, therefore, we strongly recommend its bias correction. This behavior was also observed in the application to the real data set presented, therefore, we strongly recommend that practitioners use these corrected estimators when modeling data using the BP regression model. Finally, we have applied our proposed estimators to a real data.
  
  
  
  
  
  
  
  
  
  
  9
  
  [Barreto-Souza and Vasconcellos(2011)]SouzaVasconcellos2011
  Barreto-Souza W., Vasconcellos K. L. (2011). Bias and skewness in a general
  extreme-value regression model. Computational Statistics & Data Analysis, 55(3), 1379–1393.
  
  [Bartlett(1953)]Bartlett1953
  Bartlett M. S. (1953). Approximate Confidence Intervals. Biometrika, 40(1-3), 12–19.
  
  [Bonnail et al.(2016)]Bonnail
  Bonnail E., Sarmiento A. M., DelValls, T. A. and Nieto, J.M. and Riba, I. (2016). Assessment of metal contamination, bioavailability, toxicity and bioaccumulation in extreme metallic environments (Iberian Pyrite Belt) using Corbicula fluminea. Science of the Total Environment, 544, 1031–1044.
  
  [Bourguignon et al.(2018)]bour18
  Bourguignon M., Santos-Neto M., de Castro, M. A new regression model for positive data. <https://arxiv.org/abs/1804.07734v1>, 156–193.
  
  [Bull et al.(2002)]Bull2002
  Bull S. B., Mak C., Greenwood C. M. (2002). A modified score function estimator for multinomial logistic regression in small samples. Computational Statistics & Data Analysis, 57–74.
  
  [Cox and Hinkley(1983)]cookwei:83
  Cox D., Hinkley D. (1983). Theoretical Statistic. Chapman and Hall. 
  
  [Cox and Snell(1968)]coxsnell68
  Cox D. R., Snell E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society Series B (Methodological), 30(2), 248–275.
  
  [Cysneiros et al.(2010)]Cysneiros-et-al-2010
  Cysneiros F. J., Cordeiro G. M., Cysneiros A. H. (2010). Corrected maximum likelihood estimators in heteroscedastic symmetric nonlinear models. Journal of Statistical Computation and Simulation, 80(4), 451–461.
  
  [Efron and Tibshirani(1993)]EfronTibshirani1993
  Efron B., Tibshirani R. (1993). An Introduction to the Bootstrap. Chapman and Hall. New York.
  
  [Firth(1993)]Firth1993 
  Firth D. (1993). Bias reduction of maximum likelihood estimates. Biometrika, 80(1), 27–38.
  
  [Giacomini et al.(2013)]Giacomini-et-al-2013 
  Giacomini R., Politis D. N., White H. (2013). A Warp-Speed method for conducting Monte Carlo experiments involving bootstrap estimators. Econometric Theory, 29(3), 567–589.
  
  [Kosmidis(2014)]Kosmidis2014 
  Kosmidis I. (2014). Improved estimation in cumulative link models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1), 169–196.
  
  [Kosmidis and Firth(2009)]Kosmidis2009 
  Kosmidis I., Firth D. (2009). Bias reduction in exponential family nonlinear models. Biometrika, 96(4), 793–804.
  
  [Kosmidis and Firth(2011)]Kosmidis2011 
  Kosmidis I., Firth D. (2011). Multinomial logit bias reduction via the Poisson log-linear model. Biometrika, 98(3), 755–759.
  
  [Lemonte et al.(2007)]Lemonte-et-al-2007 
  Lemonte A. J., Cribari-Neto F., Vasconcellos K. L. (2007). Improved statistical inference for the two-parameter Birnbaum-Saunders distribution. Computational Statistics & Data Analysis, 51(9), 4656–4681.
  
  [McDonald(1987)]McDonald:1987aa 
  McDonald J. (1987). Model selection: some generalized distributions. Communications in Statistics - Theory and Methods, 16, 1049–1074.
  
  [Melo et al.(2018)]melo2018 
  Melo T. F. N., Ferrari S. L. P. , Patriota A. G. (2018). Improved estimation in a general multivariate elliptical model. Brazilian Journal of Probability and Statistics, 32(1), 44–68.
  
  
  [R Core Team(2017)]r:2017
  R Core Team (2017). 
  R: A Language and Environment for Statistical Computing. 
  R Foundation for Statistical Computing, Vienna, Austria.
  <https://www.R-project.org/>
  
  
  [Sartori(2006)]Sartori2006 
  Sartori N. (2006). Bias prevention of maximum likelihood estimates for scalar skew normal and skew t distributions. Journal of Statistical Planning and Inference, 136(12), 4259–4275.
  
  [Simas et al.(2011)]Simas2011 
  Simas A. B., Rocha A. V., Barreto-Souza W. (2011). Bias-corrected estimators for dispersion models with dispersion covariates. Journal of Statistical Planning and Inference, 141(9), 3063–3074.
  
  [Stosić and Cordeiro(2009)]borkocordeiro 
  Stosić B. D., Cordeiro G. M. (2009). Using Maple and Mathematica to derive bias corrections for two parameter distributions. Journal of Statistical Computation and Simulation, 79(6), 751–767.
  
  
  
  
  
  
  
  
  
  
  
  § APPENDIX A
  
  
  
  In this Appendix, we presente the Fisher's information matrix for the BP regression model, which is expressed in matrix form as
  
      𝐊(β, ν)=(
      [ 𝐗^⊤ K_ββ𝐗 𝐗^⊤ K_βν𝐙; 𝐙^⊤ K_νβ𝐗 𝐙^⊤ K_νν𝐙;           ]),
  
  
  where
  
      K_ββ   =   { (1 + ϕ_i)^2 a_i (μ_i/η_1i)^2}, 
      K_νν = { b_i (ϕ_i/η_2i)^2},
      
      K_βν   =    K_νβ^⊤ = { (1 + ϕ_i) [a_i μ_i - ψ^(1)(μ_i(1 + ϕ_i) + ϕ_i + 2)] μ_i/η_1iϕ_i/η_2i} ,
  
  with a_i = ψ^(1)(μ_i(1+ϕ_i)) - ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) and b_i = μ_i^2 ψ^(1)(μ_i(1+ϕ_i)) - (1 + μ_i)^2 ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) + ψ^(1)(ϕ_1 + 2).
  
  
  We can rewrite the Fisher's information matrix given in (<ref>). For this, let K be a 2n × 2n matrix and 𝕏 be a 2n × (p+q) matrix defined, respectively, as
  
      K   =   (
      [ K_ββ K_βν; K_νβ K_νν ]).
  
  and
  
      𝕏   =   (
      [ X 0; 0 Z ]).
  
  
  Thus, we have that 
      𝐊(β, ν)= 𝕏^⊤K𝕏.
   
  
  
  
  
  § APPENDIX B
  
  
  
  In this Appendix, we present the cumulants and derivatives needed to obtain the second order bias of β and ν. In addiction, we describe in details how to obtain B(β) and B(ν) from Eq. (<ref>). In order to present the cumulants in a summarized form, consider the following quantities:
  
  
      a_i    =   ψ^(1)(μ_i(1+ϕ_i)) - ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2), 
      
      		b_i    =   μ_i^2 ψ^(1)(μ_i(1+ϕ_i)) - (1 + μ_i)^2 ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) + ψ^(1)(ϕ_1 + 2),
      
      		c_i    =   ψ^(2)(μ_i(1+ϕ_i)) - ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2), 
      
      		d_i    =    (1 + μ_i)^2 ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2) - μ_i^2 ψ^(2)(μ_i(1+ϕ_i)), 
      
      		e_i    =    (1 + μ_i)^3 ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2) - μ_i^3 ψ^(2)(μ_i(1+ϕ_i)) - ψ^(2)(ϕ_i + 2).
  
  So, the cumulants needed here are presented as follow:
  
      κ_rs   =    - ∑_i=1^n (1 + ϕ_i)^2 a_i (μ_i/η_1i)^2 x_irx_is,
      κ_rS   =    - ∑_i=1^n (1 + ϕ_i) [a_i μ_i - ψ^(1)(μ_i(1 + ϕ_i) + ϕ_i + 2)] μ_i/η_1iϕ_i/η_2i x_irz_iS,
      κ_RS   =    - ∑_i=1^n b_i (ϕ_i/η_2i)^2 z_iRz_iS,
      κ_rsu   =    -∑_i=1^n (1 + ϕ_i)^2 {(1 + ϕ_i)c_i (μ_i/η_1i)^3 + 3a_i μ_i/η_1i^2 μ_i/η^2_1i} x_irx_isx_iu,
      κ_rsU   =    - ∑_i=1^n (1 + ϕ_i) {2a_i + (1 + ϕ_i)c_i μ_i - (1 + ϕ_i)ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2) }(μ_i/η_1i)^2 ϕ_i/η_2i x_irx_isz_iU
         +   ∑_i=1^n (1 + ϕ_i) {ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - a_iμ_i }^2 μ_i/η^2_1iϕ_i/η_2ix_irx_isz_iU,
      κ_rSU   =   ∑_i=1^n {ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - a_iμ_i }[(1 + ϕ_i)μ_i/η_1i^2 ϕ_i/η^2_2i + 2 ( ϕ_i/η_2i)^2 μ_i/η_1i] x_irz_iSz_iU
         +   ∑_i=1^n (1 + ϕ_i) d_i μ_i/η_1i( ϕ_i/η_2i)^2 x_irz_iSz_iU,
      κ_RSU   =   ∑_i=1^n {e_i (ϕ_i/η_2i)^3 - 3b_i ^2 ϕ_i/η^2_2iϕ_i/η_2i} z_iRz_iSz_iU.
  
  
  Taking the derivative of the cumulants with respect to the model parameters, we have
  
      κ_rs^(u)   =    - ∑_i=1^n (1 + ϕ_i)^2 {(1 +  ϕ_i)c_i (μ_i/η_1i)^3 + 2a_i μ_i/η_1i^2 μ_i/η^2_1i}x_irx_isx_iu,
      κ_rs^(U)   =    - ∑_i=1^n { (1 + ϕ_i)^2[c_i μ_i - ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2)] + 2(1 + ϕ_i)a_i }(μ_i/η_1i)^2 ϕ_i/η_2ix_irx_isz_iU,
      κ_RS^(u)   =   ∑_i=1^n {d_i(1 + ϕ_i) + 2 ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - 2a_i μ_i }μ_i/η_1i(ϕ_i/η_2i)^2 z_iRz_iSx_iu,
      κ_RS^(U)   =   ∑_i=1^n {e_i( ϕ_i/η_2i)^3 - 2b_i ϕ_i/η_2i^2 ϕ_i/η^2_2i}z_iRz_iSz_iU,
      κ_rS^(u)   =    - ∑_i=1^n (1+ϕ_i) {(1+ϕ_i)μ_i c_i + a_i - (1+ϕ_i)ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2)  }(μ_i/η_1i)^2 ϕ_i/η_2i x_irx_iuz_iS
         -   ∑_i=1^n (1 + ϕ_i) {a_i μ_i - ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) }^2 μ_i/η_1i^2ϕ_i/η_2ix_irx_iuz_iS,
      κ_rS^(U)   =   ∑_i=1^n {(1+ϕ_i)d_i - a_i μ_i + ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) }μ_i/η_1i(ϕ_i/η_2i)^2 x_irz_iSz_iU
         +   ∑_i=1^n (1 + ϕ_i) {ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - a_iμ_i  }μ_i/η_1i^2 ϕ_i/η^2_2i x_irz_iSz_iU.
  
  Consider now the following diagonal matrices:
  
      M_1    =   (- (1+ϕ_i)^2/2[ (1+ϕ_i)c_i (μ_i/η_1i)^3 + a_iμ_i/η_1i^2 μ_i/η_1i^2] ), 
      
      M_2    =   ( (1+ϕ_i)/2[-a_i μ_i + ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) ] ^2 μ_i/η^2_1iϕ_i/η_2i.
         -   (1+ϕ_i)^2/2. [ c_i μ_i  - ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2) ] (μ_i/η_1i)^2 ϕ_i/η_2i),
      
      M_3    =   ( - (1+ ϕ_i)/2{[2a_i + (1+ϕ_i)c_iμ_i - (1+ϕ_i)ψ^(2)(μ_i(1+ϕ_i) + ϕ_i + 2) ] (μ_i/η_1i)^2 ϕ_i/η_2i. . 
         +   .  [ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - a_i μ_i] ^2 μ_i/η^2_1iϕ_i/η_2i}),
      
      M_4    =   (1/2{[(1+ϕ_i)d_i + 2ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - 2a_iμ_i]μ_i/η_1i(ϕ_i/η_2i)^2 . . 
         -   . .  (1 + ϕ_i)[ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) - a_i μ_i] μ_i/η_1i^2 ϕ_i/η^2_2i}),
      
      M_5    =   ( (1 + ϕ_i)/2[d_i (ϕ_i/η_2i)^2μ_i/η_1i + ψ^(1)(μ_i(1+ϕ_i) + ϕ_i + 2) -a_iμ_i]μ_i/η_1i^2 ϕ_i/η^2_2i),
      
      M_6    =   (1/2[e_i(ϕ_i/η_2i)^3 - b_iϕ_i/η_2i^2 ϕ_i/η_2i^2] ).
  
  Considering the matrices defined above, we have that m_ij represent the ith element of M_j. So, we have
  
      κ_rs^(u) - 1/2κ_rsu   =   ∑_i=1^n m_1i x_irx_isx_iu,
      κ_Rs^(u) - 1/2κ_Rsu   =   ∑_i=1^n m_2i x_isx_iuz_iR,
      κ_rS^(u) - 1/2κ_rSu   =   ∑_i=1^n m_2i x_irx_iuz_iS,
      κ_rs^(U) - 1/2κ_rsU   =   ∑_i=1^n m_3i x_irx_isz_iU,
      κ_RS^(u) - 1/2κ_RSu   =   ∑_i=1^n m_4i x_iuz_iRz_iS,
      κ_Rs^(U) - 1/2κ_RsU   =   ∑_i=1^n m_5i x_isz_iUz_iR,
      κ_rS^(U) - 1/2κ_rSU   =   ∑_i=1^n m_5i x_irz_iSz_iU,
      κ_RS^(U) - 1/2κ_RSU   =   ∑_i=1^n m_6i z_iRz_iSz_iU,
  
  We now obtain the terms from (<ref>), presenting in detail the algebra to obtain the first term of this expression. To calculate the other terms, we follow the same logic.
  
      ∑_r,s,uκ^arκ^su{κ_rs^(u) - 1/2κ_rsu}   =   ∑_r,s,u(κ^arκ^su∑_i=1^n m_1ix_irx_isx_iu) 
         =   ∑_i=1^n m_1i(∑_r x_irκ^ar) (∑_s,u x_irκ^su x_iu)
         =   ∑_i=1^n m_1i(∑_r x_irκ^ar) δ_i^⊤ (X κ^ββ X^⊤)δ_i
         =   δ_a^⊤∑_i=1^n κ^a βX^⊤δ_i m_1iδ_i^⊤ (X κ^ββ X^⊤)						 δ_i
         =   δ_a^⊤κ^a β X^⊤ M_1 P_ββ,
  
  being κ^a β the matrix κ^ββ if a=1,…,p and κ^νβ if a=p+1,…,q, δ_a (δ_i) an n × 1 vector with a one in the ath (ith) position. Also, the vector P_ββ is presented in Section <ref>. Likewise, we have the remaining quantities expressed by
  
  
      ∑_R,s,uκ^aRκ^su{κ_Rs^(u) - 1/2κ_Rsu}   =   δ_a^⊤κ^a ν Z^⊤ M_2 P_ββ,
      ∑_r,S,uκ^arκ^Su{κ_rS^(u) - 1/2κ_rSu}   =   δ_a^⊤κ^a β X^⊤ M_2 P_βν,
      ∑_r,s,Uκ^arκ^sU{κ_rs^(U) - 1/2κ_rsU}   =   δ_a^⊤κ^a β X^⊤ M_3 P_βν,
      ∑_R,S,uκ^aRκ^Su{κ_RS^(u) - 1/2κ_RSu}   =   δ_a^⊤κ^a ν Z^⊤ M_4 P_βν,
      ∑_R,s,Uκ^aRκ^sU{κ_Rs^(U) - 1/2κ_RsU}   =   δ_a^⊤κ^a ν Z^⊤ M_5 P_βν,
      ∑_r,S,Uκ^arκ^SU{κ_rS^(U) - 1/2κ_rSU}   =   δ_a^⊤κ^a β X^⊤ M_5 P_νν,
      ∑_R,S,Uκ^aRκ^SU{κ_RS^(U) - 1/2κ_RSU}   =   δ_a^⊤κ^a ν Z^⊤ M_6 P_νν,
  
  where κ^a ν is the matrix κ^βν if a=1,…,p and κ^νν if a=p+1,…,q and the vectors P_βν and P_νν were presented in Section <ref>.
  
  
  
  