Geographically Weighted Multivariate Logistic Regression Model and Its Application

This study investigates the geographically weighted multivariate logistic regression (GWMLR) model, parameter estimation, and hypothesis testing procedures. The GWMLR model is an extension to the multivariate logistic regression (MLR) model, which has dependent variables that follow a multinomial distribution along with parameters associated with the spatial weighting at each location in the study area. The parameter estimation was done using the maximum likelihood estimation and NewtonRaphson methods, and the maximum likelihood ratio test was used for hypothesis testing of the parameters. The performance of the GWMLR model was evaluated using a real dataset and it was found to perform better than the MLR model.


Introduction
Over the past decade, most research on geographically weighted regression (GWR) models has been focused on applications that contain two or more correlated responses (multivariate). Harini et al. [1,2] introduced the multivariate GWR (MGWR) model and demonstrated the parameter estimation and hypothesis test procedures using the restricted maximum likelihood estimation (RMLE) and maximum likelihood ratio test (MLRT) methods, respectively. The form and properties of the estimated errors variance-covariance parameters of the MGWR model using the MLE and weighted least squares methods were investigated [3]. Triyanto et al. [4,5] introduced the geographically weighted multivariate Poisson regression (GWMPR) model. The estimator of the GWMPR model parameters was obtained through the MLE with the Newton-Raphson iterative method, and the test statistic for hypothesis tests was determined by the MLRT method. Suyitno et al. [6] discussed the estimation of the geographically weighted trivariate Weibull regression (GWTWR) model using the MLE and Newton-Raphson methods. The geographically weighted multivariate t regression (GWMtR) model was introduced by Sugiarti et al. [7]. The MLE method and the expectationmaximization algorithm were applied to estimate the GWMtR model parameters. In [8], a new method to determine model conformity between the multivariate nonparametric truncated spline GWR model and the multivariate nonparametric truncated spline (global regression) was employed.
The responses of the multivariate GWR models in previous research were in the form of quantitative data. However, in many applications within various fields of research, the responses include not only quantitative data but also qualitative (categorical) data. Therefore, in this study, we propose the geographically weighted multivariate logistic regression (GWMLR) model. The GWMLR model is the extension of the geographically weighted bivariate logistic regression (GWBLR) proposed by Fathurahman et al. [9]. The GWMLR model has been developed from the geographically weighted logistic regression (GWLR) model proposed by Atkinson et al. [10]. The GWLR model is a combination of the GWR model [11] and the binary logistic regression model. The GWMLR model in this study is used to explain the spatial associations between two correlated categorical dependent variables with one or more independent variables, where each of the dependent variables has two categories. Similar to the methods in the works of Harini et al. [2], Triyanto et al. [4,5], Suyitno et al. [6], and Sifriyani et al. [8], the MLE and MLRT methods were used in the modeling and applying of the GWMLR model. The MLE method was used to estimate the parameters, and the statistical test for the significance of the parameters was determined by the MLRT method. The GWMLR model performance was evaluated using the factors that influence the public health development index and human development index of districts and cities in Kalimantan Island, Indonesia.

Multivariate Logistic Regression Model.
A multivariate logistic regression (MLR) explains the relationship between two or more correlated categorical dependent variables with one or more independent variables. In this study, the MLR model had two correlated categorical dependent variables, and each dependent variable had two categories. Let Y 1 and Y 2 be the two dependent variables. Y 1 and Y 2 each can have one of the two values (0 or 1). Let y = ½Y 11 Y 10 Y 01 Y 00 T be a vector of dependent variables of the MLR model. The elements of y have the probabilities of η 11 , η 10 , η 01 , and η 00 , respectively, which are presented in Table 1.

Geographically Weighted Multivariate Logistic Regression Model
The GWMLR model is an extension of the MLR model, used when the regression parameter depends on the spatial weight of all locations in the study area. The spatial weight, commonly used by the kernel functions [11,17], depends on both the Euclidean distance and an optimal bandwidth. The GWMLR model in this study is expressed as follows: Abstract and Applied Analysis where x i is a vector of independent variables at location i for i = 1, 2, ⋯, n and γ 1 ðu i , v i Þ, γ 2 ðu i , v i Þ, and γ 3 ðu i , v i Þ are the vectors of the GWMLR model parameters at location i. The vectors of independent variables and parameters at location i are η * 1 ðx i Þ and η * 2 ðx i Þ are the marginal probabilities of dependent variables at location i and are formulated as follows: ψ 2 ðx i Þ is called the odds ratio of dependent variables at location i and can be determined by where

