15 Anexo II - Exemplos de Analise Espacial

Neste anexo utilizaremos dois exemplos para estudar análise espacial. Ao final do anexo, apresentamos formas de obter os centroides para lidar com os casos em que os centroides caem fora do polígono.

ATENÇÃO

Devido a maneira que o Windows acessa as url usando a internet é necessario mudar a opção default dele para que possa usar apropiadamente os recursos https, ftps, etc… A linha abaixo deve ser utilizanda antes de usar as funções que acessam esse tipo de recurso.

options(download.file.method='libcurl',url.method='libcurl')

Carregando as bibliotecas necessárias:

library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)
library(spdep)
library(RColorBrewer)
library(colorspace)
library(spgwr)
library(RColorBrewer)

15.1 Exemplo com os dados de dengue em Dourados/MS

Neste exeplo serão utilizados os dados da monografia de Isis Rodrigues Reitman, apresentada ao Curso de Geografia da Faculdade de Ciências Humanas da Universidade Federal da Grande Douradosos/MS, em março de 2013. O título da monografia é “DISTRIBUIÇÃO ESPACIAL DOS CASOS DE DENGUE NO PERÍMETRO URBANO DE DOURADOS-MS E SUA RELAÇÃO COM OS FATORES SOCIOAMBIENTAIS E POLÍTICOS”.

Lendo a tabela da população por setor censitário e baixando os shapefiles do contorno e por setor censitário de Dourados/MS:

local <- "https://gitlab.procc.fiocruz.br/oswaldo/eco2019/raw/master/dados/"
pop2010 <- read_csv(paste0(local, "pop2010.csv"))

tmpdir <- tempdir()
download.file(paste0(local, "setores_dourados.zip"),
    destfile = paste0(tmpdir, "/dourados.zip"))

unzip(zipfile = paste0(tmpdir, "/dourados.zip"), exdir = tmpdir)
dir(tmpdir)

##  [1] "contorno.dbf"               "contorno.sbn"               "contorno.sbx"               "contorno.shp"               "contorno.shx"               "dourados.zip"               "file197ae4dc44c71"          "leaflet-providers_1.9.0.js" "olinda.dbf"                 "olinda.shp"                
## [11] "olinda.shx"                 "olinda.zip"                 "Setor_UTM_SIRGAS.dbf"       "Setor_UTM_SIRGAS.prj"       "Setor_UTM_SIRGAS.sbn"       "Setor_UTM_SIRGAS.sbx"       "Setor_UTM_SIRGAS.shp"       "Setor_UTM_SIRGAS.shx"       "vroom-197ae16fe30e1"        "vroom-197ae299fee3c"       
## [21] "vroom-197ae51a0f3c8"

setor.sf <- read_sf(paste0(tmpdir, "/Setor_UTM_SIRGAS.shp"),
    crs = 31981)
contorno.sf <- read_sf(paste0(tmpdir, "/contorno.shp"),
    crs = 31981)

Juntando com os dados de setores censitários e população:

setor.sf <- setor.sf %>%
    mutate(idsetor = as.numeric(CD_GEOCODI)) %>%
    inner_join(pop2010, by = "idsetor")

Lendo e plotando os casos de dengue georreferenciados em Dourados/MS:

casos <- read_csv(paste0(local, "dengue_dourados.csv"))
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") +
    geom_sf(data = casos.sf, color = "red", size = 1) +
    theme_void()

Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS:

lixo <- read_csv2(paste0(local, "/lixo_dourados.csv"))
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") +
    geom_sf(data = lixo.sf, color = "blue", size = 1) +
    theme_void()

Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS. Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.

lixo2.sf <- st_intersection(contorno.sf, lixo.sf)
# ou lixo2.sf <- st_intersection(setor.sf,
# lixo.sf)

ggplot(setor.sf) + geom_sf(fill = "white", color = "black") +
    geom_sf(data = lixo2.sf, color = "blue", size = 1) +
    theme_void()

Verificando as paletas de cores.

display.brewer.all()

Utilizando as informações dos casos (pontos), do lixo (ponto) e da população de cada setor censitário (mapa temático):

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf,
    size = 0.7, alpha = 0.9, aes(colour = "Caso"),
    show.legend = "point") + geom_sf(data = lixo2.sf,
    size = 1, alpha = 0.9, aes(colour = "Lixo"), show.legend = "point") +
    scale_fill_distiller(palette = "PuBu", direction = 1) +
    scale_colour_manual(values = c(Caso = "red", Lixo = "green")) +
    theme_minimal()

