Econometria Espacial
Introdução
Usamos econometria espacial quando os fenômenos estudados apresentam dependência espacial, isto é, quando o resultado observado em uma unidade geográfica está correlacionado com os resultados de unidades próximas. Essa dependência pode surgir tanto por interações diretas entre unidades vizinhas quanto pela presença de fatores não observados espacialmente correlacionados, gerando transbordamentos no espaço.
Diversos fenômenos podem exibir esse tipo de estrutura. Em eleições, por exemplo, ao analisarmos dados em nível municipal ou estadual, observamos que unidades geograficamente próximas tendem a apresentar padrões de votação semelhantes, refletindo afinidades ideológicas regionais. Outro exemplo de fenômeno com essa característica é a criminalidade, que costuma apresentar agrupamentos espaciais, com bairros ou cidades violentas frequentemente localizados próximos a outras unidades também violentas.
Em suma, segue resumidamente o que precisamos para realizar uma análise espacial
1- Arquivo shapefile das regiões (município, estado, setor censitário). Para o Brasil, é muito fácil conseguirmos, através do pacote geobr, mas nem todos países possuem pacotes equivalentes. Nesses casos, é necessário encontrar manualmente os arquivos.
2- Construir a matriz de pesos espaciais. Geralmente elas são de dois tipos, contiguidade (vizinhos que compartilham fronteira), podendo ser rook ou queen; o outro tipo de matriz de vizinhança é baseado na distância, podendo ser os k-vizinhos mais próximos ou então os vizinhos dentro de um determinado raio.
3- Antes de partirmos para a modelagem espacial, é necessário verificar se realmente é necessário. Para tanto, utiliza-se o I de Moran. Caso o valor da estatística do I de Moran seja maior que o valor esperado sob aleatoriedade, isso indica autocorrelação espacial positiva, com predomínio de padrões como Alto-Alto e Baixo-Baixo entre regiões vizinhas. Caso o valor seja menor que o esperado, há indícios de autocorrelação espacial negativa, sugerindo que regiões com do tipo Alto-Baixo e/ou Baixo-Alto. Por fim, se o valor do I de Moran for aproximadamente igual ao esperado, assume-se que há aleatoriedade espacial.
4- Identificar como os clusters se comportam localmente: I de Moran local (LISA)
5- OLS simples e teste de I de Moran nos resíduos. Caso o I de Moran dos resíduos não seja significativo, isso sugere que a especificação do modelo foi suficiente para capturar a dependência espacial presente nos dados e não há necessidade de continuar com a análise espacial.
Tomemos o seguinte modelo linear geral:
\[y=X\beta+\epsilon\]
Onde y é um vetor ×1 da variável dependente, \(X\) é uma matriz x, \(\beta\) é um vetor de parâmetros a serem estimados e \(\epsilon\) é independente e identicamente distribuído.
Em termos gerais, existem 3 tipos de correlação espacial. O primeiro delas é a autocorrelação espacial, onde a variável explicada depende do valor da variável nos vizinhos. Nesse contexto, o modelo SAR ( possui a variável dependente defasada espacialmente:
\[\begin{equation} y=\rho Wy + X\beta + \epsilon \end{equation}\]
Onde \(\rho\) é um parâmetro espacial defasado e \(W\) é a matriz x de pesos espaciais.
O segundo tipo de correlação espacial está nos erros do modelo, de tal forma que estes são autocorrelacionados espacialmente. O modelo SEM () incorpora essas defasagens espaciais do erro:
\[\begin{equation} y=X\beta+\xi \end{equation}\] \[\xi=\lambda W\xi+\epsilon\]
Por fim, podemos ter a presença de correlação espacial nas variáveis explicativas. Nesse sentido, o SDM () considera a defasagem da variável dependente e das variáveis explicativas:
\[\begin{equation} y=\rho Wy + X\beta + WZ\theta + \epsilon \end{equation}\]
Onde \(WZ\) representa as variáveis explicativas defasadas espacialmente. Ressalta-se que \(Z\) não precisa ser necessariamente igual a \(X\), apesar de ser mais comum.
Essas três formas de correlação espacial não são excludentes e podemos generalizar ainda mais, em um modelo espacial geral, dado por:
\[\begin{equation} y=\rho Wy+X\beta + WZ\theta+\xi \end{equation}\] \[\xi=\lambda W\xi + \epsilon\]
Podemos compilar as combinações na seguinte tabela:
| rho | theta | lambda | Modelo |
|---|---|---|---|
| ≠ 0 | = 0 | = 0 | SAR |
| ≠ 0 | ≠ 0 | = 0 | SDM |
| = 0 | = 0 | ≠ 0 | SEM |
| = 0 | ≠ 0 | ≠ 0 | SDEM |
| ≠ 0 | = 0 | ≠ 0 | KP |
Aplicação no R
Para aplicação, vamos fazer um exercício bem simples: regredir o percentual de votos de um partido em uma eleição em São Paulo pelo PIB per capita e Gini. As unidades geográficas serão os municípios.
Pacotes e dados
pacotes <- c("dplyr", "spdep", "sphet", "spatialreg", "tidyr", "ggplot2",
"readxl", "geobr", "lmtest", "scales", "tinytex", "purrr", "stringr", "magrittr", "spatialprobit", "Matrix", "stringr", "stringi", "viridis")
lapply(pacotes, library, character.only = TRUE)Carregando os dados e shapefile
## Using year/date 2022muni<-muni %>%
rename("municipio"="name_muni") %>%
mutate(municipio = toupper(municipio))
#municípios e votos
df<-left_join(muni, df, by ="municipio")Após organizar a base e colocar as variáveis explicativas:
## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -51.17924 ymin: -23.03648 xmax: -46.54881 ymax: -21.19703
## Geodetic CRS: SIRGAS 2000
## municipio code_state abbrev_state name_state code_region
## 1 ADAMANTINA 35 SP São Paulo 3
## 2 ADOLFO 35 SP São Paulo 3
## 3 AGUAÍ 35 SP São Paulo 3
## 4 ÁGUAS DA PRATA 35 SP São Paulo 3
## 5 ÁGUAS DE LINDÓIA 35 SP São Paulo 3
## 6 ÁGUAS DE SANTA BÁRBARA 35 SP São Paulo 3
## name_region percentual Código.município pib_pc gini_2010
## 1 Sudeste 23.78791 61018 40475.54 0.5151
## 2 Sudeste 17.75668 61034 33611.75 0.4111
## 3 Sudeste 24.74849 61050 36637.52 0.4804
## 4 Sudeste 26.31063 61077 26771.36 0.5305
## 5 Sudeste 26.97511 61093 31760.71 0.4612
## 6 Sudeste 29.57703 70190 36477.35 0.4948
## geom
## 1 MULTIPOLYGON (((-51.09171 -...
## 2 MULTIPOLYGON (((-49.61226 -...
## 3 MULTIPOLYGON (((-47.01219 -...
## 4 MULTIPOLYGON (((-46.71728 -...
## 5 MULTIPOLYGON (((-46.61212 -...
## 6 MULTIPOLYGON (((-49.29262 -...Visualização dos dados:
ggplot(df) +
geom_sf(aes(fill = percentual), color = "grey30", linewidth = 0.1) +
scale_fill_viridis_c(option = "plasma", direction = -1)+
labs(
title = "Distribuição espacial do percentual de votos"
) +
theme_minimal(base_size = 12) +
theme(
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(face = "bold", size = 14),
legend.position = "right"
)
Teste do I de Moran
#Definindo uma matriz de pesos espaciais do tipo torre
sp_rook <- poly2nb(df$geom, queen = FALSE)
rook_listw <- nb2listw(sp_rook, zero.policy = TRUE)
# Teste do I de Moran Global
moran.test(x = df$percentual, listw = rook_listw)##
## Moran I test under randomisation
##
## data: df$percentual
## weights: rook_listw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 12.116, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.3024022617 -0.0015552100 0.0006293585Como a hipótese nula do I de Moran é ausência de autocorrelação espacial e obtivemos um p-valor muito baixo, rejeitamos a hipótese nula. Além disso, o valor da estatística é maior que o valor esperado sob a hipótese nula, portanto, esperamos observar forte presença de clusters do tipo Alto-Alto e Baixo-Baixo.
LISA
Dada a presença de autocorrelação espacial, precisamos ver como esses clusters se distribuem no espaço.
# Ilha Bela é uma Ilha, precisamos remover ela para calcular o LISA
df_lisa <- df[!is.na(lag.listw(rook_listw, scale(df$percentual))), ]
# Calcular localmoran, votos padronizados e lag
locm <- localmoran(df_lisa$percentual, listw = rook_listw)
df_lisa$votos <- scale(df_lisa$percentual) %>% as.vector()
df_lisa$lag_votos <- lag.listw(rook_listw, df_lisa$votos)
# Criar variável de cluster
df_lisa$quad_sig <- NA
df_lisa[which(df_lisa$votos >= 0 &
df_lisa$lag_votos >= 0 &
!is.na(locm[, 5]) & locm[, 5] <= 0.05), "quad_sig"] <- 1
df_lisa[which(df_lisa$votos <= 0 &
df_lisa$lag_votos <= 0 &
!is.na(locm[, 5]) & locm[, 5] <= 0.05), "quad_sig"] <- 2
df_lisa[which(df_lisa$votos >= 0 &
df_lisa$lag_votos <= 0 &
!is.na(locm[, 5]) & locm[, 5] <= 0.05), "quad_sig"] <- 3
df_lisa[which(df_lisa$votos <= 0 &
df_lisa$lag_votos >= 0 &
!is.na(locm[, 5]) & locm[, 5] <= 0.05), "quad_sig"] <- 4
df_lisa[which(!is.na(locm[, 5]) & locm[, 5] > 0.05), "quad_sig"] <- 5
# Voltar a variável para a base original
df$quad_sig <- NA
df$quad_sig[match(df_lisa$`Código.município`, df$`Código.município`)] <- df_lisa$quad_sig
# Mapear
breaks <- seq(1, 5, 1)
labels <- c("Alto-Alto", "Baixo-Baixo", "Alto-Baixo", "Baixo-Alto", "Não Significante")
np <- findInterval(df$quad_sig, breaks)
colors <- c("darkred", "darkblue", "#fcbba1", "#9ecae1", "white")
par(mar = c(0, 0, 2, 0))
plot(df$geom,
col = colors[np],
border = "grey60",
lwd = 0.4)
title(main = "I de Moran Local (LISA)",
cex.main = 1.4,
font.main = 2)
legend("topright",
legend = c("Alto–Alto", "Baixo–Baixo", "Alto–Baixo", "Baixo–Alto", "Não significativo"),
fill = colors,
border = "grey60",
bty = "n",
cex = 1)
Estimação econométrica
Como mencionado, precisamos verificar se existe autocorrelação espacial nos resíduos do OLS
ols<-lm(percentual~pib_pc+gini_2010, data=df)
residuos<-resid(ols)
moran.test(x = residuos, listw = rook_listw)##
## Moran I test under randomisation
##
## data: residuos
## weights: rook_listw
## n reduced by no-neighbour observations
##
## Moran I statistic standard deviate = 12.625, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.3151503515 -0.0015552100 0.0006293325De fato, a autocorrelação espacial permaneceu presente. Portanto, seguimos para a modelagem espacial.
Spatial Autoregressive Model (SAR)
sar <- lagsarlm(percentual ~ pib_pc+gini_2010, data=df, listw = rook_listw, zero.policy = TRUE)
summary(sar)##
## Call:lagsarlm(formula = percentual ~ pib_pc + gini_2010, data = df,
## listw = rook_listw, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.96905 -5.47259 -0.93119 3.66115 40.72123
##
## Type: lag
## Regions with no neighbours included:
## 233
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1640e+00 2.6492e+00 0.4394 0.6604
## pib_pc 1.0628e-05 7.5490e-06 1.4079 0.1592
## gini_2010 2.2872e+01 5.5514e+00 4.1201 3.787e-05
##
## Rho: 0.50048, LR test value: 111.58, p-value: < 2.22e-16
## Asymptotic standard error: 0.044246
## z-value: 11.311, p-value: < 2.22e-16
## Wald statistic: 127.94, p-value: < 2.22e-16
##
## Log likelihood: -2270.775 for lag model
## ML residual variance (sigma squared): 63.307, (sigma: 7.9566)
## Number of observations: 645
## Number of parameters estimated: 5
## AIC: 4551.5, (AIC for lm: 4661.1)
## LM test for residual autocorrelation
## test value: 0.35259, p-value: 0.55265Efeitos marginais - SAR
Esses resultados são os coeficientes da regressão, que são, na maioria dos modelos espaciais, diferentes dos efeitos marginais. No entanto, \(\rho\) exerce um papel importante, o valor de 0,50 indica presença de dependência espacial na variável explicada.
Para encontrar os efeitos marginais, precisamos de:
## Impact measures (lag, exact):
## Direct Indirect Total
## pib_pc 1.128665e-05 9.990215e-06 2.127686e-05
## gini_2010 2.428879e+01 2.149888e+01 4.578766e+01
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 1.187e-05 8.481e-06 2.682e-07 2.620e-07
## gini_2010 2.420e+01 6.105e+00 1.931e-01 2.062e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -4.339e-06 6.286e-06 1.166e-05 1.76e-05 2.772e-05
## gini_2010 1.236e+01 2.002e+01 2.421e+01 2.83e+01 3.596e+01
##
## ========================================================
## Indirect:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 1.063e-05 7.869e-06 2.488e-07 2.488e-07
## gini_2010 2.173e+01 6.621e+00 2.094e-01 2.227e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -3.749e-06 5.511e-06 1.042e-05 1.60e-05 2.668e-05
## gini_2010 1.022e+01 1.722e+01 2.120e+01 2.54e+01 3.593e+01
##
## ========================================================
## Total:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 2.250e-05 1.619e-05 5.121e-07 4.806e-07
## gini_2010 4.593e+01 1.217e+01 3.847e-01 4.112e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -7.835e-06 1.212e-05 2.201e-05 3.358e-05 5.415e-05
## gini_2010 2.327e+01 3.747e+01 4.573e+01 5.355e+01 7.134e+01
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## pib_pc 8.480672e-06 7.868694e-06 1.619352e-05
## gini_2010 6.105313e+00 6.621027e+00 1.216591e+01
##
## Simulated z-values:
## Direct Indirect Total
## pib_pc 1.399608 1.350860 1.389391
## gini_2010 3.964315 3.281553 3.775354
##
## Simulated p-values:
## Direct Indirect Total
## pib_pc 0.16163 0.1767404 0.16471398
## gini_2010 7.3607e-05 0.0010324 0.00015978NOTA: p-valores estão sendo calculados por mil simulações de Monte Carlo, então podem variar um pouco a cada vez que rodamos o script.
Observe que os p-valores do pib per capita foram maiores que 10%, portanto não podemos rejeitar efeito nulo. No entanto, note que os p-valores do Gini apontam significância para o efeito direto, indireto (spillover) e total.
Dessa forma, um aumento de 0,01 no Gini do município \(i\), está associado a um aumento de 0,24 pontos percentuais no percentual de votos do partido no mesmo município. Além disso, temos ainda o efeito indireto, onde o aumento de 0,01 no Gini está associado a um aumento de 0,21 pontos percentuais no município vizinho. Por fim, o efeito agregado no sistema do aumento de 0,01 no Gini é de 0,46.
Spatial Error Model (SEM)
sem <- errorsarlm(percentual ~ pib_pc+gini_2010, data=df, listw = rook_listw, zero.policy = TRUE)
summary(sem)##
## Call:errorsarlm(formula = percentual ~ pib_pc + gini_2010, data = df,
## listw = rook_listw, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.22905 -5.48000 -0.82308 3.66159 40.61123
##
## Type: error
## Regions with no neighbours included:
## 233
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0726e+01 2.7416e+00 3.9124 9.138e-05
## pib_pc 8.8236e-06 7.9572e-06 1.1089 0.2675
## gini_2010 2.8936e+01 5.8952e+00 4.9085 9.178e-07
##
## Lambda: 0.52087, LR test value: 119.59, p-value: < 2.22e-16
## Asymptotic standard error: 0.044679
## z-value: 11.658, p-value: < 2.22e-16
## Wald statistic: 135.91, p-value: < 2.22e-16
##
## Log likelihood: -2266.769 for error model
## ML residual variance (sigma squared): 62.201, (sigma: 7.8868)
## Number of observations: 645
## Number of parameters estimated: 5
## AIC: 4543.5, (AIC for lm: 4661.1)No modelo SEM, como a dependência espacial está no erro, os coeficientes são os efeitos marginais, visto que não temos spillovers. O parâmetro \(\lambda\) de 0,52 mostra forte autocorrelação espacial positiva (e significativa) nos componentes não observados da variável dependente.
Spatial Durbin Model (SDM)
sdm <- lagsarlm(percentual ~ pib_pc+gini_2010, data=df, listw = rook_listw, type = "mixed", zero.policy = TRUE)
summary(sdm)##
## Call:lagsarlm(formula = percentual ~ pib_pc + gini_2010, data = df,
## listw = rook_listw, type = "mixed", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.21189 -5.46554 -0.76402 3.72130 40.87649
##
## Type: mixed
## Regions with no neighbours included:
## 233
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 9.0584e+00 3.7952e+00 2.3868 0.016996
## pib_pc 1.0061e-05 7.9982e-06 1.2579 0.208437
## gini_2010 2.9017e+01 5.8219e+00 4.9840 6.229e-07
## lag.pib_pc -4.7564e-06 1.5805e-05 -0.3009 0.763464
## lag.gini_2010 -2.3449e+01 8.1686e+00 -2.8707 0.004096
##
## Rho: 0.51687, LR test value: 118.14, p-value: < 2.22e-16
## Asymptotic standard error: 0.04474
## z-value: 11.553, p-value: < 2.22e-16
## Wald statistic: 133.47, p-value: < 2.22e-16
##
## Log likelihood: -2266.124 for mixed model
## ML residual variance (sigma squared): 62.142, (sigma: 7.883)
## Number of observations: 645
## Number of parameters estimated: 7
## AIC: 4546.2, (AIC for lm: 4662.4)
## LM test for residual autocorrelation
## test value: 10.626, p-value: 0.001115Efeitos marginais - SDM
## Impact measures (mixed, exact):
## Direct Indirect Total
## pib_pc 1.011828e-05 8.607849e-07 1.097906e-05
## gini_2010 2.791978e+01 -1.639643e+01 1.152336e+01
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 1.022e-05 8.143e-06 2.575e-07 2.770e-07
## gini_2010 2.745e+01 5.787e+00 1.830e-01 1.721e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -5.439e-06 4.572e-06 1.049e-05 1.536e-05 2.653e-05
## gini_2010 1.628e+01 2.361e+01 2.721e+01 3.137e+01 3.897e+01
##
## ========================================================
## Indirect:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 1.340e-06 2.875e-05 9.092e-07 6.556e-07
## gini_2010 -1.692e+01 1.443e+01 4.562e-01 4.562e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -5.519e-05 -1.819e-05 2.150e-06 0.0000218 5.537e-05
## gini_2010 -4.582e+01 -2.674e+01 -1.635e+01 -7.1285421 1.062e+01
##
## ========================================================
## Total:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 1.156e-05 3.087e-05 9.761e-07 9.761e-07
## gini_2010 1.053e+01 1.645e+01 5.203e-01 5.022e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -4.927e-05 -8.501e-06 1.118e-05 3.296e-05 6.983e-05
## gini_2010 -2.259e+01 -4.583e-01 1.115e+01 2.102e+01 4.289e+01
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## pib_pc 8.143417e-06 2.875246e-05 3.086554e-05
## gini_2010 5.786714e+00 1.442569e+01 1.645244e+01
##
## Simulated z-values:
## Direct Indirect Total
## pib_pc 1.255538 0.04661639 0.3746802
## gini_2010 4.742855 -1.17285751 0.6397997
##
## Simulated p-values:
## Direct Indirect Total
## pib_pc 0.20928 0.96282 0.7079
## gini_2010 2.1073e-06 0.24085 0.5223Spatial Durbin Error Model (SDEM)
sdem <- errorsarlm(percentual ~ pib_pc+gini_2010, data=df, listw = rook_listw, etype = "emixed", zero.policy = TRUE)
summary(sdem)##
## Call:errorsarlm(formula = percentual ~ pib_pc + gini_2010, data = df,
## listw = rook_listw, etype = "emixed", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.09845 -5.45291 -0.82333 3.66297 40.63172
##
## Type: error
## Regions with no neighbours included:
## 233
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4567e+01 5.9567e+00 2.4455 0.01447
## pib_pc 7.3973e-06 8.1234e-06 0.9106 0.36250
## gini_2010 2.8281e+01 6.0407e+00 4.6818 2.844e-06
## lag.pib_pc -1.0246e-05 2.0889e-05 -0.4905 0.62379
## lag.gini_2010 -6.5250e+00 1.0224e+01 -0.6382 0.52334
##
## Lambda: 0.5199, LR test value: 117.62, p-value: < 2.22e-16
## Asymptotic standard error: 0.044728
## z-value: 11.624, p-value: < 2.22e-16
## Wald statistic: 135.11, p-value: < 2.22e-16
##
## Log likelihood: -2266.383 for error model
## ML residual variance (sigma squared): 62.143, (sigma: 7.8831)
## Number of observations: 645
## Number of parameters estimated: 7
## AIC: 4546.8, (AIC for lm: 4662.4)Efeitos marginais - SDEM
## Impact measures (SDEM, glht, n):
## Direct Indirect Total
## pib_pc 7.397278e-06 -1.024574e-05 -2.848460e-06
## gini_2010 2.828103e+01 -6.525013e+00 2.175602e+01
## ========================================================
## Standard errors:
## Direct Indirect Total
## pib_pc 8.123350e-06 2.088894e-05 2.349117e-05
## gini_2010 6.040673e+00 1.022398e+01 1.294413e+01
## ========================================================
## Z-values:
## Direct Indirect Total
## pib_pc 0.9106192 -0.4904864 -0.1212566
## gini_2010 4.6817684 -0.6382065 1.6807637
##
## p-values:
## Direct Indirect Total
## pib_pc 0.3625 0.62379 0.903488
## gini_2010 2.8441e-06 0.52334 0.092809Kelejian-Prucha (KP)
kp <- sacsarlm(percentual ~ pib_pc+gini_2010, data=df, listw = rook_listw, zero.policy = TRUE)
summary(kp)##
## Call:sacsarlm(formula = percentual ~ pib_pc + gini_2010, data = df,
## listw = rook_listw, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.79682 -4.94136 -0.73082 3.39774 39.36817
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.2588e+01 3.9917e+00 5.6588 1.525e-08
## pib_pc 3.9307e-06 7.5203e-06 0.5227 0.6012
## gini_2010 2.6154e+01 5.5385e+00 4.7223 2.332e-06
##
## Rho: -0.44091
## Asymptotic standard error: 0.10495
## z-value: -4.2011, p-value: 2.656e-05
## Lambda: 0.75524
## Asymptotic standard error: 0.051194
## z-value: 14.752, p-value: < 2.22e-16
##
## LR test value: 129.51, p-value: < 2.22e-16
##
## Log likelihood: -2261.812 for sac model
## ML residual variance (sigma squared): 54.171, (sigma: 7.3601)
## Number of observations: 645
## Number of parameters estimated: 6
## AIC: 4535.6, (AIC for lm: 4661.1)Efeitos marginais - KP
## Impact measures (sac, exact):
## Direct Indirect Total
## pib_pc 4.067446e-06 -1.339497e-06 2.727949e-06
## gini_2010 2.706375e+01 -8.912671e+00 1.815108e+01
## ========================================================
## Simulation results ( variance matrix):
## Direct:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 3.957e-06 7.794e-06 2.465e-07 2.465e-07
## gini_2010 2.699e+01 5.530e+00 1.749e-01 1.749e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -1.143e-05 -1.351e-06 3.919e-06 8.966e-06 1.866e-05
## gini_2010 1.625e+01 2.310e+01 2.700e+01 3.080e+01 3.762e+01
##
## ========================================================
## Indirect:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc -1.251e-06 2.594e-06 8.202e-08 8.205e-08
## gini_2010 -8.829e+00 2.381e+00 7.530e-02 7.530e-02
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -6.521e-06 -2.889e-06 -1.190e-06 4.157e-07 3.684e-06
## gini_2010 -1.398e+01 -1.027e+01 -8.758e+00 -7.065e+00 -4.739e+00
##
## ========================================================
## Total:
##
## Iterations = 1:1000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## pib_pc 2.706e-06 5.286e-06 1.672e-07 1.672e-07
## gini_2010 1.816e+01 4.183e+00 1.323e-01 1.323e-01
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## pib_pc -7.489e-06 -9.209e-07 2.598e-06 6.093e-06 1.279e-05
## gini_2010 1.045e+01 1.523e+01 1.797e+01 2.093e+01 2.636e+01
##
## ========================================================
## Simulated standard errors
## Direct Indirect Total
## pib_pc 0.0000077943 2.593559e-06 5.285754e-06
## gini_2010 5.5298552934 2.381342e+00 4.182534e+00
##
## Simulated z-values:
## Direct Indirect Total
## pib_pc 0.5077169 -0.4824281 0.5119595
## gini_2010 4.8807900 -3.7075014 4.3421606
##
## Simulated p-values:
## Direct Indirect Total
## pib_pc 0.61165 0.62950185 0.60868
## gini_2010 1.0566e-06 0.00020931 1.4109e-05Escolhendo o melhor modelo
Podemos escolher o melhor modelo através do teste da razão de verossimilhança. No entanto,o teste de razão de verossimilhança não é aplicável quando os dois modelos não são aninhados, por exemplo, SAR e SEM, dessa forma, para tais casos a escolha do modelo se dá pelos valores de LL e AIC.
Hipótese nula do teste: o modelo mais simples é suficiente
Testando LR
## Likelihood ratio test
##
## Model 1: percentual ~ pib_pc + gini_2010
## Model 2: percentual ~ pib_pc + gini_2010
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -2266.8
## 2 7 -2266.1 2 1.2886 0.525## Likelihood ratio test
##
## Model 1: percentual ~ pib_pc + gini_2010
## Model 2: percentual ~ pib_pc + gini_2010
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -2266.8
## 2 6 -2261.8 1 9.9128 0.001641 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Likelihood ratio test
##
## Model 1: percentual ~ pib_pc + gini_2010
## Model 2: percentual ~ pib_pc + gini_2010
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -2266.8
## 2 7 -2266.4 2 0.7717 0.6799## Likelihood ratio test
##
## Model 1: percentual ~ pib_pc + gini_2010
## Model 2: percentual ~ pib_pc + gini_2010
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -2270.8
## 2 7 -2266.1 2 9.3015 0.009554 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## Likelihood ratio test
##
## Model 1: percentual ~ pib_pc + gini_2010
## Model 2: percentual ~ pib_pc + gini_2010
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 5 -2270.8
## 2 6 -2261.8 1 17.926 2.297e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Portanto, o melhor modelo é o KP.
Pseudo-R\(^2\)
Usamos o Pseudo-R\(^2\) pois não estamos no mundo de uma regressão linear estimada via OLS. \(\hat{y}\) não vem de uma projeção linear ortogonal, os resíduos não são i.i.d. e a variância de \(\hat{y}\) não se decompõe em variância explicada + erro. Assim, o valor do pseudo-R\(^2\) deve ser interpretado de forma diferente, não podemos mais usa-lo como a “fração da variância explicada”. O Pseudo-R\(^2\) serve apenas como indicador de quão bem \(\hat{y}\) aproxima de y.
y_obs <- df$percentual
y_pred <- kp$fitted.values
SSE <- sum((y_obs - y_pred)^2)
SST <- sum((y_obs - mean(y_obs))^2)
# Pseudo-R²
r2 <- 1 - SSE / SST
r2## [1] 0.3360343