1 Clusters de Acidentes por Affinity Propagation

library(leaflet.extras)
library(apcluster)  
library(magrittr)
library(dplyr)   
library(leaflet)
library(rgdal)
library(rgeos)
library(geojsonio)
library(mapview)
library(contoureR)
library(geosphere)
library(cluster)
library(sp)
library(rgdal)
library(RColorBrewer)
library(sp)
library(foreign)
library(car)
library(lubridate)
library(tidyr)
library(stringr)

1.1 Carga de Dados

dados = read.dbf("/Users/fagne/OneDrive/r-files/CIET/acidentes2020/_Base/acidentes_2014a2020_WGS84.dbf")
dados = dados[dados$ANO > 2014, ]
sort(unique(dados$ANO))

## [1] 2015 2016 2017 2018 2019 2020

anos = length(unique(dados$ANO))
anos

## [1] 6

class(dados)

## [1] "data.frame"

x2 <- cbind(dados$LONGITUDE, dados$LATITUDE)
x2 <- x2[complete.cases(x2), ]
dim(x2)

## [1] 75102     2

head(x2)

##           [,1]      [,2]
## [1,] -51.01213 -30.19741
## [2,] -51.03732 -30.21369
## [3,] -51.05466 -30.23513
## [4,] -51.05470 -30.23517
## [5,] -51.05470 -30.23517
## [6,] -51.05470 -30.23517

1.2 Preparação

Os dados já foram previamente clusterizados conforme o documento https://rpubs.com/fagnersutel/837844. Dessa maneira, as seçõs de preparação e treino não reprocessam os dados e apenas fazem a carga do modelo previamente obtido:

x1 <- x2

load("data/AZURE/x2-20000-99925.rda")
load("data/AZURE/apres2-20000-99925.rda")
names(x2) = c("LONGITUDE", "LATITUDE" )
head(x2)

dim(x1)

## [1] 75102     2

dim(x2)

## [1] 18776     2

1.3 Treino

plot(apres, x2)

summary(apres)

##   Length    Class     Mode 
##     3639 APResult       S4

1.4 Obtenção de Centróides

centroides = unique(apres@exemplars)
poly = data.frame()
centr_indice = 0
for (i in centroides){
  centr_indice = centr_indice + 1
  centr_lat=x2[i,1]
  centr_lon=x2[i,2]
  poly = rbind(poly, c(centr_lat, centr_lon, centr_indice))
}
names(poly) = c("Lat", "Lon", "Cluster")
#head(poly)
dim(poly)

## [1] 3639    3

exemplars = poly
frame = exemplars

1.5 Classificação Global

predict.apcluster <- function(s, exemplars, newdata){
  simMat <- s(rbind(exemplars, newdata), sel=(1:nrow(newdata)) + nrow(exemplars))[1:nrow(exemplars), ]
  unname(apply(simMat, 2, which.max))
}
resultado <- list()
dados$cluster = 0
for(i in seq(from=1, to=length(dados$ID)-1000, by=1000)){
  inicio = i
  final = i+999
  resultado = predict.apcluster(negDistMat(r=2), x2[apres@exemplars, ],  dados[inicio:final, 2:3])
  dados$cluster[inicio:final] = resultado
}
controle = length(dados$cluster)  - final
resultado = predict.apcluster(negDistMat(r=2), x2[apres@exemplars, ],  dados[(final + 1):length(dados$cluster), 2:3])
dados$cluster[(final + 1):length(dados$cluster)] = resultado
#head(dados)
#tail(dados)

1.6 Acidentes por Cluster

pal <- colorFactor(
  palette = 'Dark2',
  domain = dados$cluster
)
leaflet(dados) %>%
   addTiles(group="Mapa") %>% 
   addCircles(group="Acidentes", ~LONGITUDE, ~LATITUDE, weight = 0.1, radius=7, color=~pal(cluster),
              stroke = TRUE, fillOpacity = 0.8, popup=~paste("Cluster Nº: ", cluster,  
             "<br>Ano: ", ANO, "<br>Tipo: ", TIPO_ACID, "<br>Local: ", LOG1,  "<br>UPS: ", UPS,   sep = " ")) %>% 
   addLegend(group="Legenda", "topright", colors= "", labels=paste("Classified into ", summary(apres)[1], "clusters"), title="Accidents in Porto Alegre") %>% 
   addProviderTiles(providers$CartoDB)

1.7 Enriquecimento Informacional

Nessa etapa são calculados atributos como área de polígonos, acidentes por cluster

clusters_encontrados = sort(unique(dados$cluster))
parq = dados
poly = data.frame()
total = 0
for (i in clusters_encontrados){
  temp = parq[parq$"cluster" == i,  ]
  ch1 = convexHullAM_Indexes(temp[,2],temp[,3], includeColinear=FALSE,zeroBased = FALSE)
  poligono = temp[ch1, 2:3 ]
  area <- geosphere::areaPolygon(x = poligono)
  acidentes = nrow(temp)
  pol = temp
  coordinates(pol) = ~LONGITUDE+LATITUDE
  centr_lat=gCentroid(pol, byid=FALSE)$x
  centr_lon=gCentroid(pol, byid=FALSE)$y
  total = total+1
  if(nrow(temp) >= anos) {
    for (ii in ch1) {
    polying = temp[ii,]
    polying$area = area * 10^-6
    polying$acidentes = acidentes
    polying$centroide_lat = centr_lat
    polying$centroide_lon = centr_lon
    poly = rbind(poly, polying)
    }  
  }
}
total

