Estatística Espacial

Universidade Estadual da Paraíba - UEPB

Professor: Ricardo Alves de Olinda

Aluno: Giullber Valentim da Silva

Apresentação dos mapas de quantil e Índice de Moran

Importação dos dados necessários

Shape municípios

shape_pb = readShapePoly('PB//PB_Municipios_2021.shp', IDvar = "NM_MUN")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

head(shape_pb,6)

##                CD_MUN        NM_MUN SIGLA AREA_KM2
## Água Branca   2500106   Água Branca    PB  241.662
## Aguiar        2500205        Aguiar    PB  351.607
## Alagoa Grande 2500304 Alagoa Grande    PB  322.071
## Alagoa Nova   2500403   Alagoa Nova    PB  128.230
## Alagoinha     2500502     Alagoinha    PB  111.361
## Alcantil      2500536      Alcantil    PB  309.896

Shape mesorregiões

shape_pb_meso = readShapePoly("PB meso//PB_Mesorregioes_2021.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

head(shape_pb_meso)

##   CD_MESO           NM_MESO SIGLA
## 0    2501  Sertão Paraibano    PB
## 1    2502         Borborema    PB
## 2    2503 Agreste Paraibano    PB
## 3    2504    Mata Paraibana    PB

Banco de dados com os indicadores

Foram utilizados os dados referentes ao indicadores IDHM, Gini e Theil obtidos no banco de dados da plataforma Atlas Brasil.

dados_ind = read_xlsx('banco de dados//mun_idhm.xlsx')
kable(head(dados_ind,6),
      caption='6 primeiras linhas do banco de dados contendo os indicadores por município.', align = 'c')

6 primeiras linhas do banco de dados contendo os indicadores por município.
Município	IDHM	GINI	THEIL
Água Branca	0.572	0.54	0.56
Aguiar	0.597	0.65	0.86
Alagoa Grande	0.582	0.53	0.53
Alagoa Nova	0.576	0.54	0.55
Alagoinha	0.595	0.52	0.50
Alcantil	0.578	0.48	0.46

Manipulações

No decorrer da análise foram necessários a tranformações dentro dos shape:

Nomes dos municípios do tipo fator para o tipo caracter;
Transformar apenas a primeira letra em maiúscula dos nomes dos municípios;
Modificação do nome do município de Quixabá

shape_pb$NM_MUN = as.character(shape_pb$NM_MUN)
shape_pb$NM_MUN = str_to_title(shape_pb$NM_MUN)
shape_pb$NM_MUN[156] = "Quixabá"

Realização do match

ind <- match(shape_pb@data$NM_MUN, dados_ind$Município)
ind

##   [1]   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18
##  [19]  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36
##  [37]  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54
##  [55]  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71  72
##  [73]  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89  90
##  [91]  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107 108
## [109] 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
## [127] 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
## [145] 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
## [163] 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## [181] 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
## [199] 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
## [217] 217 218 219 220 221 222 223

Realização da concatenação

## transformando em dataFrame
dados_ind = as.data.frame(dados_ind)
##
dados_ind <- dados_ind[ind,]
row.names(dados_ind) <- shape_pb$NM_MUN
row.names(dados_ind) = row.names(shape_pb)

## Realizando a concatenação
shape_pb <- spCbind(shape_pb, dados_ind)

kable(head(shape_pb,6),
      caption='6 primeiras linhas do banco de dados contendo os indicadores por município e as informações do shape.', align = 'c')

6 primeiras linhas do banco de dados contendo os indicadores por município e as informações do shape.
	CD_MUN	NM_MUN	SIGLA	AREA_KM2	Município	IDHM	GINI	THEIL
Água Branca	2500106	Água Branca	PB	241.662	Água Branca	0.572	0.54	0.56
Aguiar	2500205	Aguiar	PB	351.607	Aguiar	0.597	0.65	0.86
Alagoa Grande	2500304	Alagoa Grande	PB	322.071	Alagoa Grande	0.582	0.53	0.53
Alagoa Nova	2500403	Alagoa Nova	PB	128.230	Alagoa Nova	0.576	0.54	0.55
Alagoinha	2500502	Alagoinha	PB	111.361	Alagoinha	0.595	0.52	0.50
Alcantil	2500536	Alcantil	PB	309.896	Alcantil	0.578	0.48	0.46

Mapas

Mapa para IDHM

Os Intervalos dos quantis foram definidos com base na classificação do IDHM, como mostrado abaixo:

[0, 0.5) = Muito baixo
[0.5, 0.6) = Baixo
[0.6, 0,7) = Médio
[0.7, 0.8) = Alto
[0.8, 1] = Muito alto

par(mar=c(1,0,1,0))

# Definindo os intevarlos dos quantis e cores associados
INT <- classIntervals(shape_pb$IDHM, style="fixed", 
                      fixedBreaks=c(0,0.500,0.600,0.700,0.800,1))
CORES1 <- brewer.pal(9, "YlGn")
COL1 <-  findColours(INT, CORES1)

# plot
plot(shape_pb, col=COL1)
plot(shape_pb_meso, add=TRUE,lwd=2)

# Legenda, título e mais informações
TB1 <- attr(COL1, "table")
names(TB1)= c('Muito Baixo','Baixo','Médio','Alto','Muito Alto')
legtext <- paste(names(TB1))
legend(-35.32445, -7.684392, fill=attr(COL1, "palette"), legend=legtext, 
       bty="n",cex=0.8, title='Legenda', title.font = 2)
title("IDHM do estado da Paraíba de acordo com o censo de 2010")
scalebar(100, xy=c(-36.37951, -8.080096), 
         type="bar", below="km",
         cex=0.8, lonlat=T,divs=4)
compassRose(-37.73203, -8.065287, cex=0.9)

Mapa para Gini

par(mar=c(1,0,2,0))

# Definindo os intevarlos dos quantis e cores associados
INT2 <- classIntervals(shape_pb$GINI, n=4, style="quantile")
CORES2 <- brewer.pal(9, "YlOrBr")
COL2 <-  findColours(INT2, CORES2)

# plot
plot(shape_pb, col=COL2)
plot(shape_pb_meso, add=TRUE,lwd=2)

# Legenda, título e mais informações
TB2 <- attr(COL2, "table")
legtext <- paste(names(TB2))
legend(-35.32445, -7.684392, fill=attr(COL2, "palette"), legend=legtext, 
       bty="n",cex=0.8, title='Legenda', title.font = 2)
title("Coeficiente de Gini por município no estado da Paraíba 
      de acordo com o censo de 2010")
scalebar(100, xy=c(-36.37951, -8.080096), 
         type="bar", below="km",
         cex=0.8, lonlat=T,divs=4)
compassRose(-37.73203, -8.065287, cex=0.9)

Mapa para Theil

par(mar=c(1,0,2,0))

# Definindo os intevarlos dos quantis e cores associados
INT3 <- classIntervals(shape_pb$THEIL, n=4, style="quantile")
CORES3 <- brewer.pal(9, "Blues")
COL3 <-  findColours(INT3, CORES3)

# plot
plot(shape_pb, col=COL3)
plot(shape_pb_meso, add=TRUE,lwd=2)

# Legenda, título e mais informações
TB3 <- attr(COL3, "table")
legtext <- paste(names(TB3))
legend(-35.32445, -7.684392, fill=attr(COL3, "palette"), legend=legtext, 
       bty="n",cex=0.8, title='Legenda', title.font = 2)
title("Índice de Theil por município no estado da Paraíba 
      de acordo com o censo de 2010")
scalebar(100, xy=c(-36.37951, -8.080096), 
         type="bar", below="km",
         cex=0.8, lonlat=T,divs=4)
compassRose(-37.73203, -8.065287, cex=0.9)

Matriz de Vizinhança

shape_pb.nb1 <- poly2nb(shape_pb)

vizinhanca = nb2listw(shape_pb.nb1, style="W",
                      zero.policy=TRUE)
vizinhanca

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 223 
## Number of nonzero links: 1206 
## Percentage nonzero weights: 2.425144 
## Average number of links: 5.408072 
## 
## Weights style: W 
## Weights constants summary:
##     n    nn  S0       S1       S2
## W 223 49729 223 90.71447 932.6326

Índice de Moran

Como a aplicação do índice de Moran é realizada para ientificar uma autocorrelação espacial, podemos utilizar as seguintes hipóteses de interesse a serem testadas:

\(H_{0}\): Não há evidências de autocorrelação espacial
\(H_{1}\): Há evidências de autocorrelação espacial

onde iremos rejeitar a \(H_{0}\) quando o p-valor < 0.05, considerando 5% como o nível de significância para este teste. Portanto:

Cálculo do índice de Moran global

Para o IDHM

Mglobal1 = moran.test(shape_pb$IDHM, listw=nb2listw(shape_pb.nb1))
Mglobal1

## 
##  Moran I test under randomisation
## 
## data:  shape_pb$IDHM  
## weights: nb2listw(shape_pb.nb1)    
## 
## Moran I statistic standard deviate = 8.0641, p-value = 3.689e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.332983019      -0.004504505       0.001751474

# Teste do pvalor
Mglobal1$p.value<0.05

## [1] TRUE

Para o índice de Gini

Mglobal2 = moran.test(shape_pb$GINI, listw=nb2listw(shape_pb.nb1))
Mglobal2

## 
##  Moran I test under randomisation
## 
## data:  shape_pb$GINI  
## weights: nb2listw(shape_pb.nb1)    
## 
## Moran I statistic standard deviate = 2.2946, p-value = 0.01088
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.091900288      -0.004504505       0.001765148

# Teste do pvalor
Mglobal2$p.value<0.05

## [1] TRUE

Para o índice de Theil

Mglobal3 = moran.test(shape_pb$THEIL, listw=nb2listw(shape_pb.nb1))
Mglobal3

## 
##  Moran I test under randomisation
## 
## data:  shape_pb$THEIL  
## weights: nb2listw(shape_pb.nb1)    
## 
## Moran I statistic standard deviate = 3.7185, p-value = 0.0001002
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.151055632      -0.004504505       0.001750074

# Teste do pvalor
Mglobal3$p.value<0.05

## [1] TRUE

O que se observa é que no cálculo do índice de Moran global, todos os testes realizados para cada indicador apresentaram p-valores menores que o nível de 5% de significância, ou seja, existem evidências que corroboram para a existência da autocorrelação espacial.

Cálculo do índice de Moran local

Para o IDHM

shape_pb.mloc1 <- localmoran(shape_pb$IDHM, listw=vizinhanca,
                            zero.policy=T, 
                            alternative = "two.sided")
head(shape_pb.mloc1,6)

##                         Ii          E.Ii      Var.Ii        Z.Ii Pr(z != E(Ii))
## Água Branca   2.343750e-01 -0.0007751543 0.057054096 0.984468723      0.3248851
## Aguiar        7.847076e-05 -0.0002762383 0.006594981 0.004367828      0.9965150
## Alagoa Grande 6.436525e-02 -0.0001011390 0.002415040 1.311810672      0.1895840
## Alagoa Nova   6.994724e-03 -0.0004296363 0.018806878 0.054137813      0.9568254
## Alagoinha     1.553492e-02 -0.0001706539 0.006198081 0.199491916      0.8418780
## Alcantil      1.569747e-01 -0.0002948333 0.016209049 1.235281123      0.2167259

Para o índice de Gini

shape_pb.mloc2 <- localmoran(shape_pb$GINI, listw=vizinhanca,
                            zero.policy=T, 
                            alternative = "two.sided")
head(shape_pb.mloc2,6)

##                       Ii          E.Ii     Var.Ii       Z.Ii Pr(z != E(Ii))
## Água Branca    0.5840943 -0.0039487993 0.28972247  1.0924916     0.27461709
## Aguiar        -0.1866800 -0.0556214351 1.25440563 -0.1170163     0.90684713
## Alagoa Grande -0.1652293 -0.0022195380 0.05288680 -0.7088266     0.47843206
## Alagoa Nova    0.1745396 -0.0039487993 0.17224596  0.4300663     0.66714739
## Alagoinha      0.3282457 -0.0009849851 0.03574512  1.7413737     0.08161810
## Alcantil       0.3857440 -0.0009938588 0.05460117  1.6550673     0.09791084

Para o índice de Theil

shape_pb.mloc3 <- localmoran(shape_pb$THEIL, listw=vizinhanca,
                            zero.policy=T, 
                            alternative = "two.sided")
head(shape_pb.mloc3,6)

##                        Ii          E.Ii      Var.Ii        Z.Ii Pr(z != E(Ii))
## Água Branca    0.72987140 -0.0034379465 0.252370690  1.45971397     0.14436870
## Aguiar         0.25503127 -0.0836600142 1.830729819  0.25031789     0.80234152
## Alagoa Grande -0.19796586 -0.0012654643 0.030182113 -1.13221895     0.25754240
## Alagoa Nova   -0.03235907 -0.0025956095 0.113373866 -0.08839484     0.92956286
## Alagoinha      0.07390497 -0.0001565683 0.005686576  0.98212618     0.32603769
## Alcantil       0.29163888 -0.0003325078 0.018279592  2.15951909     0.03080992

Amplitude do Índice de Moran local

Para o IDHM

list_w <- vizinhanca
Sd_1 <- (shape_pb$IDHM) - mean(shape_pb$IDHM)
mI_1 <- shape_pb.mloc1[, 5]
C_mI <- mI_1 - mean(mI_1) 

# Definindo os quadrantes
quadrant <- vector(mode = "numeric", length = nrow(shape_pb.mloc1))
quadrant[Sd_1 > 0 & mI_1 > 0] <- 1
quadrant[Sd_1 < 0 & mI_1 > 0] <- 2
quadrant[Sd_1 > 0 & mI_1 < 0] <- 3
quadrant[Sd_1 < 0 & mI_1 < 0] <- 4

signif <- 0.05
quadrant[shape_pb.mloc1[, 5] > signif] <- 5

#
colors <- c("red", "blue", "lightpink", "skyblue2", "white")

# plot
plot(shape_pb, col=colors[quadrant])
plot(shape_pb_meso, add=TRUE,lwd=2)

# Legenda, título e mais informações
legend(-35.32445, -7.684392, 
       legend = c("alto-alto", "baixo-baixo", "alto-baixo",
                  "baixo-alto","Não Sig."), 
       fill = colors, bty = "n", cex = 0.7, y.intersp = 1, x.intersp = 1,
       title = 'Legenda',title.font = 2)
title("Índice de Moran Local - IDHM 2010")
scalebar(100, xy=c(-36.37951, -8.080096), 
         type="bar", below="km",
         cex=0.8, lonlat=T,divs=4)
compassRose(-37.73203, -8.065287, cex=0.9)

Para o índice de Gini

Sd_2 <- (shape_pb$GINI) - mean(shape_pb$GINI)
mI_2 <- shape_pb.mloc2[, 5]
C_m2 <- mI_2 - mean(mI_2) 

# Definindo os quadrantes
quadrant2 <- vector(mode = "numeric", length = nrow(shape_pb.mloc2))
quadrant2[Sd_2 > 0 & mI_2 > 0] <- 1
quadrant2[Sd_2 < 0 & mI_2 > 0] <- 2
quadrant2[Sd_2 > 0 & mI_2 < 0] <- 3
quadrant2[Sd_2 < 0 & mI_2 < 0] <- 4

signif <- 0.05
quadrant2[shape_pb.mloc2[, 5] > signif] <- 5

#
colors2 <- c("red", "blue", "lightpink", "skyblue2", "white")

# plot
plot(shape_pb, col=colors2[quadrant2])
plot(shape_pb_meso, add=TRUE,lwd=2)

# Legenda, título e mais informações
legend(-35.32445, -7.684392, 
       legend = c("alto-alto", "baixo-baixo", "alto-baixo",
                  "baixo-alto","Não Sig."), 
       fill = colors, bty = "n", cex = 0.7, y.intersp = 1, x.intersp = 1,
       title = 'Legenda',title.font = 2)
title("Índice de Moran Local - Gini 2010")
scalebar(100, xy=c(-36.37951, -8.080096), 
         type="bar", below="km",
         cex=0.8, lonlat=T,divs=4)
compassRose(-37.73203, -8.065287, cex=0.9)

Para o índice de Theil

Sd_3 <- (shape_pb$THEIL) - mean(shape_pb$THEIL)
mI_3 <- shape_pb.mloc3[, 5]
C_m3 <- mI_3 - mean(mI_3) 

# Definindo os quadrantes
quadrant3 <- vector(mode = "numeric", length = nrow(shape_pb.mloc3))
quadrant3[Sd_3 > 0 & mI_3 > 0] <- 1
quadrant3[Sd_3 < 0 & mI_3 > 0] <- 2
quadrant3[Sd_3 > 0 & mI_3 < 0] <- 3
quadrant3[Sd_3 < 0 & mI_3 < 0] <- 4

signif <- 0.05
quadrant3[shape_pb.mloc3[, 5] > signif] <- 5

#
colors3 <- c("red", "blue", "lightpink", "skyblue2", "white")

# plot
plot(shape_pb, col=colors3[quadrant3])
plot(shape_pb_meso, add=TRUE,lwd=2)

# Legenda, título e mais informações
legend(-35.32445, -7.684392, 
       legend = c("alto-alto", "baixo-baixo", "alto-baixo",
                  "baixo-alto","Não Sig."), 
       fill = colors, bty = "n", cex = 0.7, y.intersp = 1, x.intersp = 1,
       title = 'Legenda',title.font = 2)
title("Índice de Moran Local - Theil 2010")
scalebar(100, xy=c(-36.37951, -8.080096), 
         type="bar", below="km",
         cex=0.8, lonlat=T,divs=4)
compassRose(-37.73203, -8.065287, cex=0.9)