Multivariate Gamma Regression: Parameter Estimation, Hypothesis Testing, and Its Application

Abstract

Gamma distribution is a general type of statistical distribution that can be applied in various fields, mainly when the distribution of data is not symmetrical. When predictor variables also affect positive outcome, then gamma regression plays a role. In many cases, the predictor variables give effect to several responses simultaneously. In this article, we develop a multivariate gamma regression (MGR), which is one type of non-linear regression with response variables that follow a multivariate gamma (MG) distribution. This work also provides the parameter estimation procedure, test statistics, and hypothesis testing for the significance of the parameter, partially and simultaneously. The parameter estimators are obtained using the maximum likelihood estimation (MLE) that is optimized by numerical iteration using the Berndt–Hall–Hall–Hausman (BHHH) algorithm. The simultaneous test for the model’s significance is derived using the maximum likelihood ratio test (MLRT), whereas the partial test uses the Wald test. The proposed MGR model is applied to model the three dimensions of the human development index (HDI) with five predictor variables. The unit of observation is regency/municipality in Java, Indonesia, in 2018. The empirical results show that modeling using multiple predictors makes more sense compared to the model when it only employs a single predictor.

1. Introduction

Gamma distribution is one family of continuous probability distributions and generalizations of exponential distributions [1]. Nagar, Correa, and Gupta [2] mentioned that the gamma distribution function was first introduced by Swiss mathematician Leonhard Euler (1729). Because this function is considered important, many researchers have studied and developed it. Bhattacharya [3], among others, conducted a study on testing the homogeneity of the parameters (shape and scale) of the gamma distribution. Chen and Kotz [4] conducted a study on the probability density function (pdf) of gamma distribution with three parameters (shape, scale, and location). Many researchers also study and develop bivariate gamma distribution; among others are Schickedanz and Krause [5], who conducted a study on testing scale parameters from two gamma-distributed data using the generalized likelihood ratio (GLR). Nadarajah [6] studied the types of bivariate gamma distribution. Next, Nadarajah and Gupta [7] developed two new bivariate gamma distributions based on gamma and beta random variables. In addition, Mathai and Moschopoulos [8] discussed joint densities, product moments, conditional densities, and conditional moments that were developed from two bivariate gamma distributions.

One statistical method that can be applied to analyze the data that follow gamma distribution and its predictor variables is gamma regression. Gamma regression is a type of non-linear regression. A non-linear regression contains at least one parameter with a non-linear form [9,10]. The gamma regression with multiple responses is the so-called multivariate gamma regression (MGR).

The MGR model proposed in this article is the extension of the trivariate gamma regression (TGR) proposed by Rahayu, Purhadi, Sutikno, and Prastyo [11], which describes the theory of parameter estimation and its hypothesis testing. The MGR is developed based on multivariate gamma distribution with three parameters (shape, scale, and location). The supporting references about multivariate gamma distribution were written by Mathai and Moschopoulos [12], and Vaidyanathan and Lakshmi [13]. The parameter estimation method for MGR in this study uses maximum likelihood estimation (MLE). However, the solution cannot be obtained in the closed form. Therefore, a numerical method is needed to achieve the parameter estimator value. The numerical optimization used in this study is the Berndt–Hall–Hall–Hausman (BHHH).

Based on the previously mentioned background, the aims of this study are: (i) how to construct the MGR model, (ii) how to estimate the parameters, and (iii) how to test the significance of the model as well as the significance of the individual parameter. The last objective of this work is how to apply the proposed MGR model to real data. The case study used in this study includes the factors that affect the life expectancy index (first response), education index (second response), and expenditure index (third response), the three indexes that compose the human development dimensions. The unit of observation is the regency/municipality in Java, Indonesia, in 2018. The predictor variables include the percentage of households that have a private toilet, net enrollment rate of schooling, population density, the percentage of poor people, and the unemployment rate.

The rest of the article is organized as follows. Section 2 introduces the detail of the proposed MGR model. Section 3 and Section 4 explore the data and application, respectively. The last section contains conclusions and further research.

2. Multivariate Gamma Regression Model

Suppose ${\mathit{y}}_l$ is the response variables data $\left(y_{l1},y_{l2},\dots ,y_{lk}\right)$ that follows multivariate gamma distribution and ${\mathit{x}}_l$ is the corresponding predictor variables $\left(x_{l1},x_{l2},\dots ,x_{ls}\right)$, with sample size as n observations $\left(l=1,2,\dots ,n\right).$ In this section, we discuss the construction of the MGR model, its parameter estimation, and hypothesis testing. A short explanation about univariate gamma regression is introduced to make a smooth transition into the MGR model.

According to Balakrishnan and Wang [14], a random variable $Y$ follows univariate gamma distribution with three parameters $\left(\alpha ,\gamma ,\lambda \right),$ denoted by $Y~\mathrm{Gamma}\left(\alpha ,\gamma ,\lambda \right),$ with pdf formulated in Equation (1).

$$f\left(y\right)=\{\begin{array}{cc}\hfill \frac{1}{{\gamma }^{\alpha }\Gamma \left(\alpha \right)}{\left(y-\lambda \right)}^{\alpha -1}{e}^{-\frac{y-\lambda }{\gamma }};& \alpha ,\gamma ,\lambda >0, \lambda <y<\infty ,\hfill \\ \hfill ;& \mathrm{otherwise}.\hfill \end{array}$$

(1)

If $Y~\mathrm{Gamma}\left(\alpha ,\gamma ,\lambda \right)$, then the statistics are as follows [15,16].

$\mu =E\left(y\right)=\gamma \alpha +\lambda$, $Var\left(y\right)=\alpha {\gamma }^2$, $Stdev\left(y\right)=\sqrt{\alpha }\gamma ,$ and the skewness is ${\gamma }_1=\frac{2}{\sqrt{\alpha }}.$

Mathai and Moschopoulos (1992) defined the pdf as in Equation (2) for a pair of random variables $\left(Y_1,Y_2\right)$ that follows bivariate gamma distribution as:

$$f\left(y_1,y_2\right)=\frac{{\left(y_1-{\lambda }_1\right)}^{{\alpha }_1-1}{\left(y_2-y_1-{\lambda }_2\right)}^{{\alpha }_2-1}{e}^{-\frac{y_2-{\displaystyle \sum _{i=1}^2{\lambda }_i}}{\gamma }}}{{\gamma }^{{\alpha }_{2^{*}}}{\displaystyle \mathbf{\prod }_{i=1}^2\Gamma \left({\alpha }_i\right)}},$$

(2)

with ${\alpha }_i>0, \gamma >0, {\lambda }_i\in R, {\lambda }_1<y_1<\infty , {\lambda }_2<y_2<\infty , {\alpha }_{2^{*}}={\alpha }_1+{\alpha }_2, i=1,2,$ $f\left(y_1,y_2\right)=0$ for otherwise.

The mean for $Y_1$ and $Y_2$ are $E\left(Y_1\right)=\gamma {\alpha }_1+{\lambda }_1$ and $E\left(Y_2\right)=\gamma \left({\alpha }_1+{\alpha }_2\right)+{\lambda }_1+{\lambda }_2$, while the variances are $Var\left(Y_1\right)={\gamma }^2{\alpha }_1 \mathrm{and} Var\left(Y_2\right)={\gamma }^2\left({\alpha }_1+{\alpha }_2\right).$

Suppose there are k response variables; the pdf for random variables $\left(Y_1,Y_2,\dots ,Y_k\right)$ that follow multivariate gamma distribution (Mathai and Moschopoulos, 1992) is:

$$f\left(y_1,y_2,\dots ,y_k\right)=\frac{{\left(y_1-{\lambda }_1\right)}^{{\alpha }_1-1}{\left(y_2-y_1-{\lambda }_2\right)}^{{\alpha }_2-1}\cdots {\left(y_k-y_{k-1}-{\lambda }_k\right)}^{{\alpha }_k-1}{e}^{-\frac{y_k-{\displaystyle \sum _{i=1}^k{\lambda }_i}}{\gamma }}}{{\gamma }^{{\alpha }_{k^{*}}}{\displaystyle \mathbf{\prod }_{i=1}^k\Gamma \left({\alpha }_i\right)}},$$

(3)

with ${\alpha }_i>0,\gamma >0,{\lambda }_i\in R,$ ${\lambda }_1<y_1<\infty ,{\lambda }_2<y_2<\infty ,{\lambda }_k<y_k<\infty ,{\alpha }_{k^{*}}={\alpha }_1+{\alpha }_2+\cdots +{\alpha }_k,i=1,2,\dots k,$ otherwise $f\left(y_1,y_2,\dots ,y_k\right)=0$.

The mean and variance for $Y_i$ are $E\left(Y_i\right)=\gamma {\alpha }_{i^{*}}+{\lambda }_{i^{*}}$ and $Var\left(Y_i\right)={\gamma }^2{\alpha }_{i^{*}}$ with ${\alpha }_{i^{*}}={\alpha }_1+{\alpha }_2+\cdots +{\alpha }_i$ and ${\lambda }_{i^{*}}={\lambda }_1+{\lambda }_2+\cdots +{\lambda }_i.$ The MGR model can be stated in Equation (4).

$$E\left(Y_i\right)=\gamma {\alpha }_{i^{*}}+{\lambda }_{i^{*}}=e^{x^T{\beta }_i}, i=1,2,\dots ,k,$$

(4)

with ${\alpha }_{i^{*}}={\alpha }_1+{\alpha }_2+\cdots +{\alpha }_i, {\lambda }_{i^{*}}={\lambda }_1+{\lambda }_2+\cdots +{\lambda }_i.$

The pdf for the lth observation is formulated in Equation (5) which will be used to compose the likelihood function in Equation (6).

$$f\left(y_{l1},y_{l2},\dots ,y_{lk}\right)=\frac{{\left(y_{l1}-{\lambda }_1\right)}^{{\alpha }_1-1}{\left(y_{l2}-y_{l1}-{\lambda }_2\right)}^{{\alpha }_2-1}\cdots {\left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)}^{{\alpha }_k-1}{e}^{-\frac{y_{lk}-{\displaystyle \sum _{i=1}^k{\lambda }_i}}{\gamma }}}{{\gamma }^{{\alpha }_k^{*}}\Gamma \left({\alpha }_1\right)\Gamma \left({\alpha }_2\right)\cdots \Gamma \left({\alpha }_k\right)},$$

