Continuando, mas ainda não finalizando, a reprodução do artigo baseado em minha dissertação de mestrado. Essa última regressão utilizada no artigo não fez parte da dissertação que o originou. Ela foi sugerida por um dos revisores. Na verdade, ela já havia sido sugerida quando apresentei uma versão desse artigo no congresso da ANPEC, mas, na época, acabei não fazendo.

Pacotes necessários

Para a regressão em painel utilizei dois pacotes adicionais ao spedp. O plm, que é para dados em painel e o splm, que é para dados em painel espacializados (spatial). Um pacote adicional que não havia utilizado, mas que o próprio R solicitou, foi o spatialreg.

library(spdep)
library(plm)
library(splm)
library(spatialreg)

Base de Dados

A base de dados se mantém a mesma, bem como a matriz de vizinhança utilizada.

#Leitura da Base de Dados
BASE <- read.csv("http://raw.githubusercontent.com/RodrigoAnderle/Artigos-Reproduz-veis/master/Fatores%20Espaciais%20Comuns%20(RBERU)/NE1991.csv"
                 , header = T)

#Leitura da matriz de vizinhança
R1<-read.gal("https://raw.githubusercontent.com/RodrigoAnderle/Artigos-Reproduz-veis/master/Fatores%20Espaciais%20Comuns%20(RBERU)/R1.gal"
             ,BASE$MUN) 
WR1<-nb2listw(R1) # ajuste necessário para utilização nas funções

Transformando em Painel

Para transformar os dados em painel, é necessário inserir o ano de referência, o que não constava na minha base anterior e não fazia parte do plano inicial da pesquisa. Assim, tive que readequar a base para conseguir fazer uma análise em painel.

O primeiro procedimento separa as variáveis em que o ano base é 1991. Em seguida, é adicionada a variável ANO e ajustados os nomes das variáveis. O segundo procedimento faz a mesma coisa para as variáveis com ano base 2000.

Os últimos procedimentos são para unir novamente as bases numa só, com seus respectivos anos bases. Note que o ajuste no nome das variáveis foi necessário para que essa junção fosse possível. Por fim, a BASE é transformada em uma pdata.frame, uma base em painel, em que os índices são os Municípios (MUN) e os anos (ANO).

SB1991 <-BASE[,c(2,4,6,9,11,13,15,17,19,21,23,31:39)]
  SB1991$ANO <- 1991
  colnames(SB1991)<- c("MUN","TAXA","PRODUT","ESCOL","EMED","URB",
                     "ANALF","GINI","LNPOP","INDUS","SERV",
                     "MA","PI","CE","RN","PB","PE","AL","SE","BA",
                     "ANO")

SB2000 <- BASE[,c(2,5,7,10,12,14,16,18,20,22,24,31:39)]
  SB2000$ANO <- 2000
  colnames(SB2000)<- c("MUN","TAXA","PRODUT","ESCOL","EMED","URB",
                     "ANALF","GINI","LNPOP","INDUS","SERV",
                     "MA","PI","CE","RN","PB","PE","AL","SE","BA",
                     "ANO")

BASE <- rbind(SB1991,SB2000)
BASE <- pdata.frame(BASE, index = c("MUN", "ANO"))

Painel com efeito espacial

Seriam duas as regressões possíveis. A primeira de efeitos fixos (fixed effects), ou seja, são fixados os municípios, assumindo que cada município tem um impacto determinado ao longo do tempo. E, a outra, de efeitos aleatórios (random effects), assumindo que eles não têm esse impacto. É possível testar qual modelo é mais ajustado, mas, a partir da discussão do artigo, claramente, a hipótese de efeitos fixos é a mais adequada para o problema de interesse. De toda forma, em meus rascunhos, fiz as duas e testei a de maior ajusteque é a que foi utilizada e é a que reproduzo abaixo:

#Efeitos fixos
plm_PIB<- spml(TAXA ~ PRODUT + EMED + ANALF + URB + LNPOP + INDUS +
                 SERV + GINI,data = BASE, listw = WR1, listw2 = WR1,
               model = "within",lag = T, effect = "twoways",
               spatial.error = "kkp")
summary(plm_PIB)
## Spatial panel fixed effects sarar model
##  
## 
## Call:
## spml(formula = TAXA ~ PRODUT + EMED + ANALF + URB + LNPOP + INDUS + 
##     SERV + GINI, data = BASE, listw = WR1, listw2 = WR1, model = "within", 
##     effect = "twoways", lag = T, spatial.error = "kkp")
## 
## Residuals:
##        Min.     1st Qu.      Median     3rd Qu.        Max. 
## -5.0526e-01 -1.5976e-02 -5.2042e-18  1.5976e-02  5.0526e-01 
## 
## Spatial error parameter:
##      Estimate Std. Error t-value  Pr(>|t|)    
## rho -0.619726   0.039211 -15.805 < 2.2e-16 ***
## 
## Spatial autoregressive coefficient:
##        Estimate Std. Error t-value  Pr(>|t|)    
## lambda 0.816149   0.013182  61.913 < 2.2e-16 ***
## 
## Coefficients:
##          Estimate Std. Error  t-value  Pr(>|t|)    
## PRODUT -0.1147395  0.0066562 -17.2381 < 2.2e-16 ***
## EMED    0.1097097  0.0302151   3.6310 0.0002824 ***
## ANALF  -0.0679559  0.0169690  -4.0047 6.210e-05 ***
## URB     0.0059893  0.0072540   0.8256 0.4090039    
## LNPOP   0.0087238  0.0015611   5.5884 2.292e-08 ***
## INDUS   0.0440414  0.0171928   2.5616 0.0104184 *  
## SERV    0.1142844  0.0139220   8.2089 2.232e-16 ***
## GINI   -0.1579501  0.0490972  -3.2171 0.0012950 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note que, nos resultados, o sentido do rho e do lambda estão invertidos em relação a notação que eu utilizo no artigo e nos outros textos sobre as regressões. Assim, neste modelo o rho é o Spatial error parameter e o lambda, o Spatial autoregressive coefficient.

Modelo Durbin Spatial em painel

O modelo Durbin espacial não tem uma função pre-determinada no pacote splm. Assim, foi necessário criar as defasagens espaciais manualmente para, em seguida, fazer a regressão. Estes são os procedimentos descritos a seguir.

Defasagem das explicativas

Para criar a defasagem espacial foi utilizada a função slag que cria um lag, uma defasagem da variável em questão. Isso é mais comum em modelos de séries temporais, em que esse lag é dado no tempo. No nosso caso, esse lag é “espacial”, dos vizinhos, com base em nossa matriz de vizinhança WR1.

### Inserindo defasagem espacial nas explicativas
BASE$W_PRODUT <- slag(BASE$PRODUT,WR1)
BASE$W_ESCOL <- slag(BASE$ESCOL,WR1)
BASE$W_EMED <- slag(BASE$EMED,WR1)
BASE$W_URB <- slag(BASE$URB,WR1)
BASE$W_ANALF <- slag(BASE$ANALF,WR1)
BASE$W_GINI <- slag(BASE$GINI,WR1)
BASE$W_LNPOP <- slag(BASE$LNPOP,WR1)
BASE$W_INDUS <- slag(BASE$INDUS,WR1)
BASE$W_SERV <- slag(BASE$SERV,WR1)

Modelo Durbin Spatial em painel

A regressão em si é análoga a anterior do painel espacializado com efeito fixo.A exceção são os lags da vizinhança.

#Efeitos Fixos
Bplm_PIB<- spml(TAXA ~ PRODUT + EMED + ANALF + URB + LNPOP + INDUS +
                 SERV + GINI + W_PRODUT + W_EMED + W_ANALF + W_URB +
                  W_LNPOP + W_INDUS + W_SERV + W_GINI + W_PRODUT
                ,data = BASE, listw = WR1, listw2 = WR1,
               model = "within",lag = T, effect = "twoways",
               spatial.error = "kkp")