## [1] 3603

length(unique(poly$cluster))

## [1] 2920

mean(poly$area)

## [1] 0.0151733

median(poly$area)

## [1] 0.004853773

minimoquantil = quantile(poly$area, probs = 0.01)
maximoquantil = quantile(poly$area, probs = 0.90)
quantile(poly$area, probs = c(0.01, 0.25, 0.5,0.75,0.99))

##           1%          25%          50%          75%          99% 
## 8.258996e-05 1.489416e-03 4.853773e-03 9.714150e-03 9.975409e-02

#poly = poly[(poly$area < maximoquantil) & (poly$area > minimoquantil), ]
dim(poly)

## [1] 18097    48

class(poly)

## [1] "data.frame"

pol = poly
save(poly, file = "polyoriginal.Rda")

1.8 Combinando Informações

Nessa etapa os dados de acidentes contidos no data set original são combinados com dados de clusters

dados = poly[,c(1:9, 11:13,38,41,44:48)]
names(dados) = c("ID","lat", "lon", "log1", "Log2", "Pred", "Local", "Tipo", "Via", "Data", "Dia", "Hora", 
                 "Fx_horaria","UPS", "box_id", "Area", "Acidentes", "CentLon", "CentLat")
dados$id = (dados$box_id * 11)
dados$group = dados$id
head(dados)

nrow(dados)

## [1] 18097

dadostemp = dados[, c(15:21)]
coordinates(dados)=c("lat","lon")
df = dados
length(df)

## [1] 18097

data <- data.frame(box_id=unique(df$box_id),row.names=unique(df$id))
nrow(data)

## [1] 2920

dadostemp2 = dados[!duplicated(dados$id),]
data = as.data.frame(cbind(data, dadostemp2@data))

1.9 Criação de Polígonos

São criados polígonos delimitados pelos limites externos dos clusters.

points2polygons <- function(df,data) {
  get.grpPoly <- function(group,ID,df) {
    Polygon(coordinates(df[df$id==ID & df$group==group,]))
  }
  get.spPoly  <- function(ID,df) {
    Polygons(lapply(unique(df[df$id==ID,]$group),get.grpPoly,ID,df),ID)
  }
  spPolygons  <- SpatialPolygons(lapply(unique(df$id),get.spPoly,df))
  SpatialPolygonsDataFrame(spPolygons,match.ID=T,data=data)
}
data$Log2 = NULL
spDF <- points2polygons(df,data)
spDF

## class       : SpatialPolygonsDataFrame 
## features    : 2920 
## extent      : -51.27422, -51.05466, -30.24464, -29.96328  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 19
## names       : box_id,     ID,                                    log1,  Pred,      Local,         Tipo,                                                 Via,       Data,         Dia,  Hora, Fx_horaria, UPS, box_id.1,             Area, Acidentes, ... 
## min values  :      1, 600973,                 AC B VILA NOVA BRASILIA,     0, Cruzamento, ABALROAMENTO,                                    0 R CEL CLAUDINO, 01/01/2015,     DOMINGO, 00:00,          0,   1,        1,                0,         6, ... 
## max values  :   3639, 683107, VDT COMPLEXO VIARIO TELMO T FLORES-AC B, 21010, Logradouro,   TOMBAMENTO, VDT ABDIAS DO NASCIMENTO & AV EDVALDO PEREIRA PAIVA, 31/12/2020, TERCA-FEIRA, 23:57,         23,  13,     3639, 3.97512935538404,       230, ...

length(spDF)

## [1] 2920

class(spDF)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

spDF@data$group = 1
spDF@data$box_id = NULL
dim(spDF@data)

## [1] 2920   18

dadostemp = unique(dadostemp)
spDF@data = merge(spDF@data, dadostemp, by = "box_id")
dim(spDF@data)

## [1] 2920   24

plot(spDF,col=spDF$box_id+1)

library(rgdal)
rgdal::writeOGR(obj = spDF,
                dsn = "data/myParq999-25.json",
                layer = "myParq",
                driver = "GeoJSON",
                overwrite_layer = TRUE)

Acidentes por Cluster

#carregamos os dados SpatialPolygonsDataFrame
parqs <- geojsonio::geojson_read("data/myParq999-25.json", what = "sp")
dim(parqs)

## [1] 2920   24

library(raster)
projection(parqs)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

library(mapview)
mapviewPalette(name = "Viridis")
library(RColorBrewer)
mapview(parqs, zcol = "Acidentes.x", col.regions=brewer.pal(9, "YlOrRd"))

Acidentes por m2

Acidentes por poligono

hist(parqs@data$Acidentes.x, col = "magenta")

quantile(parqs@data$Acidentes.x,  probs = c(0.1, 0.25, 0.95, 0.99, 1))

##  10%  25%  95%  99% 100% 
##    8   11   64  101  230

Acidentes por KM2

hist(parqs@data$densidade, col = "orange")

Areas

options(scipen = 999)
quantile(parqs@data$Area.x,  probs = c(0.1, 0.25, 0.95, 1))

##          10%          25%          95%         100% 
## 0.0004522779 0.0015016362 0.0261328007 3.9751293554

Locais mais densos

quantile(parqs$densidade, probs = 0.99)

##          99% 
## 296528828633

