Análise estatística espacial no R (parte 2)

Semana passada mostramos como estimar a dependência espacial de variáveis no R. Hoje, iremos continuar a falar sobre econometria espacial, apresentando como rodar modelos de regressão que levem em consideração o espaço. Para isso, utilizaremos um exemplo simples em que estimamos a relação  da taxa de homicídios com o PIB per capita.

O argumento mais importante para se estimar modelos de econometria espacial é se a suposição de independência entre as observações não é mais válida. Assim, características de um local podem influenciar nos resultados de outro.

library(geobr)
library(sidrar)
library(spdep)

library(stargazer)

set.ZeroPolicyOption(TRUE)

Utilizaremos basicamente o mesmo processo de extração de dados do post anterior, com a utilização de microrregiões como unidade geográfica. A diferença será a inclusão da variável de com os dados de taxa de homicídio.

#Dados de PIB
pib = get_sidra(5938,
variable = 37, 
geo = "MicroRegion")
#Dados de população
pop= get_sidra(202,
variable = 93,
geo = "MicroRegion") %>%
filter(`Sexo (Código)` == 0 & `Situação do domicílio (Código)` == 0)

#Juntando e criando a variável de pib per capita
df <- left_join(pib, pop, by = "Microrregião Geográfica (Código)") %>%
mutate(pibpc = (Valor.x/Valor.y) - mean(Valor.x/Valor.y)) #centralizar no 0

# colocando a taxa de homicídios
tx_homic <- ipeadatar::ipeadata("THOMIC") %>%
filter(uname == "Microregions" &
date == "2018-01-01") %>%
mutate("Microrregião Geográfica (Código)" = as.character(tcode)) %>%
rename(tx_homic = value) %>%
select(tx_homic, "Microrregião Geográfica (Código)" )

df <- left_join(df, tx_homic, by = "Microrregião Geográfica (Código)")

# mapa
mapa_micro = read_micro_region(showProgress = F)
mapa_micro$code_micro = as.character(mapa_micro$code_micro)

merged = dplyr::left_join(mapa_micro, df,
by = c("code_micro" = "Microrregião Geográfica (Código)"))

Também iremos produzir a matriz de vizinhança da mesma forma que antes. Iremos retirar da análise a microrregião de Fernando de Noronha, pois ela não apresenta vizinhos, o que impossibilita a estimação.

nb <- poly2nb(merged, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
df = df  %>%  filter(`Microrregião Geográfica.x` != "Fernando de Noronha - PE")

Existem vários modelos que se encaixam nesse arcabouço, mas os dois principais são o Autorregressivo Espacial (SAR) e o de Erros Espaciais (SEM). Ambos são modificações do modelo clássico OLS:

(1)  

O modelo SAR introduz um processo autoregressivo, semelhante ao que ocorre nos modelos de séries temporais:

(2)  

Já o modelo SEM introduz o componente do processo autoregressivo no termo de erro:

(3)  

Para determinar qual dos dois modelos rodar, é preciso estimar o teste do tipo multiplicador de Lagrange para avaliar se o Lambda e o W são diferentes de zero.

ols <- lm(tx_homic ~ pibpc, data = df)

lmLMtests <- lm.LMtests(ols,lw, test=c("RLMerr", "RLMlag"))
lmLMtests

Como em ambos os casos a hipótese nula foi rejeitada, o procedimento padrão é estimar o modelo que tenha retornado a estatística mais significativa. No nosso caso, estimaremos o modelo SAR.

sar = lagsarlm(tx_homic ~ pibpc, data = df, lw, zero.policy = TRUE)
stargazer(ols,sar)

Veja que na especificação OLS temos que a há uma relação significativa entre as duas variáveis. Quando há a inclusão do componente espacial, o parâmetro do PIB per capita deixa de ser significativo. Isso mostra que, nesse caso, estimar o modelo OLS pode gerar o chamado viés da variável omitida.