Iremos agora construir buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue?

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500)

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf,
    size = 0.7, alpha = 0.9, aes(colour = "Caso"),
    show.legend = "point") + geom_sf(data = lixo2.sf,
    size = 1, alpha = 0.9, aes(colour = "Lixo"), show.legend = "point") +
    scale_fill_distiller(palette = "PuBu", direction = 1) +
    geom_sf(data = lixo_buffer, color = "grey20", fill = "transparent",
        size = 0.4) + scale_colour_manual(values = c(Caso = "red",
    Lixo = "green")) + theme_minimal()

Represntando os casos e o lixo de forma interativa.

tm_shape(setor.sf) + tm_borders("black") + tm_shape(casos.sf) +
    tm_dots("red") + tm_shape(lixo.sf) + tm_dots("green") +
    tm_shape(lixo_buffer) + tm_borders("blue") + tmap_mode("view")

## conta casos por setor 
casos.sf$contador <- 1 
casos <- setor.sf %>%  
  st_join(casos.sf) %>% 
  filter(CLASSI_FIN == 1) %>%  ## seleciona somente os casos confirmados
  group_by(ID) %>% 
  summarise(casos = sum(contador))

st_geometry(casos) <- NULL  ## remove atributos de geometria

## numero de depositos de Lixo por setor 
lixo.sf$contador <- 1 

lixo <- setor.sf %>%  
  st_join(lixo.sf) %>% 
  group_by(ID) %>% 
  summarise(lixo = sum(contador))

st_geometry(lixo) <- NULL ## remove atributos de geometria

# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf %>% 
  left_join(lixo,by = 'ID') %>% 
  left_join(casos,by = 'ID')

Plotando o mapa temático dos casos por setor censitário:

plot(setor.sf["casos"])

Plotando o mapa temático dos pontos de coleta de lixo por setor censitário:

plot(setor.sf["lixo"])

Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário:

setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0  # Transformando os missings em zero

summary(setor.sf$tx)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    1.74    4.06    4.31   56.40

pal <- wesanderson::wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + geom_sf(aes(fill = tx), color = "black") +
    scale_fill_gradientn(colours = pal) + ggtitle("Taxa de incidência de Dengue") +
    theme_void()

15.2 Kernel por atributos

Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.

Primeiramente é necessário dissolver os poligonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

dourados.w <- as.owin(as_Spatial(dourados.contorno))  #

Extraindo os centróidas dos polígonos em Dourados/MS.

centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")

plot(centroides.sp)

Colocando os pontos no formato sp:

centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y,
    dourados.w)
plot(centroides.ppp, pch = 19, cex = 0.5)

Fazendo o kernel por atributo da taxa:

kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx,
    scalekernel = TRUE)
plot(kernel.tx)

15.3 Exemplo com os dados de dengue em Pernambuco

Nesta prática serão utilizados os dados de dengue por município no estado de Pernambuco. Os dados foram organizados por Freitas et al. (2021).

Agora vamos ler os dados e o polígono de Pernambuco.

Verificando os objetos:

head(dzc)

## # A tibble: 6 × 8
##   ID_MUNICIP  year Population dengue  zika chikungunya microcephaly live_births
##        <dbl> <dbl>      <dbl>  <dbl> <dbl>       <dbl>        <dbl>       <dbl>
## 1     260005  2013     96769.     14     0           0            0        1532
## 2     260005  2014     96831.     23     0           0            0        1526
## 3     260005  2015     96907.     54     0           0            5        1628
## 4     260005  2016     96880.     23     0          29            2        1487
## 5     260005  2017     96828.      0     0          10            0        1519
## 6     260010  2013     36182.      0     0           0            0         638

summary(dzc)

##    ID_MUNICIP          year        Population          dengue           zika        chikungunya    microcephaly    live_births   
##  Min.   :260005   Min.   :2013   Min.   :   4475   Min.   :    0   Min.   : 0.00   Min.   :   0   Min.   : 0.00   Min.   :   50  
##  1st Qu.:260419   1st Qu.:2014   1st Qu.:  14295   1st Qu.:    0   1st Qu.: 0.00   1st Qu.:   0   1st Qu.: 0.00   1st Qu.:  203  
##  Median :260822   Median :2015   Median :  22735   Median :    5   Median : 0.00   Median :   0   Median : 0.00   Median :  334  
##  Mean   :260825   Mean   :2015   Mean   :  50606   Mean   :  140   Mean   : 0.18   Mean   :  36   Mean   : 0.89   Mean   :  757  
##  3rd Qu.:261232   3rd Qu.:2016   3rd Qu.:  38049   3rd Qu.:   38   3rd Qu.: 0.00   3rd Qu.:   1   3rd Qu.: 1.00   3rd Qu.:  594  
##  Max.   :261650   Max.   :2017   Max.   :1599273   Max.   :24793   Max.   :93.00   Max.   :9209   Max.   :87.00   Max.   :23664