(5)

with ${\alpha }_i>0,\gamma >0,{\lambda }_i\in R,$ ${\lambda }_1<y_{l1}<\infty ,{\lambda }_2<y_{l2}<\infty ,{\lambda }_k<y_{lk}<\infty ,$

• ${\alpha }_1=\frac{e^{x_l^T{\beta }_1}-{\lambda }_1}{\gamma },{\alpha }_2=\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2}{\gamma },\dots ,{\alpha }_k=\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k}{\gamma },$
• ${\alpha }_{k^{*}}={\alpha }_1+{\alpha }_2+\cdots +{\alpha }_k=\frac{e^{x_l^T{\beta }_k}-{\lambda }_1-{\lambda }_2-\cdots -{\lambda }_k}{\gamma },$ otherwise $f\left(y_{l1},y_{l2},\dots ,y_{lk}\right)=0$.

Later, we discuss parameter estimation on MGR using MLE. The likelihood function constructed from Equation (5) is:

$$\begin{array}{l}L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)={\displaystyle \mathbf{\prod }_{l=1}^n f\left(y_{l1},y_{l2},\dots ,y_{lk}\right)}\hfill \\ ={\displaystyle \mathbf{\prod }_{l=1}^n\left(\frac{{\left(y_{l1}-{\lambda }_1\right)}^{{\alpha }_1-1}{\left(y_{l2}-y_{l1}-{\lambda }_2\right)}^{{\alpha }_2-1}\cdots {\left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)}^{{\alpha }_k-1}{e}^{-\frac{y_{lk}-{\displaystyle \sum _{i=1}^k{\lambda }_i}}{\gamma }}}{{\gamma }^{{\alpha }_k^{*}}\Gamma \left({\alpha }_1\right)\Gamma \left({\alpha }_2\right)\cdots \Gamma \left({\alpha }_k\right)}\right)},\hfill \end{array}$$

(6)

with values ${\alpha }_1,{\alpha }_2,{\alpha }_k,\mathrm{and} {\alpha }_{k^{*}}$ based on Equation (5).

The log-likelihood function from Equation (6) is:

$$\begin{array}{l}\log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)\hfill \\ ={\displaystyle \sum _{l=1}^n\log \left[\frac{{\left(y_{l1}-{\lambda }_1\right)}^{{\alpha }_1-1}{\left(y_{l2}-y_{l1}-{\lambda }_2\right)}^{{\alpha }_2-1}\cdots {\left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)}^{{\alpha }_k-1}{e}^{-\frac{y_{lk}-{\displaystyle \sum _{i=1}^k{\lambda }_i}}{\gamma }}}{{\gamma }^{{\alpha }_{k^{*}}} \Gamma \left({\alpha }_1\right) \Gamma \left({\alpha }_2\right)\cdots \Gamma \left({\alpha }_k\right)}\right]}.\hfill \end{array}$$

By substituting the values of ${\alpha }_1,{\alpha }_2,{\alpha }_k,\mathrm{and} {\alpha }_{k^{*}}$ according to Equation (5), the log-likelihood function is:

$$\begin{array}{l}\log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)\hfill \\ ={\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\beta }_1}-{\lambda }_1-\gamma }{\gamma }\log \left(y_{l1}-{\lambda }_1\right)}+{\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2-\gamma }{\gamma }\log \left(y_{l2}-y_{l1}-{\lambda }_2\right)}+\cdots +\hfill \\ \hspace{1em}{\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k-\gamma }{\gamma }\log \left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)}-{\displaystyle \sum _{l=1}^n\frac{y_{lk}-{\lambda }_1-{\lambda }_2-\cdots -{\lambda }_k}{\gamma }}-\\ \hspace{1em}{\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\beta }_k}-{\lambda }_1-{\lambda }_2-\cdots -{\lambda }_k}{\gamma }\log \gamma }-{\displaystyle \sum _{l=1}^n\log \Gamma \left(\frac{e^{x_l^T{\beta }_1}-{\lambda }_1}{\gamma }\right)}-{\displaystyle \sum _{l=1}^n\log \Gamma \left(\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2}{\gamma }\right)}-\cdots -\hfill \\ \hspace{1em}{\displaystyle \sum _{l=1}^n\log \Gamma \left(\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k}{\gamma }\right)}.\hfill \end{array}$$

(7)

In this article, the log value is based on e or natural logarithm. The first derivatives of the log-likelihood function for each parameter are as follows.

$$\begin{array}{c}\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial \gamma }={\displaystyle \sum _{l=1}^n\left(\frac{{\lambda }_1-e^{x_l^T{\beta }_1}}{{\gamma }^2}\log \left(y_{l1}-{\lambda }_1\right)\right)}+\\ {\displaystyle \sum _{l=1}^n\left(\frac{e^{x_l^T{\beta }_1}-e^{x_l^T{\beta }_2}+{\lambda }_2}{{\gamma }^2}\log \left(y_{l2}-y_{l1}-{\lambda }_2\right)\right)}+\cdots +{\displaystyle \sum _{l=1}^n\left(\frac{e^{x_l^T{\beta }_{k-1}}-e^{x_l^T{\beta }_k}+{\lambda }_k}{{\gamma }^2}\log \left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)\right)}+\\ {\displaystyle \sum _{l=1}^n\frac{y_{lk}}{{\gamma }^2}}-\frac{n{\lambda }_1}{{\gamma }^2}-\frac{n{\lambda }_2}{{\gamma }^2}-\frac{n{\lambda }_3}{{\gamma }^2}-(\frac{n\left(\log \gamma \right){\lambda }_1}{{\gamma }^2}-\frac{n{\lambda }_1}{{\gamma }^2}+\frac{n\left(\log \gamma \right){\lambda }_2}{{\gamma }^2}-\frac{n{\lambda }_2}{{\gamma }^2}+\frac{n\left(\log \gamma \right){\lambda }_3}{{\gamma }^2}-\frac{n{\lambda }_3}{{\gamma }^2}+\end{array}\\ {\displaystyle \sum _{l=1}^n\left(-\frac{\left(\log \gamma \right)e^{x_l^T{\beta }_k}}{{\gamma }^2}+\frac{e^{x_l^T{\beta }_k}}{{\gamma }^2}\right)})-{\displaystyle \sum _{l=1}^n\left(-\frac{1}{{\gamma }^2}\left(\Psi \left(\frac{e^{x_l^T{\beta }_1}-{\lambda }_1}{\gamma }\right)\right)\left(e^{x_l^T{\beta }_1}-{\lambda }_1\right)\right)}-\\ {\displaystyle \sum _{l=1}^n\left(-\frac{1}{{\gamma }^2}\left(\Psi \left(\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2}{\gamma }\right)\right)\left(e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2\right)\right)}-\cdots -\\ {\displaystyle \sum _{l=1}^n\left(-\frac{1}{{\gamma }^2}\left(\Psi \left(\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k}{\gamma }\right)\right)\left(e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k\right)\right)},\end{array}$$

