Avaliação da dependência espacial na modelagem do desempenho da segurança viária em zonas de tráfego

<jats:p>Uma técnica comumente usada no processo de modelagem do Desempenho da Segurança Viária (DSV), no nível de planejamento, são os Modelos Lineares Generalizados (MLG) assumindo a distribuição binomial negativa para os erros. Limitações dessa técnica, por não considerar os efeitos espaciais, têm sido contornadas com a utilização de modelos espaciais locais a partir de técnicas de regressão espacial como a Regressão de Poisson Geograficamente Ponderada (RPGP). Este trabalho tem por objetivo apresentar uma análise comparativa entre modelos de previsão de acidentes globais não espaciais e locais espaciais para a estimação do DSV agregado ao nível de zonas de tráfego em Fortaleza/CE. Foram calibrados modelos para a variável dependente acidentes totais e acidentes com vítimas e os resultados mostraram que os modelos RPGP apre-sentaram melhor desempenho que os MLG nas medidas de ajustes e na redução da autocorrelação espacial dos resíduos, sendo capazes de captar a heterogeneidade espacial da frequência dos acidentes.</jats:p>

section_content: Local spatial regression differs from global spatial regression by analyzing the relationship between variables in a specific way for each unit of study instead of combining them. In fact, regions j closer region i have greater influence in the estimates of the regression coefficients than when those regions are far apart. Having a specific adjustment for each region, the final result is a better representation of the process as a whole. The reason of to use this analysis is the violation of the assumption of stationarity demanded by global models, which allows one to attribute the same relation between variables for all units of study. The nonstationarity or spatial heterogeneity is the process where responses vary with location, area or other characteristics of the spatial units (Anselin, 1988). In this way, Geographically Weighted Regression (GWR) (see Fotheringham et al. (2002) ) is used to local modeling. 

However, this technique is used when the distribution of the data is Gaussian. In many applications, the dependent variable represents a count, which makes classic GWR inappropriate to model this type of data. Poisson and negative binomial are the most adequate distributions for the modeling of count data, and Geographically Weighted Poisson Regression (GWPR) was developed by Nakaya et al. (2005). Silva and Rodrigues (2014) developed GWR using negative binomial distribution, named GWNBR, and they showed that GWNBR incorporates GWPR. Thus, the main objective of this paper is to show a SAS ® macro to estimate the parameters of the GWNBR model.
section_name: INTRODUCTION

section_content: The most used global Negative Binomial Regression (NB-2) considers a logarithm link function. Parameterizing this model in terms of ߤ ‫ݐ/‬, where ‫ݐ‬ is an off set variable, ߤ is the predicted mean, ߙ is the parameter of overdispersion, ߚ is the parameter related to the explicative variable ‫ݔ‬, for ݇ = 1, …, ‫,ܭ‬ and ‫ݕ‬ is the j-th dependent variable for ݆ = 1, …, ݊ we have. 

where NB represents Negative Binomial. 

GWNBR is an extension of the global or non-spatial model (1), which allows the spatial variation of the parameters ߚ and ߙ. Without limiting the functional form of this variation, GWNBR produces nonparametric surfaces of the parameter estimates. This local model is described as the following:

where ‫ݑ(‬, ‫ݒ‬ ) are the locations (coordinates) of the data points j, for ݆ = 1, …, ݊. 

The parameter estimation of the global model ( 1 ) is performed interactively with the combination of the NR and IRLS methods, i. e., using ߙ ො, which is estimated by the NR method, we estimate the vector β β β β using the IRLS method. Thus, from this new β β β β, we update ߙ ො, and so forth until convergence is obtained. However, modifications in the NR and IRLS algorithms are necessary to incorporate the local variations. 

As shown in Silva and Rodrigues (2014), the analytical solution for the local log-likelihood of GWNBR is given by

where is an ݊ × ݇ matrix of the explicative variables

and ‫ݑ(‬, ‫ݒ‬ ) () is an ݊ × ݊ GLM diagonal weighting matrix for iteration m and location i

Silva and Rodrigues (2014) used the IRLS method with the observed Fisher Information Matrix, and the elements ܽ () (݆ = 1, …, ݊) of GWNBR ( 6 ) are the following:

As in GWPR (Nakaya et al., 2005), the covariance matrix of the parameter estimates can be estimated by

where

and the elements of ‫ݑ(‬, ‫ݒ‬ ) and ‫ݑ(‬, ‫ݒ‬ ) are given by ( 5 ) and (6), respectively. 

After an estimate of ‫ݑ(ࢼ‬, ‫ݒ‬ ) is obtained, the parameters ߙ will be estimated using the NR method based on the local log-likelihood given by (Silva and Rodrigues, 2014 )