PE

## Simple feature collection with 184 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -41.36 ymin: -9.483 xmax: -34.81 ymax: -7.273
## Geodetic CRS:  WGS 84
## # A tibble: 184 × 3
##    NM_MUNICIP            CD_GEOCMU                                                                                         geom
##  * <chr>                 <chr>                                                                                    <POLYGON [°]>
##  1 ABREU E LIMA          2600054   ((-35.05 -7.88, -35.02 -7.888, -35.02 -7.888, -35 -7.894, -35 -7.895, -34.98 -7.9, -34.94...
##  2 AFOGADOS DA INGAZEIRA 2600104   ((-37.56 -7.639, -37.56 -7.64, -37.56 -7.641, -37.56 -7.642, -37.56 -7.643, -37.56 -7.644...
##  3 AFRÂNIO               2600203   ((-40.86 -8.371, -40.85 -8.446, -40.85 -8.459, -40.85 -8.462, -40.85 -8.526, -40.85 -8.55...
##  4 AGRESTINA             2600302   ((-35.92 -8.366, -35.9 -8.375, -35.9 -8.375, -35.9 -8.381, -35.9 -8.384, -35.9 -8.384, -3...
##  5 ÁGUA PRETA            2600401   ((-35.47 -8.546, -35.47 -8.545, -35.47 -8.549, -35.46 -8.549, -35.46 -8.549, -35.46 -8.54...
##  6 ÁGUAS BELAS           2600500   ((-36.87 -8.887, -36.87 -8.888, -36.87 -8.893, -36.87 -8.895, -36.87 -8.905, -36.87 -8.90...
##  7 ALAGOINHA             2600609   ((-36.76 -8.599, -36.76 -8.6, -36.77 -8.603, -36.77 -8.602, -36.77 -8.602, -36.78 -8.603,...
##  8 ALIANÇA               2600708   ((-35.16 -7.503, -35.15 -7.506, -35.15 -7.507, -35.15 -7.509, -35.15 -7.51, -35.14 -7.51,...
##  9 ALTINHO               2600807   ((-35.97 -8.536, -36 -8.547, -36 -8.547, -36 -8.547, -36 -8.534, -36 -8.534, -36 -8.534, ...
## 10 AMARAJI               2600906   ((-35.41 -8.428, -35.41 -8.427, -35.41 -8.427, -35.41 -8.426, -35.41 -8.426, -35.41 -8.42...
## # … with 174 more rows

ggplot(PE) + geom_sf(size = 0.3, fill = "darkseagreen1")

Calculando a incidência de dengue acumulada 2013-2017 e verificando a distribuição:

dengue <- dzc %>%
    group_by(ID_MUNICIP) %>%
    summarise(pop = mean(Population), dengue = sum(dengue),
        inc = (dengue/pop) * 10000)

par(mfrow = c(1, 2))
hist(dengue$inc)
hist(log(dengue$inc))

par(mfrow = c(1, 1))
ggplot(dengue, aes(log(inc))) + geom_histogram(aes(y = ..density..),
    colour = "black", fill = "grey90") + geom_density(fill = "purple",
    alpha = 0.2) + theme_minimal()

Agora vamos juntar os dados de dengue com o polígono de Pernambuco. Observe que a variável ID_MUNICIP do objeto com os dados dos casos tem 6 dígitos, enquanto que a variável CD_GEOCMU do polígono tem 7 dígitos. Porém, os 6 primeiros dígitos de CD_GEOCMU são correspondentes à variável ID_MUNICIP. É comum encontrar essa situação quando trabalhamos com dados do IBGE e do SINAN.

Portanto, vamos criar a variável ID_MUNICIP no polígono, e juntar com o objeto de casos:

PE2 <- PE %>%
    rename(ID_MUNICIP = CD_GEOCMU) %>%
    mutate(ID_MUNICIP = as.numeric(str_sub(ID_MUNICIP,
        1, 6))) %>%
    left_join(dengue, by = "ID_MUNICIP")

Agora vamos ver o mapa das incidências:

g1 <- ggplot(PE2) + geom_sf(aes(fill = cut_number(inc,
    5)), size = 0.2) + scale_fill_discrete_sequential(palette = "Reds",
    name = "casos\npor 10000") + labs(title = "Incidência de dengue em Pernambuco, 2013-2017") +
    theme_void()