Model Selection
In this study, the best model was selected using the three most common information criteria, which are Akaike's information criterion (AIC), the corrected AIC (AICC), and the Bayesian information criterion (BIC). All three information criteria formulas are as follows: where Lðb γðu i , v i Þ, i = 1, 2, ::, nÞ is the log-likelihood function of an estimated model, evaluated at the maximum likelihood estimator of the parameters at all locations ðnÞ; K is the number of effective parameters in the model at all locations, where X and Wðu i , v i Þ are the matrix of independent variables and spatial weighting, respectively. The best model has the lowest values of AIC, AICC, and BIC.

Results and Discussion
5.1. Estimation of the GWMLR Model Parameters. The parameters of the GWMLR model can be obtained using the maximum likelihood method. The likelihood function is as follows: where vector of the GWMLR model parameters. Let ðη * rsi ðx i ÞÞ y rsi = ðη * rsi Þ y rsi for r, s = 0, 1 ; i = 1, 2, ⋯, n; then, the likelihood function in Equation (16) is formulated by The maximum likelihood estimator of the GWMLR model parameters can be determined by maximizing the likelihood function in Equation (17) or by maximizing the natural logarithm of the likelihood function (log-likelihood). The log-likelihood function is given by Based on the GWR method, the spatial weighting function is presented as a log-likelihood. Let w ij be the spatial weighting function for each location ðu i , v i Þ, where w ij = w j ðu i , v i Þ and i, j = 1, 2, ⋯, n. The log-likelihood function is defined as follows: 3 Abstract and Applied Analysis where w ij is a fixed kernel bi-square [16] and formulated by where d ij is the Euclidean distance from i to j, and b is called an optimal bandwidth for the parameter estimation of the model at location i. In this study, the optimal bandwidth is determined by the cross-validation (CV) method [18,19]. The formula of the CV method in this study is as follows: where y rsi is the observation of the dependent variables with category values of Y 1 = r and Y 2 = s at location i, and b η * rs,≠i ðbÞ is the estimated value of the joint probabilities of the dependent variables that have category values of Y 1 = r, Y 2 = s, and bandwidth b with location i omitted from the estimation process. The optimal bandwidth has the lowest value of CV.
Theorem 1 obtains the maximum likelihood estimator of the GWMLR model parameters. Proof. Based on the GWMLR model in Equations (6)-(8), let T are formed. We then determine the derivative of ∂ξ/∂η * . The vector of η * has four elements, whereas the vector of ξ only has three elements. To obtain a symmetrical matrix of ∂ξ/∂η * , let ξ 0 = ln η * ++ with η *  2   6  6  6  6  6  6  6  6  6  6  4   3   7  7  7  7  7  7  7  7  7  7  5 , ð22Þ , ð23Þ The log-likelihood function in Equation (9) is maximized by determining the first-order partial derivative of the likelihood function, then equating to zero. The first-order partial derivative of the log-likelihood function with respect to the parameters of γðu i , v i Þ is as follows: where Δ j and Δ * j for j = 1, 2, ⋯, n in Equation (23). The details of the first-order partial derivative of the loglikelihood function with respect to the parameters of γðu i , v i Þ in Equations (24)-(26) are presented in the appendix.
The first-order partial derivative of the log-likelihood with respect to the parameters of γðu i , v i Þ in Equations (24)-(26) produces an implicit form. This result shows that the estimator of the GWMLR model parameters cannot be obtained analytically and requires a numerical approach. The numerical approach by the Newton-Raphson method was used to obtain the maximum likelihood estimator of the GWMLR model parameters. The Newton-Raphson method requires the gradient vector and the Hessian matrix, which are formulated as follows: 4 Abstract and Applied Analysis H After obtaining the gradient vector and Hessian matrix, the Newton-Raphson iteration process is carried out with the following formula: ðu i , v i Þ and ε is a low positive number.

Hypothesis Test.
Hypothesis testing on the GWMLR model parameters was performed and included the similarity test, simultaneous test, and partial test. The similarity test was used to find the differences between the MLR and GWMLR models. The simultaneous test was used to simultaneously obtain the significant influence of the independent variables on the dependent variables. The simultaneous test was also used to obtain at least one of the independent variables that have a significant influence on the dependent variables. The partial test was used to obtain the partially significant influence of the independent variables on the dependent variables. The similarity test was conducted using the hypotheses: The statistical test is as follows: where The statistical test in Equation (32) followed an asymptotically standard normal distribution. Therefore, the null hypothesis ðH 0 Þ in Equation (30) is rejected at the level of significance ðαÞ when the value of the V statistic in Equation (32) falls into the rejection region (i.e., jVj > Z α/2 ).
The next test presented is the simultaneous test. The hypothesis of this test is formulated as follows: Theorem 2 is presented next for the simultaneous test.