temp = parqs[parqs$densidade > quantile(parqs$densidade, probs = 0.99), ]
#mapview(temp, zcol = "densidade")

1.10 Projeção de Dados

OS dados são reprojetados para Mercator afim de permitir que se calcule a distancia entre centróidos de clusters em metros.

frame = frame[, 1:2]
names(frame) = c("lon", "lat")
library(fossil)
distM <- as.matrix( fossil::earth.dist(frame))
a = as.numeric(unlist(lapply(1:nrow(distM), function(x) mean(distM[x, order(distM[x, ])[2:11]]))))
round(mean(a)*1000,0)

## [1] 318

length(a)

## [1] 3639

projection(parqs)

## [1] "+proj=longlat +datum=WGS84 +no_defs"

parqstemp = parqs
require(sf)
shape <- read_sf(dsn = ".", layer = "mercator_32722_2014_2019")
projection(shape)

## [1] "+proj=utm +zone=22 +south +datum=WGS84 +units=m +no_defs"

projection(parqstemp) = projection(shape)

ccods = coordinates(parqs)
temps = as.data.frame(ccods)
cord.dec = SpatialPoints(cbind(temps$V1, temps$V2), proj4string=CRS("+proj=longlat"))
cord.UTM <- spTransform(cord.dec, CRS("+init=epsg:32722"))
ccods = as.data.frame(cord.UTM)
points = cbind(ccods[,1],ccods[,2])
head(points)

##          [,1]    [,2]
## [1,] 480682.3 6680702
## [2,] 480644.7 6675696
## [3,] 480055.4 6677661
## [4,] 482067.8 6680599
## [5,] 480547.1 6677103
## [6,] 477890.7 6676124

library(spdep)
distNeighbors = round(mean(a)*1000,0)  #318 #400
dnb = dnearneigh(points,0,distNeighbors)
class(dnb)

## [1] "nb"

subsets = as.data.frame(matrix(dnb))
class(subsets)

## [1] "data.frame"

subsets = subsets$V1
length(subsets)

## [1] 2920

parqs$n = 1
sub = which(subsets == '0')
sub

##   [1]   21   37   45   54   70  119  131  139  162  207  214  231  232  236  313
##  [16]  366  440  458  467  475  499  515  531  532  544  546  623  632  702  712
##  [31]  726  729  735  762  786  873  891  962  988 1027 1048 1057 1066 1073 1084
##  [46] 1132 1134 1161 1173 1197 1213 1264 1294 1317 1353 1381 1438 1441 1448 1483
##  [61] 1504 1508 1551 1561 1597 1647 1656 1667 1670 1677 1696 1716 1740 1777 1779
##  [76] 1804 1814 1815 1900 1908 1919 1922 1985 1988 1993 2001 2065 2079 2091 2100
##  [91] 2147 2215 2230 2241 2257 2288 2291 2302 2306 2364 2369 2398 2411 2419 2434
## [106] 2436 2485 2508 2569 2609 2614 2636 2645 2652 2701 2740 2742 2822 2823 2837
## [121] 2898 2899 2914 2915

parqs$n[sub] = 0
length(parqs)

## [1] 2920

parqs = parqs[parqs$n > 0,]
length(parqs)

## [1] 2796

length(dnb)

## [1] 2920

ccods = ccods[-sub, ]
dim(ccods)

## [1] 2796    2

points = cbind(ccods[,1],ccods[,2])
dnb = dnearneigh(points,0,distNeighbors)
dnb

## Neighbour list object:
## Number of regions: 2796 
## Number of nonzero links: 21156 
## Percentage nonzero weights: 0.2706196 
## Average number of links: 7.566524

length(dnb)

## [1] 2796

1.11 Matriz de Vizinhanca

É criada a matriz de vizinhança de clusters para calcular suas associações espaciais.

1.11.1 Matriz Binária

W.Bin= nb2mat(neighbours = dnb, style = "B")

1.11.2 Matriz Normalizada

W.Normal= nb2mat(neighbours = dnb, style = "W")

1.12 KNN

vizinhos_4 <- knearneigh(points, k = 4)
class(vizinhos_4)

## [1] "knn"

head(vizinhos_4$nn)

##      [,1] [,2] [,3] [,4]
## [1,]  962  378 2336  994
## [2,]  631 1585  668 1930
## [3,] 1451  495 2434 1264
## [4,] 1119 1042 2057  272
## [5,] 1606 2389 1974 1510
## [6,]  168 2745 1728 1833

vizinhanca_4 <- knn2nb(vizinhos_4)
class(vizinhanca_4)

## [1] "nb"

1.13 Preparação para Análise Global e Local

mv_simpl = st_as_sf(parqs)
#plot(mv_simpl)
class(mv_simpl)

## [1] "sf"         "data.frame"

library(dplyr)
mv_simpl =  mv_simpl %>% dplyr::select(Acidentes.y)

class(mv_simpl)

## [1] "sf"         "data.frame"

mapview::mapview(mv_simpl)

sf::sf_use_s2(FALSE)#trips and tiks

## Spherical geometry (s2) switched off

mv_simpl = st_as_sf(mv_simpl)
vizinhanca_neig <- poly2nb(mv_simpl)
ShapeNEIG = parqs
ShapeNEIG$vizinhos = card(vizinhanca_neig)
ShapeNEIG <- subset(ShapeNEIG, parqs$vizinhos != 0)

1.14 Calculando o Índice de Moran Global