(8)

$$\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial {\lambda }_1}={\displaystyle \sum _{l=1}^n\left(-\frac{\log \left(y_{l1}-{\lambda }_1\right)}{\gamma }-\frac{e^{x_l^T{\beta }_1}-{\lambda }_1-\gamma }{\gamma \left(y_{l1}-{\lambda }_1\right)}\right)}+\frac{n}{\gamma }+\frac{n\left(\log \gamma \right)}{\gamma }-\\ {\displaystyle \sum _{l=1}^n\left(-\frac{1}{\gamma }\Psi \left(\frac{e^{x_l^T{\beta }_1}-{\lambda }_1}{\gamma }\right)\right)},\end{array}$$

(9)

$$\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial {\lambda }_2}={\displaystyle \sum _{l=1}^n(-\frac{\log \left(y_{l2}-y_{l1}-{\lambda }_2\right)}{\gamma }-}\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2-\gamma }{\gamma \left(y_{l2}-y_{l1}-{\lambda }_2\right)})+\frac{n}{\gamma }+\\ \frac{n\left(\log \gamma \right)}{\gamma }-{\displaystyle \sum _{l=1}^n\left(-\frac{1}{\gamma }\Psi \left(\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2}{\gamma }\right)\right)},\end{array}$$

(10)

$$\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial {\lambda }_k}={\displaystyle \sum _{l=1}^n(-\frac{\log \left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)}{\gamma }-}\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k-\gamma }{\gamma \left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)})+\\ \frac{n}{\gamma }+\frac{n\left(\log \gamma \right)}{\gamma }-{\displaystyle \sum _{l=1}^n\left(-\frac{1}{\gamma }\Psi \left(\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k}{\gamma }\right)\right)},\end{array}$$

(11)

$$\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial {\mathit{\beta }}_1}={\displaystyle \sum _{l=1}^n\left(\left(\log \left(y_{l1}-{\lambda }_1\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_1}}{\gamma }\right)}-\\ {\displaystyle \sum _{l=1}^n\left(\left(\log \left(y_{l2}-y_{l1}-{\lambda }_2\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_1}}{\gamma }\right)}-{\displaystyle \sum _{l=1}^n\left(\left(\Psi \left(\frac{e^{x_l^T{\beta }_1}-{\lambda }_1}{\gamma }\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_1}}{\gamma }\right)}-\\ {\displaystyle \sum _{l=1}^n\left(-\left(\Psi \left(\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2}{\gamma }\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_1}}{\gamma }\right)},\end{array}$$

(12)

$$\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial {\mathit{\beta }}_2}={\displaystyle \sum _{l=1}^n\left(\left(\log \left(y_{l2}-y_{l1}-{\lambda }_2\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_2}}{\gamma }\right)}-\\ {\displaystyle \sum _{l=1}^n\left(\left(\Psi \left(\frac{e^{x_l^T{\beta }_2}-e^{x_l^T{\beta }_1}-{\lambda }_2}{\gamma }\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_2}}{\gamma }\right)}-{\displaystyle \sum _{l=1}^n\left(-\left(\Psi \left(\frac{e^{x_l^T{\beta }_3}-e^{x_l^T{\beta }_2}-{\lambda }_3}{\gamma }\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_2}}{\gamma }\right)},\end{array}$$

(13)

$$\begin{array}{l}\frac{\partial \log L\left(\gamma ,{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right)}{\partial {\mathit{\beta }}_k}={\displaystyle \sum _{l=1}^n\left(\left(\log \left(y_{lk}-y_{l\left(k-1\right)}-{\lambda }_k\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_k}}{\gamma }\right)}-\\ {\displaystyle \sum _{l=1}^n\left(\left(\log \gamma \right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_k}}{\gamma }\right)}-{\displaystyle \sum _{l=1}^n\left(\left(\Psi \left(\frac{e^{x_l^T{\beta }_k}-e^{x_l^T{\beta }_{k-1}}-{\lambda }_k}{\gamma }\right)\right)\frac{{\mathit{x}}_l^T{e}^{x_l^T{\beta }_k}}{\gamma }\right)},\end{array}$$

(14)

with $\Psi \left(z\right)$ = digamma function, which is the first derivative of gamma function, formulated with $\Psi \left(z\right)=\frac{\partial \left[\log \Gamma \left(z\right)\right]}{dz}=\frac{{\Gamma }^{\prime }\left(z\right)}{\Gamma \left(z\right)}.$

A maximum likelihood (ML) can be found by setting all the derivatives above to zero and solving the system. No closed-form solution to that system can be found. A numerical method is needed to obtain the solution, i.e., parameter estimate $\widehat{\gamma },{\widehat{\lambda }}_1,{\widehat{\lambda }}_2,\dots ,{\widehat{\lambda }}_k,{\widehat{\mathit{\beta }}}_1,{\widehat{\mathit{\beta }}}_2,\dots ,{\widehat{\mathit{\beta }}}_k.$ One of the numerical techniques that can be employed is the BHHH algorithm as follows.

• Step 1. Determine the initial value for ${\widehat{\mathsf{\theta }}}^{\left(0\right)}={\left[{\widehat{\gamma }}^{\left(0\right)} {{\widehat{\lambda }}_1}^{\left(0\right)} {{\widehat{\lambda }}_2}^{\left(0\right)} \dots {{\widehat{\lambda }}_k}^{\left(0\right)} {\widehat{\mathit{\beta }}}_1^{T\left(0\right)} {\widehat{\mathit{\beta }}}_2^{T\left(0\right)}\dots {\widehat{\mathit{\beta }}}_k^{T\left(0\right)}\right]}^T,$ where ${\widehat{\gamma }}^{\left(0\right)}>0,{{\widehat{\lambda }}_1}^{\left(0\right)},{{\widehat{\lambda }}_2}^{\left(0\right)},\dots ,{{\widehat{\lambda }}_k}^{\left(0\right)}\in R$ satisfies the constraints in Equation (5), and ${\widehat{\mathit{\beta }}}_1^{T\left(0\right)},{\widehat{\mathit{\beta }}}_2^{T\left(0\right)},\dots ,{\widehat{\mathit{\beta }}}_k^{T\left(0\right)}$ are obtained from the estimate of univariate gamma regression. The Hessian $\mathbf{H}\left(\widehat{\mathsf{\theta }}\right)$ in BHHH is approximated as the negative of the sum of the outer products of the gradients of individual observations. The gradient vector $\mathbf{g}\left(\widehat{\mathsf{\theta }}\right)$ is a vector with each element consisting of the first derivative of the log-likelihood function for each of the estimated parameters.
• Step 2. Determine the tolerance limit so that the BHHH iteration process stops. In this study, the tolerance value used is $\epsilon ={10}^{-8}.$
• Step 3. Start the BHHH iteration using the following formula.

  $${\widehat{\mathsf{\theta }}}^{\left(p+1\right)}={\widehat{\mathsf{\theta }}}^{\left(p\right)}-{\mathbf{H}}^{-1}\left({\widehat{\mathsf{\theta }}}^{\left(p\right)}\right)\mathbf{g}\left({\widehat{\theta }}^{\left(p\right)}\right),$$

  (15)

  with $p=0,1,2,\dots ,p^{*}.$