Theorem 2.
The statistical test of the hypothesis in the simultaneous test is as follows: Proof. The G 2 statistic can be obtained by the maximum likelihood ratio test method. The initial step of this method determines the parameters set under the population.
Analogously, in Equations (17) and (18), the likelihood function under the population is as follows: However, the maximum likelihood estimator of the GWMLR model parameters was obtained in Theorem 1. Therefore, the maximum log-likelihood function under the 5 Abstract and Applied Analysis population is as follows: The parameters set under the null hypothesis are Analogously, in Equations (38) and (39), the likelihood function under the null hypothesis is The maximum log-likelihood function under the null hypothesis is as follows: where the joint probabilities of b η * * 11i , b η * * 10i , b η * * 01i , and b η * * 00i are obtained by with Based on the maximum likelihood ratio test method, the statistical test of the hypothesis in Equation (34) is formulated as follows: The likelihood ratio statistic ðG 2 Þ in Equation (44) has an asymptotic chi-square distribution, where the degree of freedom is the difference between the number of model parameters under the population and the number of model parameters under the null hypothesis is v = 3kn. Therefore, at an α significance level, we reject the null hypothesis when the G 2 value falls into the rejection region (i.e., G 2 > χ 2 ðα,vÞ ). The last hypothesis test of the GWMLR model parameters is the partial test. The hypothesis is The statistical test for the hypothesis in Equation (45) is given by where SEðb is the diagonal elements of −½Hðγ∧ðu i , v i ÞÞ −1 and Hðb γðu i , v i ÞÞ is derived in Equation (28). The Z statistic in Equation (47) has an asymptotic standard normal distribution. Therefore, the null hypothesis in Equation (45) is rejected when the value of the Z statistic falls into the rejection region (i.e., jZj > Z α/2 ).

5.3.
Application. The GWMLR model was applied to real data, which included the public health development index (PHDI) and the human development index (HDI) for the districts/cities in Kalimantan Island, Indonesia, in 2013. The PHDI describes the quality of health and the progress of health development of the districts/cities and provinces in Indonesia. The PHDI is used to prioritize districts/cities 6 Abstract and Applied Analysis that need assistance in health development [20]. The HDI is an index that measures the basic dimensions of human development in the districts/cities [21]. The PHDI data were provided by the Ministry of Health, Indonesia. The National Bureau of Statistics Indonesia provided the HDI data and independent variables. The variables in this study consist of two dependent variables and two independent variables. PHDI ðY 1 Þ and HDI ðY 2 Þ are dependent variables. The PHDI has two categories: 0 if the PHDI value of districts/cities is less than the PHDI value of Indonesia, and 1 if the PHDI value of districts/cities is greater than or equal to the PHDI value of Indonesia. The HDI has two categories: 0 if the HDI value of districts/cities is less than the HDI value of Indonesia, and 1 if the HDI value of districts/cities is greater than or equal to the HDI value of Indonesia. The poverty rate ðX 1 Þ and economic growth ðX 2 Þ are the independent variables. The unit observation is the districts/cities in Kalimantan Island, Indonesia, in 2013. The sample size is 55, consisting of 46 districts and 9 cities. The computation in this study is performed using MATLAB and the econometrics toolbox [22].
The implementation of the GWMLR model for the PHDI and HDI of districts/cities in Kalimantan Island began by creating a 2 × 2 contingency table for the observed frequencies of the dependent variables and for determining their proportion and correlation. The observed frequencies of the dependent variables are reported in Table 2. Table 2 shows that 13 districts/cities had PHDI and HDI values greater than or equal to the PHDI and HDI values of Indonesia, and 26 districts/cities had PHDI and HDI values less than the PHDI and HDI values of Indonesia. We also see that three districts/cities had a PHDI value greater than or equal to the PHDI value of Indonesia and an HDI value less than the HDI value of Indonesia. Finally, 13 districts/cities had a PHDI value less than the PHDI value of Indonesia and an HDI value greater than or equal to the HDI value of Indonesia. The odds ratio value of the dependent variables was 8.6667, which shows that the dependent variables were positively correlated. Therefore, the dependent variables of PHDI and HDI were appropriate for the MLR and GWMLR models.
The parameter estimation obtained a total of 55 GWMLR models. The optimal bandwidth value of the fixed kernel bisquare weighting function was 4.8572, with a CV value of 90.3673. The descriptive statistics of the maximum likelihood estimator values of the 55 GWMLR models for modeling the PHDI and HDI of districts/cities in Kalimantan Island are given in Table 3.
The similarity evaluation between the MLR and GWMLR models was carried out using the statistical test in Equation (32). The hypothesis was formulated as follows: The statistical test value was 376.0917, and the Z α/2 value at a 10% significance level was 1.6449. Therefore, the null hypothesis was rejected, and we concluded that the MLR and GWMLR models were significantly different.
The next test was the simultaneous test, and the hypothesis was formulated as follows:

Abstract and Applied Analysis
The likelihood ratio statistic value was 401.6335, and the χ 2 ð0:1;330Þ value was 363.3222. This result indicated that the likelihood ratio was statistically significant at a 10% significance level. Therefore, the null hypothesis was rejected, and we concluded that the poverty rate and economic growth were simultaneously significantly influencing the PHDI and HDI of the districts/cities in Kalimantan Island.
The partial test was used to obtain the independent variables that significantly influence the PHDI and HDI of the districts/cities in Kalimantan Island. We show the results of the parameter estimation, standard errors, and statistical test values of the partial test for the GWMLR model of Lamandau District in Table 4.
The hypothesis of the partial test for the Lamandau District GWMLR model parameters was as follows: The statistical test value of the estimated parameter value of γ 23 ðu 21 , v 21 Þ in Table 4 was not statistically significant at a 10% significance level. Therefore, we concluded that the economic growth was not a partially significant influence on the PHDI and HDI of Lamandau District.
The GWMLR model of Lamandau District is expressed as follows: The results of the estimating and hypothesis testing of the GWMLR model parameters show that not all of the independent variables had a significant influence on the PHDI and HDI of the districts/cities in Kalimantan Island. Therefore, a significant number of parameters differed for each district/city. The groupings of the independent variables that had a significant influence on the PHDI and HDI of districts/cities are presented in Table 5.
The performance of the GWMLR model was evaluated using the AIC, AICC, and BIC in Equations (13)- (15). The values of the AIC, AICC, and BIC for the MLR and GWMLR models are presented in Table 6.
The AIC, AICC, and BIC values of the GWMLR model in Table 6 are lower compared with the MLR model. This result shows that the GWMLR model is more accurate than the MLR model. Therefore, the GWMLR model is best for modeling relationships between the dependent variables (the PHDI and HDI) and the independent variables (poverty rate and the economic growth) of districts/cities in Kalimantan Island, Indonesia, in 2013.

Conclusions
The GWMLR model is capable of evaluating the relationships between two correlated categorical dependent variables with one or more independent variables that depend on the spatial weighting function at each location in the study area. The spatial weighting function is an essential tool for parameter estimation and hypothesis testing in the modeling of GWMLR. Therefore, the fixed kernel bi-square was used, as it relates to the Euclidean distance and optimal bandwidth. The cross-validation method was applied to obtain the optimal bandwidth. The GWMLR model parameters were estimated using the maximum likelihood method. The maximum likelihood estimators have an implicit form and were obtained by the Newton-Raphson method. Hypothesis testing of the GWMLR model included a similarity test, simultaneous test, and partial test. The similarity test was used to obtain a significant difference between the MLR and GWMLR models. The statistical test of the similarity test had an asymptotic standard normal distribution. The  Abstract and Applied Analysis simultaneous test was used to obtain the simultaneously significant influence of the independent variables on the dependent variables. The statistical test of the simultaneous test had an asymptotic chi-square distribution. The partial test was used to obtain the partially significant influence of the independent variables on the dependent variables. The statistical test of the partial test had an asymptotic standard normal distribution. The performance of the GWMLR model was evaluated with the factors influencing the public health development index and the human development index of districts/cities in Kalimantan Island, Indonesia, in 2013. The GWMLR model was found to be better than the MLR model in this context. Some improvements and future work on modeling of GWMLR are possible. Firstly, the spatial weighting function used in this study involves only kernel bi-square. Other spatial weighting functions of kernel functions could be used to further the GWMLR model performance, such as the Gaussian, exponential, or tri-cube functions. Secondly, this study is limited to two dependent variables that only have two categories. Having more dependent variables with more than two categories, which could be either multinomial or ordinal, should also be considered for future work.