g2 <- ggplot(PE2, aes(inc)) + geom_density(fill = "red",
    alpha = 0.2) + geom_vline(xintercept = quantile(PE2$inc,
    probs = seq(0, 1, by = 0.2)), color = "red", linetype = "dashed") +
    theme(plot.margin = margin(3, 1, 1, 1, "cm"))

ggpubr::ggarrange(g1, g2, nrow = 1, widths = c(6.5,
    3.5))

15.4 Matriz de Vizinhança

Iremos precisar da coordenadas dos centróides:

centroides <- st_centroid(PE)
coordenadas <- st_coordinates(centroides)

Criando e verificando a vizinhança de Pernambuco:

# Convertendo para Spatial Polygons para usar
# poly2nb
PE.sp <- as_Spatial(PE)

# Criando a vizinhança
viz <- poly2nb(PE.sp)
viz

## Neighbour list object:
## Number of regions: 184 
## Number of nonzero links: 948 
## Percentage nonzero weights: 2.8 
## Average number of links: 5.152

Fazendo o mapa com a vizinhança:

# Transformando a vizinhança em linhas para
# plotar
viz.sf <- nb2lines(viz, coords = coordenadas, as_sf = TRUE)
viz.sf <- st_set_crs(viz.sf, st_crs(PE))

# Plotando o mapa de conectividade por
# contiguidade
ggplot() + geom_sf(data = PE, fill = "salmon", size = 0.2,
    color = "white") + geom_sf(data = centroides, size = 1) +
    geom_sf(data = viz.sf, size = 0.2) + ggtitle("Vizinhança por conectividade") +
    ylab("Latitude") + xlab("Longitude") + theme_minimal()

15.5 Autocorrelação Espacial

Verificando a correlação da incidência de dengue em Pernambuco.