summary(Bplm_PIB)
## Spatial panel fixed effects sarar model
##  
## 
## Call:
## spml(formula = TAXA ~ PRODUT + EMED + ANALF + URB + LNPOP + INDUS + 
##     SERV + GINI + W_PRODUT + W_EMED + W_ANALF + W_URB + W_LNPOP + 
##     W_INDUS + W_SERV + W_GINI + W_PRODUT, data = BASE, listw = WR1, 
##     listw2 = WR1, model = "within", effect = "twoways", lag = T, 
##     spatial.error = "kkp")
## 
## Residuals:
##        Min.     1st Qu.      Median     3rd Qu.        Max. 
## -4.9256e-01 -1.6030e-02  1.0408e-16  1.6030e-02  4.9256e-01 
## 
## Spatial error parameter:
##      Estimate Std. Error t-value  Pr(>|t|)    
## rho -0.608342   0.039013 -15.593 < 2.2e-16 ***
## 
## Spatial autoregressive coefficient:
##        Estimate Std. Error t-value  Pr(>|t|)    
## lambda 0.807225   0.013152  61.377 < 2.2e-16 ***
## 
## Coefficients:
##            Estimate Std. Error  t-value  Pr(>|t|)    
## PRODUT   -0.1436938  0.0077005 -18.6603 < 2.2e-16 ***
## EMED      0.1350571  0.0333658   4.0478 5.171e-05 ***
## ANALF    -0.0155788  0.0253487  -0.6146  0.538831    
## URB       0.0240038  0.0088884   2.7006  0.006922 ** 
## LNPOP     0.0096959  0.0018156   5.3404 9.274e-08 ***
## INDUS     0.0586002  0.0204855   2.8606  0.004229 ** 
## SERV      0.0787714  0.0182440   4.3177 1.577e-05 ***
## GINI     -0.0140869  0.0607553  -0.2319  0.816644    
## W_PRODUT  0.0717246  0.0096072   7.4657 8.285e-14 ***
## W_EMED   -0.0329885  0.0588687  -0.5604  0.575225    
## W_ANALF  -0.0478223  0.0348084  -1.3739  0.169481    
## W_URB    -0.0199880  0.0140558  -1.4220  0.155013    
## W_LNPOP  -0.0043202  0.0024792  -1.7426  0.081406 .  
## W_INDUS  -0.0431155  0.0322902  -1.3352  0.181795    
## W_SERV    0.0274103  0.0266508   1.0285  0.303715    
## W_GINI   -0.2564099  0.0922404  -2.7798  0.005439 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Considerações

Mais uma vez faço a ressalva que não pretendo discutir os resultados neste texto. A discussão já se encontra no artigo. Assim, este texto se propõe a disponibilizar a programação utilizada e a base de dados, permitindo a reprodução.

Para além destes textos que fiz (ver links úteis), tenho mais um que pretendo fazer que trata da apresentação dos resultados. Na dissertação e no artigo, a escolha foi por utilizar, majoritariamente, tabelas. Como pôde ser visto nos resultados acima, isso não torna muito claro os objetivos e os resultados do artigo. Por isso que pretendo fazer outro texto tratando somente da apresentação desses resultados.

Observação: fiz uma correção no dia 21/06/2019, pois a transformação da base de dados para o formato em painel estava errada. Também adicionei a explicação de que, nestes modelos da função spml , os termos rho e lambda tem sentido inverso ao que eu utilizo no artigo, apesar deles especificarem que o rho é o “Spatial error parameter” e o lambda o “Spatial autoregressive coefficient”, isto pode confundir em uma leitura rápida - ainda mais quando você não fez nenhuma anotação a respeito :(.