• Step 4. The iteration stops at the $p^{*}-th$ iteration if it satisfies $\Vert {\widehat{\mathsf{\theta }}}^{\left(p^{*}+1\right)}-{\widehat{\mathsf{\theta }}}^{\left(p*\right)}\Vert \le \epsilon$. When converging, the last iteration produces an estimator value for each parameter.

The null hypothesis on the MGR model is $H_0:{\beta }_{11}={\beta }_{21}=\cdots ={\beta }_{s1}={\beta }_{12}={\beta }_{22}=\cdots ={\beta }_{s2}=\cdots =$ ${\beta }_{1k}={\beta }_{2k}=\cdots ={\beta }_{sk}=0$ and alternative hypothesis H₁. At least one ${\beta }_{qi}\ne 0,$ with $q=1,2,\dots ,s,i=1,$ $2,\dots ,k.$ $\Omega =\left\{\gamma ,{\lambda }_1,\dots ,{\lambda }_k,{\mathit{\beta }}_1,{\mathit{\beta }}_2,\dots ,{\mathit{\beta }}_k\right\}$ is the set of parameters under the population. The $\omega =\{\gamma ,{\lambda }_1,$ $\dots ,{\lambda }_k,{\beta }_{01},{\beta }_{02},\dots ,{\beta }_{0k}\}$ is the set of parameters under the null hypothesis. The first derivatives of the log-likelihood function for each parameter under the null hypothesis are provided in Appendix A.

Proposition 1.

If $\Omega$ is a set of parameters under the population, $\omega$ is a set of parameters under the null hypothesis, and the hypothesis being used is the simultaneous test of MGR model, then the test statistic is $G^2=-2\log \Lambda =2\log L(\widehat{\Omega })-2\log L(\widehat{\omega }).$

A Corollary of Proposition 1:

The hypothesis being used in the simultaneous test of the MGR model in Section 2 can be stated in the following form: the null hypothesis is ${\mathit{\beta }}^{*}={\mathbf{0}}_{\left(sk\times 1\right)}$ and the alternative hypothesis is ${\mathit{\beta }}^{*}\ne {\mathbf{0}}_{\left(sk\times 1\right)},$ with ${\mathit{\beta }}^{*}={\left[\begin{array}{cccccc}{\mathit{\beta }}_1^{*T}& {\mathit{\beta }}_2^{*T}& \cdots & {\mathit{\beta }}_i^{*T}& \cdots & {\mathit{\beta }}_k^{*T}\end{array}\right]}^T$ and ${\mathit{\beta }}_i^{*}=\left[\begin{array}{cccc}{\beta }_{1i}& {\beta }_{2i}& \cdots & {\beta }_{si}\end{array}\right]$ for $i=1,2,\dots ,k$.

It is noted that ${\widehat{\mathsf{\theta }}}_{\Omega }$ and ${\widehat{\mathsf{\theta }}}_{\omega }$ are estimators that maximize the likelihood and the log-likelihood functions under the population and under the null hypothesis. The principle of the MLE method is to maximize the likelihood functions [17]. The following are test statistics for the hypothesis being used in the simultaneous test of the MGR model in Section 2.

$$\Lambda =\frac{L\left(\widehat{\omega }\right)}{L\left(\widehat{\Omega }\right)}<{\Lambda }_0,$$

(16)

where ${\Lambda }_0$ is a constant value between $0<{\Lambda }_0\le 1$.

$L\left(\widehat{\omega }\right)$ and $L\left(\widehat{\Omega }\right)$ in Equation (16) are:

$$L\left(\widehat{\omega }\right)={\displaystyle \mathbf{\prod }_{l=1}^n\left(\frac{{\left(y_{l1}-{\widehat{\lambda }}_1\right)}^{{\widehat{\alpha }}_{11}-1}{\left(y_{l2}-y_{l1}-{\widehat{\lambda }}_2\right)}^{{\widehat{\alpha }}_{22}-1}\cdots {\left(y_{lk}-y_{l\left(k-1\right)}-{\widehat{\lambda }}_k\right)}^{{\widehat{\alpha }}_{kk}-1}{e}^{-\frac{y_{lk}-{\displaystyle \sum _{i=1}^k{\widehat{\lambda }}_i}}{\widehat{\gamma }}}}{{\widehat{\gamma }}^{{\widehat{\alpha }}_{k{k}^{*}}}\Gamma \left({\widehat{\alpha }}_{11}\right)\Gamma \left({\widehat{\alpha }}_{22}\right)\cdots \Gamma \left({\widehat{\alpha }}_{kk}\right)}\right)},$$

and

$$L\left(\widehat{\Omega }\right)={\displaystyle \mathbf{\prod }_{l=1}^n\left(\frac{{\left(y_{l1}-{\widehat{\lambda }}_1\right)}^{{\widehat{\alpha }}_1-1}{\left(y_{l2}-y_{l1}-{\widehat{\lambda }}_2\right)}^{{\widehat{\alpha }}_2-1}\cdots {\left(y_{lk}-y_{l\left(k-1\right)}-{\widehat{\lambda }}_k\right)}^{{\widehat{\alpha }}_k-1}{e}^{-\frac{y_{lk}-{\displaystyle \sum _{i=1}^k{\widehat{\lambda }}_i}}{\widehat{\gamma }}}}{{\widehat{\gamma }}^{{\widehat{\alpha }}_{k^{*}}}\Gamma \left({\widehat{\alpha }}_1\right)\Gamma \left({\widehat{\alpha }}_2\right)\cdots \Gamma \left({\widehat{\alpha }}_k\right)}\right)},$$

(17)

with ${\widehat{\alpha }}_{11}=\frac{e^{{\widehat{\beta }}_{01}}-{\widehat{\lambda }}_1}{\widehat{\gamma }},{\widehat{\alpha }}_{22}=\frac{e^{{\widehat{\beta }}_{02}}-e^{{\widehat{\beta }}_{01}}-{\widehat{\lambda }}_2}{\widehat{\gamma }},\dots ,{\widehat{\alpha }}_{kk}=\frac{e^{{\widehat{\beta }}_{0k}}-e^{{\widehat{\beta }}_{0\left(k-1\right)}}-{\widehat{\lambda }}_k}{\widehat{\gamma }},$

• ${\widehat{\alpha }}_{k{k}^{*}}={\widehat{\alpha }}_{11}+{\widehat{\alpha }}_{22}+\cdots +{\widehat{\alpha }}_{kk}=\frac{e^{{\widehat{\beta }}_{0k}}-{\widehat{\lambda }}_1-{\widehat{\lambda }}_2-\cdots -{\widehat{\lambda }}_k}{\widehat{\gamma }},$
• $\begin{array}{l}{\widehat{\alpha }}_1=\frac{e^{x_l^T{\widehat{\beta }}_1}-{\widehat{\lambda }}_1}{\widehat{\gamma }},{\widehat{\alpha }}_2=\frac{e^{x_l^T{\widehat{\beta }}_2}-e^{x_l^T{\widehat{\beta }}_1}-{\widehat{\lambda }}_2}{\widehat{\gamma }},\dots ,{\widehat{\alpha }}_k=\frac{e^{x_l^T{\widehat{\beta }}_k}-e^{x_l^T{\widehat{\beta }}_{k-1}}-{\widehat{\lambda }}_k}{\widehat{\gamma }},\mathrm{and}\\ {\widehat{\alpha }}_{k^{*}}={\widehat{\alpha }}_1+{\widehat{\alpha }}_2+\cdots +{\widehat{\alpha }}_k=\frac{e^{x_l^T{\widehat{\beta }}_k}-{\widehat{\lambda }}_1-{\widehat{\lambda }}_2-\cdots -{\widehat{\lambda }}_k}{\widehat{\gamma }}.\end{array}$

Based on Equation (17), $\frac{L\left(\widehat{\omega }\right)}{L\left(\widehat{\Omega }\right)}$ is difficult to simplify. To simplify the calculation, the test statistics in Equation (16) are expressed in a form equivalent to:

$${(\Lambda )}^{-2}={\left(\frac{L\left(\widehat{\omega }\right)}{L\left(\widehat{\Omega }\right)}\right)}^{-2}={\left(\frac{L\left(\widehat{\Omega }\right)}{L\left(\widehat{\omega }\right)}\right)}^2.$$

(18)

The application of natural logarithms in Equation (18) obtains the following test statistics.

$$G^2=-2\log \Lambda =-2\log \left(\frac{L\left(\widehat{\omega }\right)}{L\left(\widehat{\Omega }\right)}\right)=2\log \left(\frac{L\left(\widehat{\Omega }\right)}{L\left(\widehat{\omega }\right)}\right)=2\log L\left(\widehat{\Omega }\right)-2\log L\left(\widehat{\omega }\right),$$

(19)

with $\log L\left({\widehat{\Omega }}_{MGR}\right)={\displaystyle \sum _{l=1}^n\log }\left(f\left(y_{l1},y_{l2},\dots ,y_{lk}|{\widehat{\Omega }}_{MGR}\right)\right)$

$$\begin{array}{l}\log L\left({\widehat{\Omega }}_{MGR}\right)={\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\widehat{\beta }}_1}-{\widehat{\lambda }}_1-\widehat{\gamma }}{\widehat{\gamma }}\log \left(y_{l1}-{\widehat{\lambda }}_1\right)}+{\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\widehat{\beta }}_2}-e^{x_l^T{\widehat{\beta }}_1}-{\widehat{\lambda }}_2-\widehat{\gamma }}{\widehat{\gamma }}\log \left(y_{l2}-y_{l1}-{\widehat{\lambda }}_2\right)}+\cdots +\\ \begin{array}{l}{\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\widehat{\beta }}_k}-e^{x_l^T{\widehat{\beta }}_{k-1}}-{\widehat{\lambda }}_k-\widehat{\gamma }}{\widehat{\gamma }}\log \left(y_{lk}-y_{l\left(k-1\right)}-{\widehat{\lambda }}_k\right)}-{\displaystyle \sum _{l=1}^n\frac{y_{lk}-{\widehat{\lambda }}_1-{\widehat{\lambda }}_2-\cdots -{\widehat{\lambda }}_k}{\widehat{\gamma }}}-\\ {\displaystyle \sum _{l=1}^n\frac{e^{x_l^T{\widehat{\beta }}_k}-{\widehat{\lambda }}_1-{\widehat{\lambda }}_2-\cdots -{\widehat{\lambda }}_k}{\widehat{\gamma }}\log \widehat{\gamma }}-{\displaystyle \sum _{l=1}^n\log \Gamma \left(\frac{e^{x_l^T{\widehat{\beta }}_1}-{\widehat{\lambda }}_1}{\widehat{\gamma }}\right)}-{\displaystyle \sum _{l=1}^n\log \Gamma \left(\frac{e^{x_l^T{\widehat{\beta }}_2}-e^{x_l^T{\widehat{\beta }}_1}-{\widehat{\lambda }}_2}{\widehat{\gamma }}\right)}-\cdots -\end{array}\\ {\displaystyle \sum _{l=1}^n\log \Gamma \left(\frac{e^{x_l^T{\widehat{\beta }}_k}-e^{x_l^T{\widehat{\beta }}_{k-1}}-{\widehat{\lambda }}_k}{\widehat{\gamma }}\right)},\end{array}$$

$$\log L\left({\widehat{\omega }}_{MGR}\right)={\displaystyle \sum _{l=1}^n\log }\left(f\left(y_{l1},y_{l2},\dots ,y_{lk}|{\widehat{\omega }}_{MGR}\right)\right)$$

$$\begin{array}{cc}\hfill \log L\left({\widehat{\omega }}_{MGR}\right)=& {\displaystyle \sum _{l=1}^n\frac{e^{{\widehat{\beta }}_{01}}-{\widehat{\lambda }}_1-\widehat{\gamma }}{\widehat{\gamma }}\log \left(y_{l1}-{\widehat{\lambda }}_1\right)}+{\displaystyle \sum _{l=1}^n\frac{e^{{\widehat{\beta }}_{02}}-e^{{\widehat{\beta }}_{01}}-{\widehat{\lambda }}_2-\widehat{\gamma }}{\widehat{\gamma }}\log \left(y_{l2}-y_{l1}-{\widehat{\lambda }}_2\right)}+\cdots +\hfill \\ \hfill & {\displaystyle \sum _{l=1}^n\frac{e^{{\widehat{\beta }}_{0k}}-e^{{\widehat{\beta }}_{0\left(k-1\right)}}-{\widehat{\lambda }}_k-\widehat{\gamma }}{\widehat{\gamma }}\log \left(y_{lk}-y_{l\left(k-1\right)}-{\widehat{\lambda }}_k\right)}-{\displaystyle \sum _{l=1}^n\frac{y_{lk}-{\widehat{\lambda }}_1-{\widehat{\lambda }}_2-\cdots -{\widehat{\lambda }}_k}{\widehat{\gamma }}}-\hfill \\ \hfill & \frac{e^{{\widehat{\beta }}_{0k}}-{\widehat{\lambda }}_1-{\widehat{\lambda }}_2-\cdots -{\widehat{\lambda }}_k}{\widehat{\gamma }}\log \widehat{\gamma }-\log \Gamma \left(\frac{e^{{\widehat{\beta }}_{01}}-{\widehat{\lambda }}_1}{\widehat{\gamma }}\right)-\log \Gamma \left(\frac{e^{{\widehat{\beta }}_{02}}-e^{{\widehat{\beta }}_{01}}-{\widehat{\lambda }}_2}{\widehat{\gamma }}\right)-\cdots -\hfill \\ \hfill & \log \Gamma \left(\frac{e^{{\widehat{\beta }}_{0k}}-e^{{\widehat{\beta }}_{0\left(k-1\right)}}-{\widehat{\lambda }}_k}{\widehat{\gamma }}\right).\hfill \end{array}$$

Proposition 2.

Based on Proposition 1, the distribution of test statistics $G^2$ is Chi-square with sk degrees of freedom, which can be written as follows.

$$G^2=2\log L(\widehat{\Omega })-2\log L(\widehat{\omega })\stackrel{d}{\to }{\chi }_{sk}^2,n\to \infty .$$

A Corollary of Proposition 2:

If ${\widehat{\mathsf{\theta }}}_{\Omega }$ is an estimator that maximizes the likelihood and the log-likelihood functions under the population, ${\widehat{\mathsf{\theta }}}_{\omega }$ is an estimator that maximizes the likelihood and the log-likelihood functions under the null hypothesis, based on Equation (19), so:

$$\begin{array}{cc}\hfill G^2& =2\log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-2\log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)\hfill \\ \hfill & =2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-\log L\left({\mathsf{\theta }}_{\omega }\right)\right)-2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)-\log L\left({\mathsf{\theta }}_{\omega }\right)\right).\hfill \end{array}$$

(20)

$\mathrm{Log}L\left({\mathsf{\theta }}_{\omega }\right)$ function can be approached by Taylor’s second-degree expansion around ${\widehat{\mathsf{\theta }}}_{\Omega }$ as follows.

$$\log L\left({\mathsf{\theta }}_{\omega }\right)\approx \log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)+\mathbf{g}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\Omega }\right)-\frac{1}{2}{\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\Omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\Omega }\right),$$

(21)

with $\mathbf{g}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)={\frac{\partial \log L\left({\mathsf{\theta }}_{\Omega }\right)}{\partial {\mathsf{\theta }}_{\Omega }}|}_{{\mathsf{\theta }}_{\Omega }={\widehat{\mathsf{\theta }}}_{\Omega }}=\mathbf{0}$ and $\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)=-{\frac{{\partial }^2\log L\left({\mathsf{\theta }}_{\Omega }\right)}{\partial {\mathsf{\theta }}_{\Omega }\partial {\left({\mathsf{\theta }}_{\Omega }\right)}^T}|}_{{\mathsf{\theta }}_{\Omega }={\widehat{\mathsf{\theta }}}_{\Omega }}.$

Because $\mathbf{g}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)=\mathbf{0},$ then Equation (21) becomes:

$$\begin{array}{l}\log L\left({\mathsf{\theta }}_{\omega }\right)\approx \log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-\frac{1}{2}{\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\Omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\Omega }\right)\\ 2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-\log L\left({\mathsf{\theta }}_{\omega }\right)\right)\approx {\left({\widehat{\mathsf{\theta }}}_{\Omega }-{\mathsf{\theta }}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\widehat{\mathsf{\theta }}}_{\Omega }-{\mathsf{\theta }}_{\omega }\right).\end{array}$$

(22)

$\mathrm{Log}L\left({\mathsf{\theta }}_{\omega }\right)$ function can be approached by Taylor’s second-degree expansion around ${\widehat{\mathsf{\theta }}}_{\omega }$ as follows.

$$\log L\left({\mathsf{\theta }}_{\omega }\right)\approx \log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)+\mathbf{g}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\omega }\right)-\frac{1}{2}{\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\omega }\right).$$

(23)

Because $\mathbf{g}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)=\mathbf{0},$ then Equation (23) becomes:

$$\begin{array}{l}\log L\left({\mathsf{\theta }}_{\omega }\right)\approx \log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)-\frac{1}{2}{\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\mathsf{\theta }}_{\omega }-{\widehat{\mathsf{\theta }}}_{\omega }\right)\\ 2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)-\log L\left({\mathsf{\theta }}_{\omega }\right)\right)\approx {\left({\widehat{\mathsf{\theta }}}_{\omega }-{\mathsf{\theta }}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\widehat{\mathsf{\theta }}}_{\omega }-{\mathsf{\theta }}_{\omega }\right).\end{array}$$

(24)

Based on Equations (22) and (24), the test statistics on Equation (20) can be stated as follows.

$$G^2=2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-\log L\left({\mathsf{\theta }}_{\omega }\right)\right)-2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)-\log L\left({\mathsf{\theta }}_{\omega }\right)\right)$$

(25)

$$G^2\approx {\left({\widehat{\mathsf{\theta }}}_{\Omega }-{\mathsf{\theta }}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\widehat{\mathsf{\theta }}}_{\Omega }-{\mathsf{\theta }}_{\omega }\right)-{\left({\widehat{\mathsf{\theta }}}_{\omega }-{\mathsf{\theta }}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\widehat{\mathsf{\theta }}}_{\omega }-{\mathsf{\theta }}_{\omega }\right).$$

Equation (25) can be simplified by outlining the quadratic form of ${\left({\widehat{\mathsf{\theta }}}_{\Omega }-{\mathsf{\theta }}_{\omega }\right)}^T\left[\mathbf{I}\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)\right]\left({\widehat{\mathsf{\theta }}}_{\Omega }-{\mathsf{\theta }}_{\omega }\right)$, so we obtained:

$$2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-\log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)\right)\approx {\widehat{\mathit{\beta }}}^{*T}\left(\left[{\mathbf{I}}_{11}\right]-\left[{\mathbf{I}}_{12}\right]{\left[{\mathbf{I}}_{22}\right]}^{-1}\left[{\mathbf{I}}_{21}\right]\right){\widehat{\mathit{\beta }}}^{*}\approx {\widehat{\mathit{\beta }}}^{*T}{\left[{\mathbf{I}}^{11}\right]}^{-1}{\widehat{\mathit{\beta }}}^{*}.$$

(26)

From Equation (26), this can be obtained:

$${\widehat{\mathit{\beta }}}^{*}\stackrel{d}{\to }N\left(\mathbf{0},{\left[{\mathbf{I}}^{11}\right]}_{\left(sk\times sk\right)}\right),n\to \infty ,$$

(27)

$${\left[{\mathbf{I}}^{11}\right]}^{-\frac{1}{2}}{\widehat{\mathit{\beta }}}^{*}\stackrel{d}{\to }N\left(\mathbf{0},{\mathbf{I}}_{sk}\right).$$

(28)

Based on Equation (28), the quadratic form given by Equation (26) distributed Chi-square with sk degrees of freedom is:

$$\begin{array}{cc}\hfill 2\left(\log L\left({\widehat{\mathsf{\theta }}}_{\Omega }\right)-\log L\left({\widehat{\mathsf{\theta }}}_{\omega }\right)\right)\approx & {\left[{\left[{\mathbf{I}}^{11}\right]}^{-\frac{1}{2}}{\widehat{\mathit{\beta }}}^{*}\right]}^T\left[{\left[{\mathbf{I}}^{11}\right]}^{-\frac{1}{2}}{\widehat{\mathit{\beta }}}^{*}\right]\hfill \\ \hfill & ={\mathbf{z}}^T\mathbf{z}\stackrel{d}{\to }{\chi }_{sk}^2,n\to \infty ,\hfill \end{array}$$

(29)

with $\mathbf{z}={\left[{\mathbf{I}}^{11}\right]}^{-\frac{1}{2}}{\widehat{\mathit{\beta }}}^{*}\stackrel{d}{\to }N\left(\mathbf{0},{\mathbf{I}}_{sk}\right),n\to \infty .$

sk is a vector dimension ${\mathit{\beta }}^{*}$ or the difference between the number of parameter sets under the population with the number of parameter sets under the null hypothesis, symbolized by $n(\Omega )-n\left(\omega \right).$

Proposition 3.

The critical area for testing the hypothesis of the MGR model regression parameters simultaneously with regard to Equation (16) is:

$$\begin{array}{cc}\hfill \alpha & =P\left(\Lambda <{\Lambda }_0\right)\hfill \\ \hfill & =P\left(-2\log \Lambda >-2\log {\Lambda }_0\right)\hfill \\ \hfill & =P\left(G^2>c_1\right), \mathrm{with} c_1=-2\log {\Lambda }_0\hfill \\ \hfill & =P\left(G^2>{\chi }_{\alpha ,sk}^2\right)\hfill \\ \hfill & =P\left(G^2>{\chi }_{\alpha ,n(\Omega )-n\left(\omega \right)}^2\right).\hfill \end{array}$$

(30)

Based on Proposition 2 and Proposition 3, the decision to reject the null hypothesis is made if $G^2>{\chi }_{\alpha ;df}^2$, with $df=n(\Omega )-n\left(\omega \right),$ $n(\Omega )$ is the number of parameters under the population, and $n\left(\omega \right)$ is the number of parameters under the null hypothesis.

The null hypothesis for the partial test is $H_0:{\beta }_{qi}=0$, whereas the alternative is $H_1:{\beta }_{qi}\ne 0$, with $q=1,2,\dots ,s,i=1,2,\dots k$. According to Pawitan [18], the test statistic is stated in Equation (31).

$$Z=\frac{{\widehat{\beta }}_{qi}}{SE({\widehat{\beta }}_{qi})},$$

(31)

with $SE({\widehat{\beta }}_{qi})=\sqrt{\widehat{\mathrm{var}}({\widehat{\beta }}_{qi})}$. The $\widehat{\mathrm{var}}({\widehat{\beta }}_{qi})$ is diagonal elements that correspond to the $-{\mathbf{H}}^{-1}\left(\widehat{\mathsf{\theta }}\right)$ matrix. The null hypothesis is rejected if $\left|Z\right|>Z_{\alpha /2}.$

3. Data and Method

The parameter estimation and hypothesis testing on MGR were done based on the following steps. The MGR model was specified based on the pdf in Equation (5) for n observations, l = 1, 2, …, n, to construct the likelihood and the log-likelihood functions. The first derivative of the log-likelihood function for each parameter was computed, then equalized to zero. If the solutions were closed-form, then the parameter estimators were obtained. Otherwise, numerical optimization was needed. As shown in the previous section, the solution for parameter optimization was not closed-form, such that the BHHH algorithm was employed in this work.

The overall test for MGR’s significance was done using the maximum likelihood ratio test (MLRT). The test statistic was formulated in Equation (19). Meanwhile, the partial test for individual parameter significance in MGR was done using the Wald test [18]. Its test statistics are provided in Equation (31). The proposed MGR model, along with its parameter estimation and hypothesis testing, was applied on real data as an application of this study.

This study used secondary data obtained from Statistics Indonesia. The data used were three response variables, i.e., the life expectancy index, education index, and expenditure index, with six predictor variables: percentage of households that have a private toilet, net enrollment rate of schooling, population density, percentage of poor people, and unemployment rate. The data were observed for 119 regencies/municipalities in Java, Indonesia, in the year 2018.