pesos.viz <- nb2listw(viz)
moran.test(PE2$inc, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  PE2$inc  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 3.2, p-value = 0.0006
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##          0.140386         -0.005464          0.002022

Plotando o correlograma:

correl <- sp.correlogram(viz, PE2$inc, order = 8, method = "I")
correl

## Spatial correlogram for PE2$inc 
## method: Moran's I
##          estimate expectation  variance standard deviate Pr(I) two sided   
## 1 (184)  0.140386   -0.005464  0.002022             3.24          0.0012 **
## 2 (184)  0.008862   -0.005464  0.000977             0.46          0.6467   
## 3 (184)  0.050621   -0.005464  0.000724             2.08          0.0372 * 
## 4 (184)  0.017289   -0.005464  0.000643             0.90          0.3694   
## 5 (184) -0.048713   -0.005464  0.000640            -1.71          0.0873 . 
## 6 (184) -0.014771   -0.005464  0.000629            -0.37          0.7106   
## 7 (184)  0.031309   -0.005464  0.000617             1.48          0.1387   
## 8 (184) -0.020329   -0.005464  0.000669            -0.57          0.5655   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

plot(correl)

Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local:

PE2$pval <- localmoran(PE2$inc, pesos.viz)[, 5]

tm_shape(PE2) + tm_polygons(col = "pval", title = "p-valores",
    breaks = c(0, 0.01, 0.05, 0.1, 1), border.col = "white",
    palette = "-Oranges") + tm_scale_bar(width = 0.15) +
    tm_layout(frame = FALSE)

Moran Local (Lisa Map) da incidência:

resI <- localmoran.sad(lm(PE2$inc ~ 1), 1:length(viz),
    viz, style = "W")
summary(resI)[1:5, ]

##     Local Morans I Stand. dev. (N) Pr. (N) Saddlepoint Pr. (Sad) Expectation Variance Skewness Kurtosis Minimum Maximum       omega     sad.r     sad.u
## 1 1      -0.122210        -0.26537  0.6046    -0.44298    0.6711   -0.005464   0.1935 -0.07213    8.624  -40.97   39.97 -0.00303917 -0.259013 -0.271654
## 2 2      -0.170546        -0.33470  0.6311    -0.53962    0.7053   -0.005464   0.2433 -0.06434    8.624  -45.87   44.87 -0.00331024 -0.323850 -0.347289
## 3 3       0.072842         0.11165  0.4556     0.18386    0.4271   -0.005464   0.4919 -0.04525    8.622  -65.02   64.02  0.00085244  0.110805  0.111706
## 4 4      -0.007457        -0.00453  0.5018    -0.02005    0.5080   -0.005464   0.1935 -0.07213    8.624  -40.97   39.97 -0.00005565 -0.004505 -0.004506
## 5 5       0.299109         0.88311  0.1886     1.11828    0.1317   -0.005464   0.1189 -0.09200    8.627  -32.23   31.23  0.00941159  0.788180  1.022395

PE2$MoranLocal <- summary(resI)[, 1]

ggplot(PE2) + geom_sf(aes(fill = MoranLocal), color = "black",
    size = 0.2) + scale_fill_continuous_diverging(palette = "Blue-Red 3",
    mid = 0) + ggtitle("Moran local") + theme_void()

15.6 Modelos Linear, CAR e GWR

Agora vamos utilizar o IDH para explorar os modelos.

Primeiro, vamos visualizar o IDH por município e comparar com o mapa de incidência de dengue.

idh <- read_csv("https://gitlab.procc.fiocruz.br/oswaldo/eco2019/-/raw/master/dados/idh_pernambuco.csv")
head(idh)

## # A tibble: 6 × 2
##   ID_MUNICIP  IDHM
##        <dbl> <dbl>
## 1     260005 0.679
## 2     260010 0.657
## 3     260020 0.588
## 4     260030 0.592
## 5     260040 0.553
## 6     260050 0.526

PE2 <- PE2 %>%
    left_join(idh)

g3 <- ggplot(PE2) + geom_sf(aes(fill = IDHM), size = 0.2) +
    scale_fill_continuous_sequential(palette = "Viridis",
        name = "IDH") + labs(title = "IDH, 2010") +
    theme_void()

ggpubr::ggarrange(g1, g3, nrow = 2)

Ajustando o modelo de regressão linear simples.

dengue.lm <- lm(log(inc + 1) ~ IDHM, data = PE2)
summary(dengue.lm)

## 
## Call:
## lm(formula = log(inc + 1) ~ IDHM, data = PE2)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.740 -1.001  0.065  1.142  3.538 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    0.576      1.414    0.41    0.684  
## IDHM           5.474      2.369    2.31    0.022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.46 on 182 degrees of freedom
## Multiple R-squared:  0.0285, Adjusted R-squared:  0.0232 
## F-statistic: 5.34 on 1 and 182 DF,  p-value: 0.0219

Checando os residuos para verificar a presença de autocorrelação.

dengue.lm$lmresid <- residuals(dengue.lm)
moran.test(dengue.lm$lmresid, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  dengue.lm$lmresid  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = 3.8, p-value = 0.00007
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##          0.174454         -0.005464          0.002233

Ajustando o modelo CAR.

dengue.car <- spautolm(log(inc + 1) ~ IDHM, data = PE2,
    listw = pesos.viz, family = "CAR")
summary(dengue.car)

## 
## Call: spautolm(formula = log(inc + 1) ~ IDHM, data = PE2, listw = pesos.viz,     family = "CAR")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -3.957162 -0.944274  0.069071  1.082266  3.519477 
## 
## Coefficients: 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  -1.7043     1.5392 -1.1073 0.2681615
## IDHM          9.1806     2.5812  3.5567 0.0003756
## 
## Lambda: 0.4945 LR test value: 9.016 p-value: 0.0026767 
## Numerical Hessian standard error of lambda: 0.1806 
## 
## Log likelihood: -324.8 
## ML residual variance (sigma squared): 1.943, (sigma: 1.394)
## Number of observations: 184 
## Number of parameters estimated: 4 
## AIC: 657.7

Checando os residuos para verificar a presença de autocorrelação.

dengue.car$carresid <- residuals(dengue.car)
moran.test(dengue.car$carresid, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  dengue.car$carresid  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = -1.5, p-value = 0.9
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         -0.078140         -0.005464          0.002231

Ajustando o modelo GWR.

Precisamos estimar a largura de banda “ideal”.

GWRbanda <- gwr.sel(log(inc + 1) ~ IDHM, data = PE2,
    coords = coordenadas, adapt = T)

## Adaptive q: 0.382 CV score: 381.1 
## Adaptive q: 0.618 CV score: 387.6 
## Adaptive q: 0.2361 CV score: 373.9 
## Adaptive q: 0.1459 CV score: 364.4 
## Adaptive q: 0.09017 CV score: 356 
## Adaptive q: 0.05573 CV score: 355.6 
## Adaptive q: 0.06991 CV score: 356.2 
## Adaptive q: 0.03444 CV score: 364.3 
## Adaptive q: 0.0476 CV score: 357.6 
## Adaptive q: 0.06115 CV score: 355.4 
## Adaptive q: 0.06449 CV score: 355.4 
## Adaptive q: 0.06187 CV score: 355.4 
## Adaptive q: 0.06224 CV score: 355.4 
## Adaptive q: 0.06215 CV score: 355.4 
## Adaptive q: 0.0622 CV score: 355.4 
## Adaptive q: 0.06211 CV score: 355.4 
## Adaptive q: 0.06215 CV score: 355.4

GWRbanda

## [1] 0.06215

Ajustando o modelo:

dengue.gwr <- gwr(log(inc + 1) ~ IDHM, data = PE2,
    coords = coordenadas, adapt = GWRbanda, hatmatrix = TRUE,
    se.fit = TRUE)
dengue.gwr

## Call:
## gwr(formula = log(inc + 1) ~ IDHM, data = PE2, coords = coordenadas, 
##     adapt = GWRbanda, hatmatrix = TRUE, se.fit = TRUE)
## Kernel function: gwr.Gauss 
## Adaptive quantile: 0.06215 (about 11 of 184 data points)
## Summary of GWR coefficient estimates at data points:
##               Min. 1st Qu. Median 3rd Qu.  Max. Global
## X.Intercept. -9.94   -4.51  -1.31    2.25  6.44   0.58
## IDHM         -4.01    3.20   8.74   13.56 22.62   5.47
## Number of data points: 184 
## Effective number of parameters (residual: 2traceS - traceS'S): 23.35 
## Effective degrees of freedom (residual: 2traceS - traceS'S): 160.6 
## Sigma (residual: 2traceS - traceS'S): 1.358 
## Effective number of parameters (model: traceS): 16.77 
## Effective degrees of freedom (model: traceS): 167.2 
## Sigma (model: traceS): 1.331 
## Sigma (ML): 1.269 
## AICc (GWR p. 61, eq 2.33; p. 96, eq. 4.21): 649.4 
## AIC (GWR p. 96, eq. 4.22): 626.6 
## Residual sum of squares: 296.3 
## Quasi-global R2: 0.2549

Colocando a saída do modelo dentro de um dataframe:

results <- as.data.frame(dengue.gwr$SDF)
head(results)

##   sum.w X.Intercept.   IDHM X.Intercept._se IDHM_se   gwr.e  pred pred.se localR2 X.Intercept._se_EDF IDHM_se_EDF pred.se.1      X      Y
## 1 14.88       0.7623  5.026           3.386   5.226 -1.6280 4.175  0.3164 0.07349               3.455       5.332    0.3228 -35.02 -7.879
## 2 14.64      -0.5301  7.412           4.686   7.820  1.3988 4.339  0.5436 0.07056               4.781       7.978    0.5547 -37.62 -7.737
## 3 18.21      -3.7654 12.055           4.043   6.769  0.7456 3.323  0.2477 0.33591               4.125       6.906    0.2528 -41.01 -8.619
## 4 23.48      -3.8527 12.241           3.466   5.954  1.4915 3.394  0.2101 0.22686               3.536       6.075    0.2144 -35.95 -8.457
## 5 19.53      -6.4564 15.852           3.695   6.265  0.5644 2.310  0.3167 0.22434               3.770       6.392    0.3231 -35.48 -8.729
## 6 20.67       5.7463 -2.824           2.921   5.135  1.1993 4.261  0.2978 0.17956               2.980       5.239    0.3039 -37.04 -9.087

Checando os residuos para verificar a presença de autocorrelação para o modelo GWR:

results$residuos <- log(PE2$inc + 1) - results$pred
moran.test(results$residuos, pesos.viz)

## 
##  Moran I test under randomisation
## 
## data:  results$residuos  
## weights: pesos.viz    
## 
## Moran I statistic standard deviate = -0.064, p-value = 0.5
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##         -0.008494         -0.005464          0.002231

Definindo as paletes de cores para a construção dos mapas.

pal <- brewer.pal(n = 10, name = "BuPu")

Mapa dos coeficientes de regressão para a variável IDHM.

PE2$coef.IDH <- results$IDHM
PE2$localR2 <- results$localR2
PE2$pred.gwr <- results$pred

map.idh <- ggplot(PE2) + geom_sf(aes(fill = coef.IDH),
    color = "black", size = 0.2) + scale_fill_gradientn(colours = pal) +
    ggtitle("Distribuição dos coeficientes da variável IDH") +
    theme_void()

Mapeando os coeficientes de regressão para a variável IDHM por significância através do teste de Wald.

# Calculando a estatística de wald
PE2$wald.teste <- abs(results$IDHM/results$IDHM_se)
# Dicotomizando a estatística de wald
PE2$wald.teste <- ifelse(PE2$wald.teste < 2, 0, 1)


map.wald <- ggplot(PE2) + geom_sf(aes(fill = factor(wald.teste)),
    color = "black", size = 0.2) + scale_fill_manual(values = c("white",
    "purple"), labels = c("< 2", ">= 2"), name = "Wald") +
    ggtitle("Coef. IDH significativos") + theme_void()

gridExtra::grid.arrange(map.idh, map.wald, nrow = 2)

Mapa dos coeficientes de determinação regionalizados ( $R^2$ local).

pal <- brewer.pal(n = 10, name = "YlOrRd")
ggplot(PE2) + geom_sf(aes(fill = localR2), color = "black",
    size = 0.2) + scale_fill_gradientn(colours = pal) +
    ggtitle("R² local") + theme_void()

Mapeando os valores observados e preditos dos modelos ajustados.

# quantile(log(PE2$inc+1),seq(0,1,by=0.2))
breaks <- c(-1, quantile(log(PE2$inc + 1), seq(0, 1,
    by = 0.2)))
labels <- c("[0]", "(0-2.6]", "(2.6-3.4]", "(3.4-4.4]",
    "(4.4-5.2]", "> 5.2")

PE2$brks <- cut(log(PE2$inc + 1), include.lowest = FALSE,
    right = TRUE, breaks = breaks, labels = labels)

bruta <- ggplot(PE2) + geom_sf(aes(fill = brks), color = "black",
    size = 0.2) + ggtitle("Log da Incidência") + scale_fill_discrete_sequential(palette = "Heat2",
    c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2, p2 = 1.5,
    na.value = "grey75", drop = FALSE, name = "") +
    theme_void()


PE2$brks.lm <- cut(dengue.lm$fitted.values, lowest = FALSE,
    right = TRUE, breaks = breaks, labels = labels)


pred.lm <- ggplot(PE2) + geom_sf(aes(fill = brks.lm),
    color = "black", size = 0.2) + ggtitle("Predito - LM") +
    scale_fill_discrete_sequential(palette = "Heat2",
        c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2,
        p2 = 1.5, na.value = "grey75", drop = FALSE,
        name = "") + theme_void()

PE2$brks.car <- cut(dengue.car$fit$fitted.values, lowest = FALSE,
    right = TRUE, breaks = breaks, labels = labels)


pred.car <- ggplot(PE2) + geom_sf(aes(fill = brks.car),
    color = "black", size = 0.2) + ggtitle("Predito - CAR") +
    scale_fill_discrete_sequential(palette = "Heat2",
        c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2,
        p2 = 1.5, na.value = "grey75", drop = FALSE,
        name = "") + theme_void()

PE2$brks.gwr <- cut(results$pred, lowest = FALSE, right = TRUE,
    breaks = breaks, labels = labels)


pred.gwr <- ggplot(PE2) + geom_sf(aes(fill = brks.car),
    color = "black", size = 0.2) + ggtitle("Predito - GWR") +
    scale_fill_discrete_sequential(palette = "Heat2",
        c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2,
        p2 = 1.5, na.value = "grey75", drop = FALSE,
        name = "") + theme_void()

ggpubr::ggarrange(bruta, pred.lm, pred.car, pred.gwr,
    ncol = 1, common.legend = TRUE, legend = "right")

Verificando a distribuição dos resíduos:

vioplot::vioplot(dengue.lm$residuals, dengue.car$fi$residuals,
    results$residuos, names = c("LM", "CAR", "GWR"),
    col = "orange")
title("Gráficos violinos da distribuição dos resíduos")

Mapeando os resíduos dos modelos ajustados:

breaks <- c(-4, -1, 0, 1.2, 3.6)
labels <- c("< -1", "-1 a 0", "0 a 1.2", "> 1.2")

PE2$brks.res.lm <- cut(dengue.lm$residuals, include.lowest = FALSE,
    right = TRUE, breaks = breaks, labels = labels)

res.lm <- ggplot(PE2) + geom_sf(aes(fill = brks.res.lm),
    color = "black", size = 0.2) + ggtitle("Resíduos - LM") +
    scale_fill_discrete_sequential(palette = "Purple-Yellow",
        c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2,
        p2 = 1.5, na.value = "grey75", drop = FALSE,
        name = "") + theme_void()

PE2$brks.res.car <- cut(dengue.car$fit$residuals, include.lowest = FALSE,
    right = TRUE, breaks = breaks, labels = labels)

res.car <- ggplot(PE2) + geom_sf(aes(fill = brks.res.car),
    color = "black", size = 0.2) + ggtitle("Resíduos - CAR") +
    scale_fill_discrete_sequential(palette = "Purple-Yellow",
        c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2,
        p2 = 1.5, na.value = "grey75", drop = FALSE,
        name = "") + theme_void()


PE2$brks.res.gwr <- cut(results$residuos, include.lowest = TRUE,
    right = TRUE, breaks = breaks, labels = labels)

res.gwr <- ggplot(PE2) + geom_sf(aes(fill = brks.res.gwr),
    color = "black", size = 0.2) + ggtitle("Resíduos - GWR") +
    scale_fill_discrete_sequential(palette = "Purple-Yellow",
        c1 = 80, c2 = 30, l1 = 30, l2 = 90, p1 = 0.2,
        p2 = 1.5, na.value = "grey75", drop = FALSE,
        name = "") + theme_void()

ggpubr::ggarrange(res.lm, res.car, res.gwr, ncol = 1,
    common.legend = TRUE, legend = "right")

15.7 Obtendo centroides

Aqui vamos utilizar três formas de obter os centroides dos polígonos do estado do Rio de Janeiro como exemplo.

Obtendo o polígono dos municípios do estado do Rio de Janeiro utilizando o pacote geobr.

library(geobr)

RJ <- read_municipality('RJ')

ggplot(RJ) + geom_sf() + theme_void()

Veja o que acontece se tentamos obter o contorno do estado do Rio de Janeiro através da função st_union:

cont <- st_union(RJ)

ggplot(cont) + geom_sf() + theme_void()

Muitos polígonos do IBGE e de outras fontes podem possuir problemas nas geometrias das linhas, de forma que a função st_union não consegue dissolver corretamente os polígonos.

Porém, podemos simplesmente pegar o polígono do estado e usar como contorno:

cont_RJ <- read_state('RJ')

ggplot(cont_RJ) + geom_sf() + theme_void()

Agora vamos obter os centroides.

Método 1

Obtendo o centroides através da função st_centroid. O argumento of_largest_polygon = TRUE força o centroide ficar dentro do maior polígono do município, útil para casos com ilhas.

centroides1 <- st_centroid(RJ, of_largest_polygon = TRUE)

ggplot() + geom_sf(data = cont_RJ) + geom_sf(data = centroides1,
    col = "red") + geom_sf(data = RJ, fill = "transparent") +
    theme_void()

Método 2

Método ponto na superfície.

centro2 <- RJ %>%
    mutate(x = map_dbl(geom, ~st_point_on_surface(.x)[[1]]),
        y = map_dbl(geom, ~st_point_on_surface(.x)[[2]]))

centroides2 <- st_drop_geometry(centro2)[, c("x", "y")]

ggplot() + geom_sf(data = cont_RJ) + geom_point(data = centroides2,
    aes(x = x, y = y), col = "seagreen") + geom_sf(data = RJ,
    fill = "transparent") + theme_void()

Método 3

Usando as coordenadas das sedes do município, disponíveis no arquivo CADMUN:

cadmun <- read_csv2("https://gitlab.procc.fiocruz.br/oswaldo/eco2019/-/raw/master/dados/CADMUN.csv") %>%
    # Filtrando o estado do RJ e linhas sem
    # coordenadas geográficas
filter(UFCOD == "33" & LATITUDE != 0) %>%
    select(MUNCODDV, LATITUDE, LONGITUDE)

ggplot() + geom_sf(data = cont_RJ) + geom_point(data = cadmun,
    aes(x = LONGITUDE, y = LATITUDE), col = "darkblue") +
    geom_sf(data = RJ, fill = "transparent") + theme_void()

Comparando os três métodos

ggplot() + geom_sf(data = cont_RJ) + geom_sf(data = centroides1,
    col = "red", alpha = 0.7) + geom_point(data = centroides2,
    aes(x = x, y = y), col = "seagreen", alpha = 0.7) +
    geom_point(data = cadmun, aes(x = LONGITUDE, y = LATITUDE),
        col = "darkblue", alpha = 0.7) + geom_sf(data = RJ,
    fill = "transparent") + theme_void()

Ainda é possível alterar manualmente as coordenadas dos centroides. Um site útil para isto é: https://www.geoplaner.com/.

15.8 Bibliografia sugerida

Trevor C. Bailey , Anthony C. Gatrell Routledge. Interactive Spatial Data Analysis. 1995.

Lance A. Waller, Carol A. Gotway. Applied Spatial Statistics for Public Health Data. Wiley-Interscience 1St ed. 2004.

Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer. Applied Spatial Data Analysis with R. 2nd ed. 2013

Online:

Geocomputation with R