Os índices de autocorrelção espacial global calculados pelos testes de normalidade e permutação.

1.14.1 Pelo teste de Normalidade

moran.test(parqs$Acidentes.y,listw=nb2listw(dnb, style = "W"), randomisation= FALSE)

## 
##  Moran I test under normality
## 
## data:  parqs$Acidentes.y  
## weights: nb2listw(dnb, style = "W")    
## 
## Moran I statistic standard deviate = 14.716, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1803498726     -0.0003577818      0.0001507801

1.14.2 Pelo teste de Permutação ou Teste de pseudo-significˆancia

 moran.test(parqs$Acidentes.y,listw=nb2listw(dnb, style = "W"), randomisation= TRUE)

## 
##  Moran I test under randomisation
## 
## data:  parqs$Acidentes.y  
## weights: nb2listw(dnb, style = "W")    
## 
## Moran I statistic standard deviate = 14.763, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1803498726     -0.0003577818      0.0001498236

1.14.3 Por simulação de Monte-Carlo

moran.mc(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"), nsim=999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W")  
## number of simulations + 1: 1000 
## 
## statistic = 0.18035, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

1.14.4 Pelo teste de Permutação

Diferente dos demais testes globais o teste para o EBI é exclusivo para taxas e tem-se apenas a opção de teste da permutação

#EBImoran.mc(parqs$Acidentes.y,parqs$Area.y,
 #           nb2listw(dnb, style="B", zero.policy=TRUE), nsim=999, zero.policy=TRUE)

1.14.5 Por simulação de Monte-Carlo

shapeCG.p=parqs$Acidentes.y/parqs$Area.y
#moran.mc(shapeCG.p, nb2listw(dnb, style="B", zero.policy=TRUE),
#         nsim=999, zero.policy=TRUE)

1.15 Calculando a Estatística C de Geary Global

1.15.1 Pelo teste de Normalidade

