Code
library(tidyverse)
library(readxl)
library(plotly)
library(spdep)
library(spatialreg)
library(sphet)
library(lmtest)
library(car)
`%ni%` <- negate(`%in%`)Econometria Espacial - Resultados Preliminares
Nesse estudo estamos interessados em entender o impacto da qualidade do saneamento básico na taxa de mortalidade infantil e se existe alguma dependência espacial nessa relação.
Para termos um modelo mais bem especificado e uma maior probabilidade de encontrar o efeito real do saneamento básico na taxa de mortalidade infantil, adicionamos covariáveis relacionadas a saúde (acesso a água), educação (nota no IDEB) e economia (PIB per capita e Salário Médio Mensal). Todas as informações são referentes ao ano de 2020, com exceção do IDEB e do salário médio mensal, que são referentes ao ano anterior.
1 Importação e Limpeza dos Dados
Primeiro, importamos as bases de dados. A base que denominamos por “agua_esgoto” foi extraída do SNIS, enquanto a base “dados” foi extraída através de técnicas de web scrrapping do site IBGE Cidades.
Code
agua_esgoto <- read_excel('agua_esgoto.xlsx')
dados <- read_excel('indicadores_munic2.xlsx')No código abaixo, fazemos apenas algumas manipulações nos dados para garantir que eles estão no padrão adequado para modelagem.
Code
agua_esgoto <- agua_esgoto %>%
dplyr::select(c("Código do IBGE", "IN015_AE - Índice de coleta de esgoto", "IN055_AE - Índice de atendimento total de água")) %>%
set_names(c('codigo_munic', 'esgoto_2020', 'agua_2020')) %>%
mutate(esgoto_2020 = as.numeric(gsub(',', '.', esgoto_2020)),
agua_2020 = as.numeric(gsub(',', '.', agua_2020)))
dados <- dados %>%
pivot_wider(names_from = indicador, values_from = valor) %>%
rename(densd_demo_2010 = `Densidade demográfica [2010]`,
pessoal_ocupado_2020 = `Pessoal ocupado [2020]`,
popul_ocupada_2020 = `População ocupada [2020]`,
percnt_menor_meio_sal_2010 = `Percentual da população com rendimento nominal mensal per capita de até 1/2 salário mínimo [2010]`,
ideb_inicio_fundam_2019 = `IDEB – Anos iniciais do ensino fundamental (Rede pública) [2019]`,
ideb_fim_fundam_2019 = `IDEB – Anos finais do ensino fundamental (Rede pública) [2019]`,
percent_recet_extern_2015 = `Percentual das receitas oriundas de fontes externas [2015]`,
total_receita_2017 = `Total de receitas realizadas [2017]`,
total_despesas_2017 = `Total de despesas empenhadas [2017]`,
intern_diarreia_2016 = `Internações por diarreia [2016]`,
esgoto_adequado_2010 = `Esgotamento sanitário adequado [2010]`,
arboriz_vias_pub_2010 = `Arborização de vias públicas [2010]`,
urbaniz_vias_pub_2010 = `Urbanização de vias públicas [2010]`,
populacao_2010 = `População no último censo [2010]`,
sal_medio_mensal_2020 = `Salário médio mensal dos trabalhadores formais [2020]`,
tx_escolariz_6a14_2010 = `Taxa de escolarização de 6 a 14 anos de idade [2010]`,
pib_per_capita_2019 = `PIB per capita [2019]`,
mortal_infantil_2020 = `Mortalidade Infantil [2020]`,
area_unid_territ_2021 = `Área da unidade territorial [2021]`,
populacao_risco_2010 = `População exposta ao risco [2010]`) %>%
mutate(densd_demo_2010 = gsub('_', '.', gsub('\\.', '', gsub(',', '_', gsub(' hab/km²', '', densd_demo_2010)))),
pessoal_ocupado_2020 = gsub('\\.', '', gsub(' pessoas', '', pessoal_ocupado_2020)),
popul_ocupada_2020 = gsub(',', '.', gsub(' %', '', popul_ocupada_2020)),
percnt_menor_meio_sal_2010 = gsub(',', '.', gsub(' %', '', percnt_menor_meio_sal_2010)),
ideb_inicio_fundam_2019 = gsub(',', '.', ideb_inicio_fundam_2019),
ideb_fim_fundam_2019 = gsub(',', '.', ideb_fim_fundam_2019),
percent_recet_extern_2015 = gsub(',', '.', gsub(' %', '', percent_recet_extern_2015)),
total_receita_2017 = gsub('_', '.', gsub('\\.', '', gsub(',', '_', gsub(' R\\$ \\(×1000)', '', total_receita_2017)))),
total_despesas_2017 = gsub('_', '.', gsub('\\.', '', gsub(',', '_', gsub(' R\\$ \\(×1000)', '', total_despesas_2017)))),
intern_diarreia_2016 = gsub(',', '.', gsub(' internações por mil habitantes', '', intern_diarreia_2016)),
esgoto_adequado_2010 = gsub(',', '.', gsub(' %', '', esgoto_adequado_2010)),
arboriz_vias_pub_2010 = gsub(',', '.', gsub(' %', '', arboriz_vias_pub_2010)),
urbaniz_vias_pub_2010 = gsub(',', '.', gsub(' %', '', urbaniz_vias_pub_2010)),
populacao_2010 = gsub('\\.', '', gsub(' pessoas', '', populacao_2010)),
sal_medio_mensal_2020 = gsub(',', '.', gsub(' salários mínimos', '', sal_medio_mensal_2020)),
tx_escolariz_6a14_2010 = gsub(',', '.', gsub(' %', '', tx_escolariz_6a14_2010)),
pib_per_capita_2019 = gsub('_', '.', gsub('\\.', '', gsub(',', '_', gsub(' R\\$', '', pib_per_capita_2019)))),
mortal_infantil_2020 = gsub(',', '.', gsub(' óbitos por mil nascidos vivos', '', mortal_infantil_2020)),
area_unid_territ_2021 = gsub('_', '.', gsub('\\.', '', gsub(',', '_', gsub(' km²', '', area_unid_territ_2021)))),
populacao_risco_2010 = gsub('\\.', '', gsub(' pessoas', '', populacao_risco_2010))) %>%
mutate(across(!c(uf, nome_munic, codigo_munic), as.numeric)) %>%
left_join(agua_esgoto, by = 'codigo_munic')
dados <- dados %>%
dplyr::select(c("uf", "codigo_munic", "mortal_infantil_2020", "sal_medio_mensal_2020",
"pib_per_capita_2019", "esgoto_2020", "agua_2020", "ideb_fim_fundam_2019"))Por fim, apresentamos a estrutura dos nossos dados.
tibble [5,570 x 8] (S3: tbl_df/tbl/data.frame)
$ uf : chr [1:5570] "RO" "RO" "RO" "RO" ...
$ codigo_munic : num [1:5570] 1100015 1100379 1100403 1100346 1100023 ...
$ mortal_infantil_2020 : num [1:5570] 6.01 NA 4.59 9.52 7.3 ...
$ sal_medio_mensal_2020: num [1:5570] 1.9 1.9 1.9 1.6 2 1.9 1.8 2 1.9 2.2 ...
$ pib_per_capita_2019 : num [1:5570] 21601 21802 16729 18066 23908 ...
$ esgoto_2020 : num [1:5570] NA NA NA 95.48 2.47 ...
$ agua_2020 : num [1:5570] 47.15 23.28 6.37 60.48 84.7 ...
$ ideb_fim_fundam_2019 : num [1:5570] 5.2 4.5 4.4 4.7 5.1 5 4.6 4.4 5.3 4.4 ... uf codigo_munic mortal_infantil_2020
Length:5570 Min. :1100015 Min. : 1.20
Class :character 1st Qu.:2512126 1st Qu.: 9.09
Mode :character Median :3146280 Median : 13.33
Mean :3253591 Mean : 16.03
3rd Qu.:4119190 3rd Qu.: 19.87
Max. :5300108 Max. :142.86
NA's :1570
sal_medio_mensal_2020 pib_per_capita_2019 esgoto_2020 agua_2020
Min. :0.900 Min. : 4483 Min. : 0.00 Min. : 0.00
1st Qu.:1.700 1st Qu.: 10516 1st Qu.: 36.88 1st Qu.: 52.87
Median :1.900 Median : 18182 Median : 71.31 Median : 75.25
Mean :1.985 Mean : 24547 Mean : 63.70 Mean : 70.62
3rd Qu.:2.200 3rd Qu.: 29883 3rd Qu.: 94.98 3rd Qu.: 93.48
Max. :6.000 Max. :464884 Max. :100.00 Max. :100.00
NA's :2774 NA's :233
ideb_fim_fundam_2019
Min. :1.900
1st Qu.:4.100
Median :4.700
Mean :4.613
3rd Qu.:5.100
Max. :7.800
NA's :283 2 Qualidade e disponibilidade dos dados
Primeiro, apresentamos a proporção dos municípios em cada estado que possuem todos os dados necessários para nossa análise.
Code
dados_na <- dados %>%
mutate(num_na = as.numeric(rowSums(!is.na(.)) == ncol(.))) %>%
group_by(uf) %>%
summarise(prop_na = sum(num_na) / length(num_na))
fig <- plot_ly(dados_na, x = ~reorder(uf, -prop_na), y = ~prop_na, type = 'bar',
marker = list(color = '#3885b5',
line = list(color = '#3885b5',
width = 1.5)))
fig %>% layout(title = list(text = 'Proporção de dados mompletos', y = 0.97),
yaxis = list(title = 'Proporção', tickformat = ".0%"),
xaxis = list(title = ''))Pelo gráfico acima, podemos ver quão severa é nossa perda de dados para garantir que temos informações de todos os indicadores. Ademais, podemos perceber que essa perda de dados não é aleatória. De fato, observamos uma densidade baixa de municípios com todos os dados no norte, principalmente, mas em outras regiões também, como Nordeste.
Entretanto, o Sudeste parece apresentar uma maior concentração de municípios com todos os dados disponíveis. Por esse motivo, iremos focar nossa análise nessa região. Fazemos isso ao custo de uma generalização dos resultados do nosso modelo, mas, concomitantemente, garantimos uma maior validade interna.
Code
dados <- dados %>%
dplyr::filter(uf %in% c('MG', 'SP', 'RJ', 'ES')) %>%
drop_na()3 Análise Econométrica
3.1 Shapefile e Matriz de Pesos
Primeiro importamos o shapefile do Brasil, extraído do IBGE, e garantimos que temos, na mesma ordem, os mesmos municípios em ambas as bases de dados.
Code
db <- st_read('BR_Municipios_2021.shp')Reading layer `BR_Municipios_2021' from data source
`C:\Users\Pedro\Desktop\UFMG20221\Econometria Espacial\saude\BR_Municipios_2021.shp'
using driver `ESRI Shapefile'
Simple feature collection with 5572 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -73.99045 ymin: -33.75118 xmax: -28.84764 ymax: 5.271841
Geodetic CRS: SIRGAS 2000Code
db <- db %>%
dplyr::filter(CD_MUN %in% dados$codigo_munic) CD_MUN NM_MUN SIGLA AREA_KM2
Length:1029 Length:1029 Length:1029 Min. : 3.565
Class :character Class :character Class :character 1st Qu.: 213.881
Mode :character Mode :character Mode :character Median : 407.427
Mean : 644.934
3rd Qu.: 737.986
Max. :10727.097
geometry
MULTIPOLYGON :1029
epsg:4674 : 0
+proj=long...: 0
Agora apresentamos uma representação visual dos nossos dados. Como podemos observar, após restringirmos nossa amostra para os municípios do sudeste, temos uma densidade bem maior de dados e os municípios sem dados estão menos clusterizados. Ou seja, nessa subamostra, a hipótese nula de aleatoriedade dos valores faltantes parece ser mais difícil de ser rejeitada que no caso do Brasil como um todo.
Code
plot(st_geometry(db))
Nesse momento, construímos a matriz de pesos. Como ainda temos uma presença relevante de valores em branco, utilizamos uma matriz de pesos baseada no inverso da distância.
Code
coords <- st_centroid(st_geometry(db))
# Garantimos que todos os municipios tenham pelo menos um vizinho
k1 <- knn2nb(knearneigh(coords))
critical.threshold <- max(unlist(nbdists(k1,coords)))
nb.dist.band <- dnearneigh(coords, 0, critical.threshold+1)
distances <- nbdists(nb.dist.band,coords)
# https://spatialanalysis.github.io/lab_tutorials/Spatial_Weights_as_Distance_Functions.html#inverse-distance-weights
invd1 <- lapply(distances, function(x) (1/(x/100)))
W <- nb2listw(nb.dist.band, glist = invd1, style = "W")Characteristics of weights list object:
Neighbour list object:
Number of regions: 1029
Number of nonzero links: 73136
Percentage nonzero weights: 6.907175
Average number of links: 71.07483
Link number distribution:
1 3 4 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
1 2 1 1 2 1 3 4 4 2 3 2 6 7 9 8 6 1 6 10
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
6 10 5 2 4 6 7 7 13 1 8 5 14 6 4 9 11 6 8 14
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
14 14 10 11 11 7 15 3 11 6 5 12 3 11 9 7 9 10 10 5
63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82
8 10 7 15 12 6 13 11 14 17 10 21 14 19 23 17 16 14 19 14
83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
16 14 14 12 10 12 9 12 14 9 7 9 7 6 5 8 3 7 5 2
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 122 123
3 5 4 3 1 4 5 1 4 2 4 2 7 4 1 1 5 1 3 3
124 125 126 127 129 131 132 133 134 135 136 137 138 140 141 142 143 144 145 146
3 6 6 2 2 3 6 3 3 1 3 3 3 5 2 2 5 4 3 5
147 148 149 150 151 152 153 154 157
3 5 5 2 3 2 3 2 2
1 least connected region:
523 with 1 link
2 most connected regions:
836 884 with 157 links
Weights style: W
Weights constants summary:
n nn S0 S1 S2
W 1029 1058841 1029 59.13711 4136.9443.2 Testes estatísticos e regressões
Começamos primeiro com uma regressão linear simples e, em seguida, testes de heterocedasticidade e multicolinearidade.
Na regressão abaixo, o intercepto, salário médio mensal, esgoto (10%) e água foram significativos. O sinal de todos esses coeficientes está dentro do esperado, com exceção do relacionado à cobertura de esgoto.
Code
ols <- lm(mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019, dados)
summary(ols)
Call:
lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 + pib_per_capita_2019 +
esgoto_2020 + agua_2020 + ideb_fim_fundam_2019, data = dados)
Residuals:
Min 1Q Median 3Q Max
-18.095 -5.907 -2.368 2.598 129.526
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.374e+01 3.038e+00 7.813 1.38e-14 ***
sal_medio_mensal_2020 -2.383e+00 8.134e-01 -2.929 0.00348 **
pib_per_capita_2019 -1.615e-05 1.090e-05 -1.481 0.13884
esgoto_2020 3.332e-02 1.420e-02 2.346 0.01915 *
agua_2020 -8.291e-02 1.859e-02 -4.459 9.16e-06 ***
ideb_fim_fundam_2019 8.005e-02 6.285e-01 0.127 0.89867
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.46 on 1023 degrees of freedom
Multiple R-squared: 0.06078, Adjusted R-squared: 0.05619
F-statistic: 13.24 on 5 and 1023 DF, p-value: 1.608e-12Pelo teste de Breusch-Pagan, apresentado abaixo, podemos perceber que temos heterocedasticidade. Entretanto, esse é um problema menor, dado que ele reduz a eficiência do estimador, mas não o torna inconsistente.
O teste VIF, por sua vez, indica que não temos multicolinearidade.
[1] "BP TEST:"
studentized Breusch-Pagan test
data: ols
BP = 12.159, df = 5, p-value = 0.03268[1] "VIF:"sal_medio_mensal_2020 pib_per_capita_2019 esgoto_2020
1.632346 1.405139 1.117357
agua_2020 ideb_fim_fundam_2019
1.312591 1.155142 Agora, entrando nos testes espaciais, pelo I de Moran Global para resíduos de regressão, rejeitamos a hipótese nula de aleatoriedade espacial. Assim, podemos concluir que existe uma dependência espacial nos dados.
Code
lm.morantest(ols, W)
Global Moran I for regression residuals
data:
model: lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019,
data = dados)
weights: W
Moran I statistic standard deviate = 2.1123, p-value = 0.01733
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
1.381570e-02 -1.494265e-03 5.253205e-05 Agora, através dos testes diagnósticos de multiplicador de Lagrange para dependência espacial, consideramos que a melhor forma funcional para modelar essa dependência é o modelo de erro espacial.
Code
lm.LMtests(ols, W, test = 'all')
Lagrange multiplier diagnostics for spatial dependence
data:
model: lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019,
data = dados)
weights: W
LMerr = 3.4176, df = 1, p-value = 0.06451
Lagrange multiplier diagnostics for spatial dependence
data:
model: lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019,
data = dados)
weights: W
LMlag = 1.5622, df = 1, p-value = 0.2113
Lagrange multiplier diagnostics for spatial dependence
data:
model: lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019,
data = dados)
weights: W
RLMerr = 3.3286, df = 1, p-value = 0.06809
Lagrange multiplier diagnostics for spatial dependence
data:
model: lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019,
data = dados)
weights: W
RLMlag = 1.4732, df = 1, p-value = 0.2248
Lagrange multiplier diagnostics for spatial dependence
data:
model: lm(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020 + ideb_fim_fundam_2019,
data = dados)
weights: W
SARMA = 4.8908, df = 2, p-value = 0.08669Antes de partimos para a modelagem em si, primeiro temos que verificar se nossa variável dependente segue uma normal para podermos utilizar modelos baseados em Máxima Verossimilhança.
Como podemos ver pelo teste de Shapiro, a hipótese nula de normalidade é rejeitada. Desse modo, não podemos utilizar modelos baseados em Máxima Verossimilhança.
Code
shapiro.test(dados$mortal_infantil_2020)
Shapiro-Wilk normality test
data: dados$mortal_infantil_2020
W = 0.73427, p-value < 2.2e-16Dado que não podemos utilizar modelos baseados em ML, partimos direto para o modelo de erro espacial por GMM.
Code
errormodel <- spreg(mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020,
dados, W, model = 'error', het = TRUE)
Generalized stsls
Call:
spreg(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020, data = dados,
listw = W, model = "error", het = TRUE)
Residuals:
Min. 1st Qu. Median Mean 3rd Qu. Max.
-18.308 -5.899 -2.295 0.005 2.593 129.713
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 2.4387e+01 2.0000e+00 12.1937 < 2.2e-16 ***
sal_medio_mensal_2020 -2.3224e+00 6.9784e-01 -3.3280 0.0008747 ***
pib_per_capita_2019 -1.6119e-05 7.6166e-06 -2.1163 0.0343201 *
esgoto_2020 3.2787e-02 1.5130e-02 2.1670 0.0302370 *
agua_2020 -8.7244e-02 2.0915e-02 -4.1714 3.027e-05 ***
rho 1.8073e-01 9.9610e-02 1.8144 0.0696232 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Entretanto, como a omissão da variável Lag no modelo pode gerar um viés de variável omitida, estimamos esse modelo mesmo os testes tendo indicado que não existe dependência espacial do tipo Lag.
Como podemos ver no resultado abaixo, a variável Lag foi significativa. Por esse motivo e pelos expostos acima, consideramos esse como sendo o modelo que melhor descreve os dados.
Code
lagmodel <- spreg(mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020,
dados, W, model = 'lag', het = TRUE)
Stsls
Call:
spreg(formula = mortal_infantil_2020 ~ sal_medio_mensal_2020 +
pib_per_capita_2019 + esgoto_2020 + agua_2020, data = dados,
listw = W, model = "lag", het = TRUE)
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-17.4266 -6.1373 -2.1779 3.1221 129.5482
Coefficients:
Estimate Std. Error t-value Pr(>|t|)
(Intercept) 1.2488e+01 6.1082e+00 2.0444 0.0409162 *
sal_medio_mensal_2020 -1.6428e+00 8.3465e-01 -1.9683 0.0490365 *
pib_per_capita_2019 -1.7592e-05 7.9377e-06 -2.2163 0.0266712 *
esgoto_2020 3.4844e-02 1.5196e-02 2.2929 0.0218527 *
agua_2020 -7.8572e-02 2.0555e-02 -3.8225 0.0001321 ***
lambda 6.6255e-01 3.1271e-01 2.1187 0.0341128 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1