4. Application on Human Development Dimensions Data

First, testing the gamma distribution was done using the Kolmogorov–Smirnov (KS) test. The null hypothesis is the data that follows the gamma distribution against the alternative hypothesis that data does not follow the gamma distribution. The test statistic value of the KS test for each response variable is presented in

Table 1. Gamma distribution test with Kolmogorov–Smirnov (KS) for $\alpha =0.05$.

Response	${\mathit{D}}_{\mathit{n}}$	${\mathit{D}}_{\left(0,05\right)}$	p-Value
Y1	0.118	0.124	0.066
Y2	0.107	0.124	0.123
Y3	0.065	0.124	0.667

. In this paper, the goodness of fit is done univariately as the test for multivariate gamma distribution is not available yet. The test for that is another extensive work that is not covered in this paper. Once each response follows gamma distribution, we assume the multiresponses data follow a multivariate gamma distribution. This assumption is the limitation of this work, such that the proposed model can be applied to real data without delay.

Each response variable has $D_n<D_{\left(0.05\right)}$ and p-value > α. The test concludes not to reject the null hypothesis, meaning that the data of life expectancy index (Y₁), the education index (Y₂), and the expenditure index (Y₃) follow the gamma distribution. Therefore, as our research limitation, as mentioned previously, the three response variables are assumed to follow MG distribution.

To support our assumption, we calculated the correlation between the pair of the response variables to show there are dependencies among responses. The correlation coefficients for each pair are as follows: (i) Y₁ and Y₂ is 0.398 with p-value close to zero, (ii) Y₁ and Y₃ is 0.324 with p-value close to zero, (iii) Y₂ and Y₃ is 0.818 with p-value close to zero. The correlation coefficient between education index (Y₂) and expenditure index (Y₃) is stronger than the other pairs. To find out whether there is dependency among the response variables, one can use Bartlett’s test of sphericity so that the data are feasible for multivariate analysis. This test has statistic value ${\chi }^2=148.735$ and p-value = 2.22 × 10⁻¹⁶. The ${\chi }^2>{{\chi }^2}_{3;0,05}$ (or 7.815) and p-value < α, and alpha is 0.05. The decision is to reject the null hypothesis (Pearson correlation matrix not equal to an identity matrix), which means the correlation between the response variables is significant in the multivariate sense. Therefore, the data analysis needs to be done in a multivariate way using the MGR model.

We also tested the multicollinearity among the predictor variables. The variance inflation factor (VIF) value for each predictor variable is 1.358 (for X₁), 1.350 (X₂), 1.560 (X₃), 1.849 (X₄), and 1.211 (X₅). The VIF value for each of the predictor variables is less than ten which shows there is no multicollinearity among the predictor variables.

In

Table 2. Description of data.

Variables	Mean	SD	Coefficient of Variation	Min	Max
Life expectancy index (Y1)	0.806	0.042	5.200	0.680	0.890
Education index (Y2)	0.632	0.077	12.210	0.470	0.850
Expenditure index (Y3)	0.735	0.067	9.140	0.620	0.960
Percentage of households that have a private toilet (X1)	80.215	9.710	12.100	37.820	98.010
Net enrollment rate of schooling (X2)	63.579	10.650	16.750	34.220	89.460
Population density (X3)	3298	4493	136.250	278	19757
Percentage of poor people (X4)	9.623	4.211	43.750	1.680	21.210
Unemployment rate (X5)	5.337	2.430	45.530	1.430	12.770

, the mean values for response variables Y₁, Y₂, and Y₃ are 0.806, 0.632, and 0.735. Although Y₁, Y₂, and Y₃ have mean values that do not differ greatly, they are not necessarily of the same quality; it depends on the size of the spread of the data. One measure of data distribution that can be used is the coefficient of variation (CoV). The CoV for Y₁, Y₂, and Y₃ are 5.200, 12.210, and 9.140, respectively. The CoV for education index (Y₂) is the highest among others, which means that the variable is more heterogeneous. The CoV for predictor variables X₁, X₂, X₃, X₄, and X₅ are, respectively, 12.100, 16.750, 136.250, 43.750, and 45.530. The CoV for population density (X₃) is the highest among other predictor variables as its range is also the biggest one.

The dependency between response and predictor variables can be shown visually by the matrix plot, as exhibited in Figure 1. The correlation between X₃ and X₄ (−0.585) is stronger than the correlation between X₄ and the other predictor variables, even stronger than other pairs. The correlation between X₁ and X₅ (−0.017) is weakest compared to the correlation of other couples. There are indications that the relationship is non-linear between X₃ with the response variable and the other predictor variables. For the correlation between response and predictor variables, log(Y₁) has the strongest correlation with X₁ (0.434) compared to other predictors. The log(Y₂) and X₃ have the strongest correlation (0.705) compared with other predictors, while the correlation of log(Y₃) and X₃ is the strongest one (0.744). This value shows that log expenditure index and population density has the strongest relationship among other pairs.

To find out which predictor variables significantly predicted response variables, we employed the MGR model.

Table 3. Parameter estimation of multivariate gamma regression (MGR) model with a single predictor.

Parameter	Estimate	Standard Error	z	p-Value
$\left(Y_1,Y_2,Y_3\right)~X_1$
${\beta }_{01}$	−0.422284	0.646166	−0.654	0.513
${\beta }_{11}$	0.002758	0.007970	0.346	0.729
${\beta }_{02}$	−0.762038	1.120002	−0.680	0.496
${\beta }_{12}$	0.004297	0.016611	0.259	0.796
${\beta }_{03}$	−0.437765	0.760203	−0.576	0.565
${\beta }_{13}$	0.002339	0.011250	0.208	0.835
$\left(Y_1,Y_2,Y_3\right)~X_2$
${\beta }_{01}$	−0.305102	0.518863	−0.588	0.557
${\beta }_{21}$	0.001896	0.008076	0.235	0.814
${\beta }_{02}$	−0.881184	0.575235	−1.532	0.126
${\beta }_{22}$	0.007363	0.008563	0.860	0.390
${\beta }_{03}$	−0.430093	0.857503	−0.502	0.616
${\beta }_{23}$	0.002789	0.014685	0.190	0.849
$\left(Y_1,Y_2,Y_3\right)~X_3$
${\beta }_{01}$	−0.205148	0.000421	−486.716	0.000 *
${\beta }_{31}$	0.000002	0.000054	0.029	0.977
${\beta }_{02}$	−0.483898	0.000262	−1846.353	0.000 *
${\beta }_{32}$	0.000020	0.000021	0.930	0.352
${\beta }_{03}$	−0.306820	0.000209	−1465.698	0.000 *
${\beta }_{33}$	0.000016	0.000014	1.132	0.258
$\left(Y_1,Y_2,Y_3\right)~X_4$
${\beta }_{01}$	−0.166290	0.581543	−0.286	0.775
${\beta }_{41}$	−0.002307	0.019731	−0.117	0.907
${\beta }_{02}$	−0.238049	1.200557	−0.198	0.843
${\beta }_{42}$	−0.019448	0.046895	−0.415	0.678
${\beta }_{03}$	−0.136018	1.044931	−0.130	0.896
${\beta }_{43}$	−0.012480	0.035650	−0.350	0.726
$\left(Y_1,Y_2,Y_3\right)~X_5$
${\beta }_{01}$	−0.168892	0.555574	−0.304	0.761
${\beta }_{51}$	−0.006666	0.057555	−0.116	0.908
${\beta }_{02}$	−0.460389	0.862713	−0.534	0.594
${\beta }_{52}$	0.006020	0.092371	0.065	0.948
${\beta }_{03}$	−0.279763	0.831607	−0.336	0.737
${\beta }_{53}$	0.003974	0.040046	0.099	0.921

* Significant at $\alpha =10\%$.

presents the ML estimates of the MGR model with a single predictor and their corresponding standard errors, z score, and p-value. Every single predictor does not affect any response variables. Only the intercepts when the MGR model employs X₃ as a single predictor are significant.

The MGR model with a single predictor (for example, the X₅) for the life expectancy index, education index, and expenditure index is obtained as follows.