geary.test(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"), randomisation= FALSE)

## 
##  Geary C test under normality
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W") 
## 
## Geary C statistic standard deviate = 12.658, p-value <
## 0.00000000000000022
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.8383709118      1.0000000000      0.0001630341

1.15.2 Pelo teste de Permutação

geary.test(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"), randomisation=TRUE)

## 
##  Geary C test under randomisation
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W") 
## 
## Geary C statistic standard deviate = 9.8068, p-value <
## 0.00000000000000022
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.8383709118      1.0000000000      0.0002716369

1.15.3 Por simulação de Monte-Carlo

geary.mc(parqs$Acidentes.y, listw=nb2listw(dnb, style = "W"),nsim=999)

## 
##  Monte-Carlo simulation of Geary C
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "W") 
## number of simulations + 1: 1000 
## 
## statistic = 0.83837, observed rank = 1, p-value = 0.001
## alternative hypothesis: greater

1.16 Calculando Índice de Getis e Ord Global

Getis-Ord é um indicador que mede a concentração local de uma variável de atributo distribuída espacialmente

globalG.test(parqs$Acidentes.y, nb2listw(dnb, style="B"))

## 
##  Getis-Ord global G statistic
## 
## data:  parqs$Acidentes.y 
## weights: nb2listw(dnb, style = "B") 
## 
## standard deviate = 20.789, p-value < 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##  0.003831266871821  0.002707164080555  0.000000002923849

1.17 Getis e Ord Local

localG(parqs$Acidentes.y, nb2listw(dnb, style="B"), zero.policy=NULL, spChk=NULL, return_internals=FALSE)[1:50]

##  [1]  1.20567047 -1.17236886  0.29069182  0.57129456  1.07015366  3.88489654
##  [7]  0.91119347  2.69173479  0.02341446 -0.60222652  0.89671505  1.50199457
## [13] -0.87392495  2.15133971 -1.41091289 -0.56588860 -0.98421860 -0.20538193
## [19]  1.81277911 -0.39462222  0.77441828 -0.87195422  0.35958863 -0.69840062
## [25] -0.48536819  0.67982984 -0.75169103  1.47100259  3.42517188 -0.24570449
## [31]  1.46896849  7.58442241  0.14615659 -1.33404304  0.13050851 -1.63473713
## [37] -0.93653276 -1.37241645  0.62667345 -0.37019731  0.35957715  0.38261616
## [43]  1.53682970  1.82531545  1.63305841  0.21829896  2.66991973  4.56722261
## [49] -0.78337975  1.09555459

1.18 Moran Local

Todas as análises feitas até o momento foram de escala global. No entanto, é necessário que seja feita também uma análise local do estudo. Essa análise pode ser feita pelo índice local de autocorrelaçãoo espacial (LISA). Para isso é preciso calcular o índice de Moran local.

ShapePB.mloc <- localmoran(parqs$Acidentes.y, listw=nb2listw(dnb, style="W")) 
head(ShapePB.mloc, 10)

##              Ii           E.Ii      Var.Ii        Z.Ii Pr(z != E(Ii))
## 1   0.853795286 -0.00309849209 0.505122518  1.20567047   0.2279445748
## 2   0.305858094 -0.00024466748 0.068171986  1.17236886   0.2410489969
## 3  -0.077455853 -0.00030361317 0.070441964 -0.29069182   0.7712870280
## 4  -0.170942033 -0.00019207747 0.089330760 -0.57129456   0.5677999911
## 5  -0.076385746 -0.00001644005 0.005092668 -1.07015366   0.2845501473
## 6   0.355458144 -0.00003306175 0.008373342  3.88489654   0.0001023735
## 7  -0.105686234 -0.00007244147 0.013434456 -0.91119347   0.3621934429
## 8  -0.344749487 -0.00008840850 0.016395322 -2.69173479   0.0071081446
## 9   0.005516032 -0.00021500881 0.059909933  0.02341446   0.9813196728
## 10 -0.305436991 -0.00064308250 0.256149086 -0.60222652   0.5470233634

write.csv2(ShapePB.mloc, file = "ShapePBmloc.csv")

1.19 Mapa das probabilidades (Signific?ncias do I de Moral Local)

Por meio dos valor-p do éndice de Moran local é possível construir um mapa de probabilidades.

library(classInt)
INT4 <- classIntervals(ShapePB.mloc[,5], style="fixed", 
                       fixedBreaks=c(0,0.01, 0.05, 0.10))
CORES.4 <- c(rev(brewer.pal(3, "Reds")), brewer.pal(3, "Blues"))
COL4 <- findColours(INT4, CORES.4)
parqs$COL = COL4  
parqs$p_valor = ifelse(parqs$COL == "#DE2D26", "[0,0.01)", ifelse(parqs$COL == "#EEE5E4", "[0.01,0.05)", "[0.05,0.1]"))
plot(parqs, col=COL4)
title("P-valores do I de Moran Local por Distäncia de Centróides")
TB4 <- attr(COL4, "table")
legtext <- paste(names(TB4))
legend("bottomright", fill=attr(COL4, "palette"), legend=legtext, 
       bty="n", cex=0.7, y.inter=0.7)

APLs = parqs
colnames(APLs@data)[28] = "P-value"
mapview(APLs, zcol = "P-value", col.regions=c("red", "orange", "gray"))

temp = parqs[parqs$p_valor != "[0.05,0.1]", ]
mapview(temp, zcol = "p_valor", col.regions=c("red", "orange"))

1.19.1 Montando matrix W de vizinhança

ShapeCG.nb1.mat <- nb2mat(dnb)

1.19.2 Incidência de acidentes padronizada

Acidentes_SD <- scale(parqs$Acidentes.y)

1.19.3 Média das incidências de acedentes padronizada

Acidentes_W <- ShapeCG.nb1.mat %*% Acidentes_SD

2 Diagrama de espalhamento de Moran

plot(Acidentes_SD, Acidentes_W,xlab="Normalized Values (Z)",ylab="Average of Normalized Neighbors (WZ)")
abline(v=0, h=0)
title("Diagrama de Espalhamento de Moran por Distancia de Centróides")

Q <- vector(mode = "numeric", length = nrow(ShapePB.mloc))
Q[(Acidentes_SD>0  & Acidentes_W > 0)] <- 1            
Q[(Acidentes_SD<0  & Acidentes_W < 0)] <- 2
Q[(Acidentes_SD>=0 & Acidentes_W < 0)] <- 3
Q[(Acidentes_SD<0  & Acidentes_W >= 0)]<- 4
signif=0.05
parqs$Q = Q

2.0.0.1 Quadrantes

table(Q)

## Q
##    1    2    3    4 
##  631 1134  364  667

3 Mapa LISA

Q[ShapePB.mloc[,5]>signif]<-5
CORES.5 <- c("blue", "green" , "red", "yellow", "gray", rgb(0.95,0.95,0.95))
#CORES.5 <- c(1:5, rgb(0.95,0.95,0.95))
parqs$cores5Q = CORES.5[Q]
plot(parqs, col=CORES.5[Q])
title("Mapa LISA por Distancia Centroides")
legend("bottomright", c("Q1(+/+)", "Q2(-/-)", "Q3(+/-)", "Q4(-/+)","NS"), 
       fill=CORES.5)

CORES.5[Q][1:5]

## [1] "gray" "gray" "gray" "gray" "gray"

head(CORES.5[Q])

## [1] "gray" "gray" "gray" "gray" "gray" "blue"

parqs$cores5 =  ifelse(parqs$cores5Q == "blue", "H-H", ifelse(parqs$cores5Q == "green", "L-L", 
            ifelse(parqs$cores5Q == "red", "H-L", ifelse(parqs$cores5Q == "yellow", "L-H", "N-S"))))
APLs = parqs
colnames(APLs@data)[31] = "Relationships"
mapview(APLs, zcol = "Relationships", col.regions=c("red", "orange", "green", "yellow", "grey"))

APLs = parqs[parqs$cores5 == "H-H", ]
colnames(APLs@data)[31] = "Relationships"
mapview(APLs, zcol = "Relationships", col.regions=c("red", "orange", "green", "yellow", "grey"))

mapview(APLs, zcol = "Relationships", col.regions=c("red", "orange", "green", "yellow", "grey"))

mapview(temp, zcol = "p_valor", col.regions=c("red", "blue"))

4 Análise longitudinal

aa = parq[parq$cluster == 2699,]
aa = parq
library(lubridate)
aa$date = dmy(aa$DATA)
library(changepoint)
library(tidyverse)

df = aa %>%
  group_by(ANO, MES) %>%
  summarise(Total=n(),
            UPS=sum(UPS)) 
head(df)

4.1 Variação Anual

dfanual = aa %>%
  group_by(ANO) %>%
  summarise(Total=n(),
            UPS=sum(UPS)) 

dfanual = dfanual %>% fill(ANO, .direction = "down") %>%
  mutate(VariacaoAbsoluta = as.numeric(formatC((Total - lag(Total))*100/lag(Total), digits = 2)),
         VariacaoUPS = as.numeric(formatC((UPS - lag(UPS))*100/lag(UPS), digits = 2)))
dfanual

dfanual <- dfanual %>% 
  group_by(ANO) %>% 
  mutate(
    VariacaoAcumuladaAcidentes = cumprod((VariacaoAbsoluta/100)+1)-1,
    VariacaoAcumuladaUPS = cumprod((UPS/100)+1)-1
  )
dfanual

4.2 Variação Acumulada

glimpse(dfanual)

## Rows: 6
## Columns: 7
## Groups: ANO [6]
## $ ANO                        <int> 2015, 2016, 2017, 2018, 2019, 2020
## $ Total                      <int> 19641, 12661, 11524, 10073, 12385, 8818
## $ UPS                        <int> 41265, 30893, 27416, 23881, 30501, 23762
## $ VariacaoAbsoluta           <dbl> NA, -36, -9, -13, 23, -29
## $ VariacaoUPS                <dbl> NA, -25, -11, -13, 28, -22
## $ VariacaoAcumuladaAcidentes <dbl> NA, -0.36, -0.09, -0.13, 0.23, -0.29
## $ VariacaoAcumuladaUPS       <dbl> 412.65, 308.93, 274.16, 238.81, 305.01, 237~

dfanual = dfanual %>% fill(ANO, .direction = "down") %>%
  mutate(Variacao = formatC((Total - lag(Total))*100/lag(Total), digits = 2))

dfanual

df$date = as.Date(paste(df$ANO, df$MES, "01", sep = "-"))
df$date = as.character(df$date)
df$ANO = NULL
df$MES = NULL
df$UPS = NULL
data_ini = as.Date(paste(year(min(df$date)), "01-01",sep = "-"))
data_fim = as.Date(paste(year(max(df$date)), "12-01",sep = "-"))
df = df %>%
  mutate(date = as.Date(date), 
         date = as.Date(format(date, "%Y-%m-01"))) %>%
  tidyr::complete(date = seq(data_ini, data_fim, "1 month"))
df$Total =  df$Total %>% replace_na(0)

df$date = as.Date(df$date)
passageiros = ts(df$Total, start=c(min(df$date)), frequency=12)
tamanho = length(passageiros)
passageiros = as.data.frame(t(matrix(passageiros, 12)))
maximo = (nrow(passageiros)*ncol(passageiros))
start = maximo-tamanho
start = ncol(passageiros)-start
passageiros[nrow(passageiros), start:ncol(passageiros)] = 0
names(passageiros) = paste(c("jan", "fev", "mar", "abr", "mai", "jun", "jul", "ago", "set", "out", "nov", "dez"), rep("01", 12), sep = "-")
passageiros$ano = year(data_ini):year(data_fim)
passageiros

library(reshape2)
test_data_long <- melt(passageiros, id="ano")  
test_data_long$date <- paste(test_data_long$ano, test_data_long$variable, sep = "-")
test_data_long$variable <- as.Date(parse_date(test_data_long$date,"%Y-%b-%d",locale=locale("pt")))
test_data_long <- test_data_long[order(test_data_long$variable),]
test_data_long$ano = NULL
test_data_long$date = NULL
test_data_long = test_data_long[!is.na(test_data_long$value), ]
test_data_long$value = as.numeric(test_data_long$value)
test_data_long$value = round(test_data_long$value, 0)
head(test_data_long)

library(data.table)
require(scales)
p = ggplot(test_data_long, aes(x=variable, y=value))+
  geom_line(size=.8) + 
  scale_x_date(breaks = date_breaks("6 months"),
               labels = date_format("%Y/%m"))+
  theme(axis.text.x=element_text(angle=45, hjust=1),
        plot.title = element_text(size=10, face='bold'))+
  labs(x='Ano', y='Acidentes (Mês)',
       title='Acidentes',
       caption='Fonte: EPTC')
p

4.3 Análise Temporal

passageiros = as.ts(read.zoo(test_data_long, FUN = as.yearmon))
plot(decompose(passageiros, type = "m"))

4.4 Change Point Pattern

nomes = c("x", "y")
df = test_data_long
names(df) = nomes

fit_changepointBINSEG = cpt.mean(df$y, method = "BinSeg", Q=5)

fit_changepointPELT = cpt.mean(df$y, method = "PELT")

fit_changepointAMOC = cpt.mean(df$y, method = "AMOC")

fit_changepointSEGN = cpt.mean(df$y, method = "SegNeigh", pen.value = 0.05, penalty = "AIC")

source("metodoChangePoint.R", local = knitr::knit_global())
#apls = temp$box_id
apls = APLs$box_id
aplssub = apls
length(APLs)

## [1] 174

length(aplssub)

## [1] 174

source("metodoChangePoint.R", local = knitr::knit_global())
seriesAPLs = data.frame()
par(mfrow=c(1,1))
for (i in 1:length(aplssub)) {
  aAPL = metodoChangePoint(aplssub[i], "PELT", pen.valor =  0.5, penalte =   "AIC")
  pontual = aAPL@dataframelocal
  pontual$id = aplssub[i]
  seriesAPLs = rbind(seriesAPLs, pontual)
  df = aAPL@dataframelocal
  if (aAPL@cptns == 1) {
    pontos = aAPL@cptRes2@cpts
    dfstart = df[1:(aAPL@cptRes2@cpts/2), ]
    dfend = df[((aAPL@cptRes2@cpts/2)+1):aAPL@cptRes2@cpts, ]
    dfstart$periodo = "Start"
    dfend$periodo = "End"
    juncao = rbind(dfstart, dfend)
    anova_a <- aov(y ~ periodo, data = juncao)
    aov = as.numeric(unlist(summary(anova_a)[[1]])[9])
    medias <- with(juncao,tapply(y, periodo, mean)) 
    shapiro=  as.numeric(shapiro.test(resid(anova_a))[[2]]) #P-valor > 0.05 os dadosapresentam Distribuição Nornal
    levene = as.numeric(leveneTest(resid(anova_a)~periodo, juncao, center=mean)[[3]][[1]]) # P-valor > 0.05 as varianias dos grupos sao homogeneas
    mannwhitney = wilcox.test(y ~ periodo, data = juncao)$p.value # P-valor > 0.05 as médianas são iguais
  } else if (aAPL@cptns == 2) {
    pontos = aAPL@cptRes2@cpts
    dfstart = df[1:pontos[1],]
    dfend = df[(pontos[1]+1):pontos[length(pontos)],]
    dfstart$periodo = "Start"
    dfend$periodo = "End"
    juncao = rbind(dfstart, dfend)
    anova_a <- aov(y ~ periodo, data = juncao)
    aov = as.numeric(unlist(summary(anova_a)[[1]])[9])
    medias <- with(juncao,tapply(y, periodo, mean)) 
    shapiro=  as.numeric(shapiro.test(resid(anova_a))[[2]]) #P-valor > 0.05 os dadosapresentam Distribuição Nornal
    levene = as.numeric(leveneTest(resid(anova_a)~periodo, juncao, center=mean)[[3]][[1]]) # P-valor > 0.05 as varianias dos grupos sao homogeneas
    mannwhitney = wilcox.test(y ~ periodo, data = juncao)$p.value # P-valor > 0.05 as médianas são iguais
  } else if (aAPL@cptns > 2) {
    pontos = aAPL@cptRes2@cpts
    dfstart = df[1:pontos[1],]
    dfend = df[(pontos[length(pontos)-1]+1):pontos[length(pontos)],]
    dfstart$periodo = "Start"
    dfend$periodo = "End"
    juncao = rbind(dfstart, dfend)
    anova_a <- aov(y ~ periodo, data = juncao)
    aov = as.numeric(unlist(summary(anova_a)[[1]])[9])
    medias <- with(juncao,tapply(y, periodo, mean)) 
    shapiro=  as.numeric(shapiro.test(resid(anova_a))[[2]])
    levene = as.numeric(leveneTest(resid(anova_a)~periodo, juncao, center=mean)[[3]][[1]])
    mannwhitney = wilcox.test(y ~ periodo, data = juncao)$p.value
  }else{
    aov = 99991
    shapiro=  9991
    levene = 9991
    mannwhitney = 9991
  }
  estac = ts(df$y, frequency = 12, start = c(2015, 01))
  estac = Box.test(estac,type="Ljung-Box")[[3]]
  df1 <- data.frame(Cluster = aplssub[i], 
                    Acidentes = aAPL@acidentes, 
                    VarTotal=aAPL@vartotal,
                    VarUPS=aAPL@varups, 
                    VarCPT=aAPL@varcpt, 
                    AOV=aov, 
                    Levene=levene, 
                    Shappiro=shapiro, 
                    MannWhitney = mannwhitney,
                    Pontos = length(pontos))
  resumo = rbind(resumo, df1)
    if (i < 6) {
      plot(aAPL@cptRes2, main = paste("ID: ", aplssub[i], ", Freq: ", aAPL@acidentes, ", Acid.: ", round(aAPL@vartotal, 2), ", UPS: ", round(aAPL@varups,2), ", CPT: ", round(aAPL@varcpt,2), ", AOV: ", round(aov,6), ", MH: ", round(mannwhitney,6), sep = ""), ylab="Frequency", xlab="Month")    
    }
   # source("plotarIndividuais.R", local = knitr::knit_global()) # Plotar graficos por cluster
  
}

4.4.1 Tabela Resumo

head(resumo)

nrow(resumo)

## [1] 174

length(aplssub)

## [1] 174

length(APLs)

## [1] 174

4.4.1.1 Resultadosde ANOVA e Mann-Whitney

resumo$naov = ifelse(resumo$AOV >= 0.05, "Não há Diferença entre médias", "Há Diferença entre médias")
resumo$nmw = ifelse(resumo$MannWhitney >= 0.05, "Não há Diferença entre médianas", "Há Diferença entre médianas")
estacionarios = resumo[is.na(resumo$VarCPT), ]

Nestacionarios = resumo[!is.na(resumo$VarCPT), ]
nrow(estacionarios)

## [1] 138

nrow(Nestacionarios)

## [1] 36

Foram encontrados entre os APLs 138 locais com média estacioária e 36 não estacionários.

4.4.2 Estacionários

4.4.2.1 ANOVA

table(estacionarios$naov)

## 
##     Há Diferença entre médias Não há Diferença entre médias 
##                            29                           109

table(estacionarios$naov)[2]/nrow(estacionarios)

## Não há Diferença entre médias 
##                     0.7898551

4.4.2.2 Mann Whitney

table(estacionarios$nmw)

## 
##     Há Diferença entre médianas Não há Diferença entre médianas 
##                              36                             102

table(estacionarios$nmw)[2]/nrow(estacionarios)

## Não há Diferença entre médianas 
##                       0.7391304

4.4.3 Não Estacionários

4.4.3.1 ANOVA

table(Nestacionarios$naov)

## 
##     Há Diferença entre médias Não há Diferença entre médias 
##                            33                             3

table(Nestacionarios$naov)[1]/nrow(Nestacionarios)

## Há Diferença entre médias 
##                 0.9166667

4.4.3.2 Mann Whitney

table(Nestacionarios$nmw)

## 
##     Há Diferença entre médianas Não há Diferença entre médianas 
##                              33                               3

table(Nestacionarios$nmw)[1]/nrow(Nestacionarios)

## Há Diferença entre médianas 
##                   0.9166667

4.4.4 Tabelas de Contingencia

4.4.4.1 ANOVA

resumo$desfechoaov = ifelse((is.na(resumo$VarCPT)&(resumo$naov == "Não há Diferença entre médias")), 0,1)
resumo$desfechomw = ifelse(!is.na(resumo$VarCPT)&(resumo$nmw == "Há Diferença entre médianas"), 1,0)
resumo$predito = ifelse(is.na(resumo$VarCPT), 0,1)
previsoes = resumo$predito
aov = resumo$desfechoaov 
amw = resumo$desfechomw
library(caret)

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

cm = confusionMatrix(table(previsoes, aov, dnn=c("Predito", "Atual")))
cm

## Confusion Matrix and Statistics
## 
##        Atual
## Predito   0   1
##       0 109  29
##       1   0  36
##                                           
##                Accuracy : 0.8333          
##                  95% CI : (0.7695, 0.8854)
##     No Information Rate : 0.6264          
##     P-Value [Acc > NIR] : 0.000000001876  
##                                           
##                   Kappa : 0.6087          
##                                           
##  Mcnemar's Test P-Value : 0.000000199858  
##                                           
##             Sensitivity : 1.0000          
##             Specificity : 0.5538          
##          Pos Pred Value : 0.7899          
##          Neg Pred Value : 1.0000          
##              Prevalence : 0.6264          
##          Detection Rate : 0.6264          
##    Detection Prevalence : 0.7931          
##       Balanced Accuracy : 0.7769          
##                                           
##        'Positive' Class : 0               
## 

source('matriz_confusao.R')
matriz(cm)

4.4.4.2 Mann-Whitney

cm = confusionMatrix(table(previsoes, amw, dnn=c("Predito", "Atual")))
cm

## Confusion Matrix and Statistics
## 
##        Atual
## Predito   0   1
##       0 138   0
##       1   3  33
##                                            
##                Accuracy : 0.9828           
##                  95% CI : (0.9504, 0.9964) 
##     No Information Rate : 0.8103           
##     P-Value [Acc > NIR] : 0.000000000001532
##                                            
##                   Kappa : 0.9458           
##                                            
##  Mcnemar's Test P-Value : 0.2482           
##                                            
##             Sensitivity : 0.9787           
##             Specificity : 1.0000           
##          Pos Pred Value : 1.0000           
##          Neg Pred Value : 0.9167           
##              Prevalence : 0.8103           
##          Detection Rate : 0.7931           
##    Detection Prevalence : 0.7931           
##       Balanced Accuracy : 0.9894           
##                                            
##        'Positive' Class : 0                
## 

matriz(cm)

4.4.4.3 Resumo

resumo[, c("VarCPT","AOV", "naov", "nmw", "desfechoaov", "desfechomw", "predito")]

save(resumo, file = "resumo.Rda")

resumo$VarCPT = ifelse(is.na(resumo$VarCPT), 0, resumo$VarCPT)
resumo$diferenca = resumo$VarCPT - (-59.61)
resumo$apl = ifelse(resumo$diferenca > 0, 1, 0)
resumo$apl = ifelse(is.na(resumo$diferenca), 1, resumo$apl) #data.frame
#resumo = resumo[resumo$apl ==1, ]
filtrar = resumo$Cluster #SpatialPolygonsDataFram
ccp <- resumo
nrow(ccp)

## [1] 174

row.names(ccp) <- row.names(APLs@data)
require(maptools)    

## Loading required package: maptools

## Checking rgeos availability: TRUE

APLs <- spCbind(APLs, ccp)
APLs$VarCPT = ifelse(is.na(APLs$VarCPT), 0, APLs$VarCPT)
length(APLs)

## [1] 174

APLs = APLs[APLs$apl ==1, ]
length(APLs)

## [1] 152

5 APLs

library(tidyverse)
dftemp = as.data.frame(cbind(1:nrow(APLs@data), APLs$diferenca))
dftemp$V2 = as.double(dftemp$V2)
dftemp = dftemp %>%
    mutate(quantile = ntile(V2, 20))
APLs$Quantis = dftemp$quantile
mapview(APLs, zcol = "diferenca", col.regions=c("orange", "red"))

mapview(APLs, zcol = "diferenca", col.regions=c("yellow", "green", "red"))

mapview(APLs, zcol = "Quantis", col.regions=c("yellow", "green", "red"))

## 
## Decreasing    Growing Stationary 
##         11          3        138