Maximizing the local log-likelihood (11) using the univariate NR method, we obtain

where ܷ () and ‫ܪ‬ () are the first and second derivates of the local log-likelihood with respect to ‫ݎ‬ (), i. e.,

where ߰(. ) and ߰′(. ) are the digamma and trigamma functions, respectively, which are given by

Using delta method, Silva and Rodrigues (2014) found that, for a function ݃(. (Casella and Berger, 2011). If, under certain conditions, the estimators that maximize the local likelihood are asymptotically normal, unbiased and consistent (Staniswalis, 1989), then

because ‫)‪‬ݎ(ݎܸܽ‬ = ‫ܪ/1-‬, where ‫ܪ‬ is given in ( 14 ), and

Thus, for each regression point i, the NR and IRLS algorithms are used alternately until the parameter estimates achieve convergence. 

To complete the fitting of the model, it is necessary to estimate the bandwidth of the chosen kernel function, which helps determine the weights ‫݀(ݓ‬ ). One possibility could be to estimate it such that it minimizes the corrected AIC criterion (AICc):

where k is the effective number of parameters and ‫,ࢼ(ܮ‬ ࢻ) is the log-likelihood of GWNBR shown in Equation ( 11 ). 

The effective number of parameters of GWNBR can be written as ݇ = ݇ ଵ + ݇ ଶ, where ݇ ଵ and ݇ ଶ are the effective number of parameters due to ࢼ and ࢻ, respectively. Following the method developed by Nakaya et al. ( 2005 ), ݇ ଵ is given by the trace of the matrix, which is given by elements

where is the j-th row of. 

However, to date, it has not been possible to estimate ݇ ଶ, i. e., the contribution of the surface of ࢻ on the effective number of parameters of the model. Consequently, Silva and Rodrigues (2014) opted to estimate the bandwidth using the cross-validation criterion given by (Fotheringham et al., 2002) :

where ‫ݕ‬ ො ஷ (ܾ) is the estimated value for point j, omitting the observation j and b is the bandwidth. 

Note that the indetermination of ݇ ଶ does not prevent the adjustment of GWNBR. However, this indetermination makes it difficult to compare models because the complexity of the GWNBR, which is given by the effective number of parameters, is unknown.
section_name: A BASIC OUTLINE OF GWNBR

section_content: To avoid the difficult associate with the estimation of ݇ ଶ, Silva and Rodrigues (2014) proposed the Geographically Weighted Negative Binomial Regression with ߙ global methodology, which is namely GWNBRg. In this model, the spatial variation is allowed only for ‫ݑ(ࢼ‬, ‫ݒ‬ ), i. e.,

where the parameters are the same as in Equation ( 2 ). 

In the GWNBRg model, the estimation of the parameter ߙ is made globally, i. e., Silva and Rodrigues (2014) assumed that all of the parameters in the model are stationary, and we estimate a global overdispersion ߙ ො to be used in the local estimates ‫ݑ(ࢼ‬, ‫ݒ‬ ). Consequently, they proposed that the estimate of ߙ in the GWNBRg model will be the same as that obtained through non-spatial (or global) negative binomial regression. 

The parameters ‫ݑ(ࢼ‬, ‫ݒ‬ ) are estimated using the IRLS method, as in the GWNBR model described earlier, assuming that ߙ ො = ߙ ො for all i. Note that it is not necessary to alternate the NR and IRLS methods, because once ߙ is estimated globally, they estimate ‫ݑ(ࢼ‬, ‫ݒ‬ ) for each regression point i using only the IRLS method. 

Because there is no spatial variation for ߙ, its contribution to the effective number of the parameters in the model is the unity, i. e., ݇ ଶ = 1. Consequently, the bandwidth can be found using the AIC criterion ( 16 ), where Cleveland and Devlin (1988).
section_name: GWNBRG MODEL

section_content: The SAS ® macros use the IML (Interactive Matrix Language) Procedure and the parameters are described as follows. SAS ® %golden macro is used to the Golden Section Search (find the optimal bandwidth) and %gwnbr macro is used to estimated GWNBR model. 

%golden (data=,y=,x=,lat=,long=,method=,type=,gwr=,offset=,out=) ;

%gwnbr (data=,y=,x=,lat=,long=,h=,grid=,latg=,longg=,gwr=,method=,alphag=, offset=,geocod=,out=) ;

Input arguments of %golden and %gwnbr macros are defined as follows:

• data = specifies the data set to be analyzed. 

• y = specifies the response or dependent variable. 

• x = specifies the independent or explicative variables. 

• lat = specifies the latitude or y axis variable. 

• long = specifies the longitude or x axis variable. 

• method = specifies the method to be used to find the bandwidth: FIXED, ADAPTIVE1 (1 bandwidth related to the number of neighbors ) or ADAPTIVEN ( n bandwidths). 

• type = specifies the statistic to be used to find the bandwidth: AIC, CV (Cross Validation) or DEV (Deviance). 

• gwr = specifies the model: GLOBAL (GWNBRg), LOCAL (GWNBR), POISSON. 

• offset = specifies the offset variable. 

• out = specifies output data set. 

• h = specifies the bandwidth. 

• grid = specifies the data set with coordinates to be estimate the parameters. 

• latg = specifies the latitude or y axis variable of the grid data set. 

• longg = specifies the longitude or x axis variable of the grid data set. 

• alphag= specifies the global alpha value to be used in GWNBRg. 

• geocod = specifies the geocoding variable.
section_name: SAS ® MACRO

section_content: To illustrate how GWNBR works, we will use the data set presented by Silva and Rodrigues (2014) about the vehicles used for road freight transportation in Espirito Santo, Brazil. Also we use the Poisson and negative binomial models, in their global and spatial forms. The results of golden section search and of the fit are in Table 1. %gwnbr (data=data_gwnbr,y=fleet,x=industry,lat=y,long=x,h=9. 3796134, gwr=poisson,method=fixed,geocod=geocod,out=gwr) ;
section_name: ILLUSTRATION

section_content: As shown in Table 1, the PR and GWPR models were the worst models. The bandwidth found for GWPR with the best ‫,ࢼ(ܮ‬ ࢻ) is inappropriate because there are only between 1 and 9 points for each regression. The problem detected in the estimation of the bandwidth was a clue of the lack of fit obtained with the Poisson distribution. Thus, we can conclude that the fleet of vehicles used in road freight transportation probably exhibits overdispersion. Therefore, the negative binomial distribution is the best approach for the modeling of these variables. 

As an example, the surface of the parameters for GWNBR with bandwidth equal to 53. 2 km and their standard errors are shown in Figure 1. The greater values for the intercept (ܾ ) are concentrated in the Vitoria metropolitan region (southeast side), which is the capital of the state. This finding reflects the greater amount of vehicles located in that place. In contrast, the parameter estimates for ܾ ଵ are smaller in that area because there are many industries and the link function is exponential.
section_name: Model

section_content: Output 1 shows the output of %gwnbr macro, where we can see, besides descriptive statistics, the actual alpha-level, t-critical and the number of parameter estimated by GWNBR model, following Silva and Fotheringham (2015), and pseudo ܴ ଶ and adjusted ܴ ଶ measures, considering Deviance (pctdev and adjpctdev, respectively) and considering likelihood (pctll and adjpctll, respectively), following Cameron and Windmeijer (1996). print gwr method ll dev pctdev adjpctdev pctll adjpctll npar aic aicc bic; create _res_ from res[colname={"_id_" "xcoord" "ycoord" "yobs" "yhat" "res" "resraw"}]; append from res; stat=ll|| dev|| pctdev || adjpctdev|| pctll || adjpctll || npar|| aic|| aicc|| bic; create _stat_ from stat[colname={"l1" "dev" "pctdev" "adjpctdev" "pctll" "adjpctll" "npar" "aic" "aicc" "bic"}]; _colname1_={"Intercept" &x}; _label_=repeat("std_",ncol(x))//repeat("tstat_",nco l(x))//repeat("probt_",ncol(x)); _colname_={"&geocod" "x" "y"}||_colname1_||concat(_label_,repeat(_colname1_, 3))`; call change(_colname_, "_ ", "_"); call change(_colname_, "_ ", "_"); create _parameters_ from _bistdt_[colname=_colname_]; append from _bistdt_;

close _parameters_; _sig_=j(n,ncol(x),"not significant at 90%"); v1=npar; do i=1 to n;
section_name: Figure 1. Surface of the parameter estimates and standard errors obtained with the GWNBR model

section_content: nome={"id" "geocod" "x" "y" "b" "sebi" "tstat" "probtstat"}; (data=,y=,x=,lat=,long=,h=,grid=,latg=,longg=, gwr=,method=,alphag=,offset=,geocod=,rep= else if method= "adaptiven" then do;

if type="cv" then do; if _probt_[i,j]<0. 01*(ncol(x)/v1) then _sig_[i,j]="significant at 99%"; else if _probt_[i,j]<0. 05*(ncol(x)/v1) then _sig_[i,j]="significant at 95%"; else if _probt_[i,j]<0. 1*(ncol(x)/v1) then _sig_[i,j]="significant at 90%"; else _sig_[i,j]="not significant at 90%";