$$\begin{array}{l}{\widehat{\mu }}_{l1}=\exp \left(-0.168892-0.006666X_{l5}\right),\\ {\widehat{\mu }}_{l2}=\exp \left(-0.460389+0.006020X_{l5}\right),\\ {\widehat{\mu }}_{l3}=\exp \left(-0.279763+0.003974X_{l5}\right).\end{array}$$

As summarized in Table 3, it is shown that all predictor variables are not significant. For comparison, we also did MGR modeling with multiple predictors.

Table 4. Parameter estimation of MGR model with multiple predictors.

Parameter	Estimate	Standard Error	z	p-Value
Life expectancy index (Y1)
${\beta }_{01}$	−0.353421	0.000119	−2969.383	0.000 **
${\beta }_{11}$	0.002005	0.002987	0.671	0.502
${\beta }_{21}$	0.000653	0.003494	0.187	0.852
${\beta }_{31}$	0.000004	0.000029	0.124	0.902
${\beta }_{41}$	−0.000469	0.002831	−0.166	0.868
${\beta }_{51}$	−0.009502	0.005252	−1.809	0.070 **
Education index (Y2)
${\beta }_{02}$	−0.606408	0.000120	−5069.156	0.000 **
${\beta }_{12}$	0.000207	0.003800	0.055	0.956
${\beta }_{22}$	0.004706	0.004952	0.950	0.342
${\beta }_{32}$	0.000012	0.000018	0.677	0.498
${\beta }_{42}$	−0.011729	0.004125	−2.843	0.004 **
${\beta }_{52}$	−0.011194	0.006636	−1.687	0.092 **
Expenditure index (Y3)
${\beta }_{03}$	−0.274026	0.000101	−2723.533	0.000 **
${\beta }_{13}$	0.000779	0.002337	0.333	0.739
${\beta }_{23}$	0.000298	0.002621	0.114	0.910
${\beta }_{33}$	0.000013	0.000014	0.943	0.346
${\beta }_{43}$	−0.006542	0.003125	−2.093	0.036 **
${\beta }_{53}$	−0.007994	0.003900	−2.050	0.040 **

** Significant at $\alpha =10\%$.

presents the ML estimates of the MGR model with multiple predictors along with their corresponding standard errors, z score, and p-value.

The estimate of the scale parameter is 0.649423, with its standard error 0.000028. The estimate of ${\lambda }_1$, the location parameter for Y₁, is 0.670845 (standard error 0.006884); meanwhile, the estimate for ${\lambda }_2$ is −0.309362, with standard error 0.006507, and for ${\lambda }_3$ is 0.000468 (standard error 0.006530). The significant parameters are the scale parameter $\gamma$, the location parameter for Y₁ and Y₂, respectively, and ${\lambda }_1$ and ${\lambda }_2$, as their p-values are less than $\alpha =10\%.$ The estimate of each parameter corresponding to each predictor is summarized in Table 4. Therefore, the MGR model for the life expectancy index, education index, and expenditure index is obtained as follows.

$$\begin{array}{l}{\widehat{\mu }}_{l1}=\exp \left(-0.353421+0.002005X_{l1}+0.000653X_{l2}+0.000004X_{l3}-0.000469X_{l4}-0.009502X_{l5}\right)\\ {\widehat{\mu }}_{l2}=\exp \left(-0.606408+0.000207X_{l1}+0.004706X_{l2}+0.000012X_{l3}-{0.011729}_{l4}-0.011194X_{l5}\right)\\ {\widehat{\mu }}_{l3}=\exp \left(-0.274026+0.000779X_{l1}+0.000298X_{l2}+0.000013X_{l3}-0.006542X_{l4}-0.007994X_{l5}\right)\end{array}$$

The Akaike information criterion (AIC) value is −63.903, and the corrected Akaike information criterion (AICc) value is −53.361. To know the average squared difference between the estimated and the actual values, one can use the mean square error (MSE). The MSEs for the life expectancy index, education index, and expenditure index are 0.001, 0.002, and 0.003, respectively. As the MSE is an unbiased estimator of variance, the MSE value is expected to be not much different from the variance of each response variable, i.e., 0.002 (expectancy index), 0.006 (education index), and 0.004 (expenditure index).

We can perform the simultaneous test for the model’s significance using Wilk’s likelihood ratio statistics derived based on the MLRT. The test statistic value is 46.682, and the value of the Chi-square table with 15 degrees of freedom and $\alpha =10\%$ is 22.307. The test statistical value is larger than the value of the Chi-square table; therefore, the decision is to reject the null hypothesis. It means that the five predictor variables have a significant effect on the response variables simultaneously. To find out the predictor variables that partially affect the response variable, one can use test statistics in Equation (31). From Table 4, it can be seen that the significant predictor variable that influences the life expectancy index is the unemployment rate (X₅); meanwhile the education index and expenditure index are significantly affected by the percentage of poor people (X₄) and unemployment rate (X₅).

Based on the results of MGR modeling with a single predictor (Table 3) and multiple predictors (Table 4), it can be determined the differences in the coefficient signs only happen for X₅ in response to Y₂ and Y₃, as shown in

Table 5. The difference in coefficient signs and significance of the parameters in the MGR model.

Response Variables	Predictor Variables	MGR Modeling
		Multiple Predictors	Single Predictor
			X1	X2	X3	X4	X5
Y1	X1	+	+
	X2	+		+
	X3	+			+
	X4	−				−
	X5	− ***					−
Y2	X1	+	+
	X2	+		+
	X3	+			+
	X4	− ***				−
	X5	− ***					+
Y3	X1	+	+
	X2	+		+
	X3	+			+
	X4	− ***				−
	X5	− ***					+

*** Significant at $\alpha =10\%$.

. We can also find the supports of this evidence from the matrix plot in Figure 1, that individually, X₅ has a negative relationship with Y₁, while it has positive dependencies with Y₂ and Y₃. On the other side, the X₄ has a stronger negative individual relationship with all responses. Therefore, when X₄ and X₅ are used as predictors together in the MGR model, the sign of X₅ changes as there is a significant correlation (−0.391 with p-value < 0.05) between X₄ and X₅, where X₄ affects the response of Y₂ and Y₃ is stronger than X₅.

.

Recall the VIF value for X₅ is 1.211, which is small. This value means that there is a weak relationship between X₅ and (X₁, X₂, X₃, and X₄). However, there is a significant correlation between X₄ alone and X₅. In the MGR with multiple predictors, the positive sign for X₅ will not change if X₄ also has a positive sign for responses X₂ and X₃. Unfortunately, that is not the case. The correlation between X₄ and X₂ has a different sign compared with the correlation between X₅ and X₂. The sign of X₄ and X₅ in MGR with multiple predictors can change depending on its correlation with the response variable. The same explanation pertains to response X₃.

Life expectancy index (Y₁) has a negative association with the percentage of poor people (X₄), even though it is not significant for regency/municipality in Java. This finding means that an increase in life expectancy index is not affected by the percentage of poor people in Java. The education index (Y₂) and expenditure index (Y₃) have a significant negative dependency on the percentage of poor people in Java.

The predictions resulting from the MGR model are expected to be close to the actual values. The closer those two values, the narrower the spread, as displayed in Figure 2. It can be seen that fitting values for Y₂ and Y₃ are better than those of Y₁. This result is also supported by significant predictors, as reported in Table 4. The life expectancy index has one significant predictor, while the other two responses have two significant predictors that increase their coefficients of determination.

5. Conclusions

The proposed MGR model has been developed along with its parameter estimation and hypothesis testing. The solution of parameter estimation using MLE is not closed-form such that it is optimized numerically using the BHHH algorithm. The MLRT and Wald tests are employed for testing the model’s significance and the individual parameter, respectively. The proposed MGR model is applied to model the three dimensions of the human development index (HDI) with five predictor variables. The empirical results show that modeling using multiple predictors makes more sense compared to the model when it only employs a single predictor. When multiple predictors are used in the MGR model, there is a possibility that the sign of a particular parameter changes compared to when it is employed alone. This is a common problem that arises in modeling caused by collinearity among predictors. This issue can be overcome in future work.