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("~/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] 16095    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] 16095

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

## [1] 16095

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

## [1] 2585

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    : 2585 
## extent      : -51.26617, -51.08242, -30.23944, -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, 601034,                 AC B VILA NOVA BRASILIA,     0, Cruzamento, ABALROAMENTO,                                    0 R CEL CLAUDINO, 01/01/2015,     DOMINGO, 00:00,          0,   1,        1, 8.29790161279961e-05,         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/2017, TERCA-FEIRA, 23:57,         23,  13,     3639,   0.0171522759325802,       230, ...

length(spDF)

## [1] 2585

class(spDF)

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

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

## [1] 2585   18

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

## [1] 2585   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] 2585   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   12   65  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.0004675219 0.0013221140 0.0139097764 0.0171522759

Locais mais densos

quantile(parqs$densidade, probs = 0.99)

##          99% 
## 117879065279

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] 2585

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

##  [1]   20   35   41   99  107  115  142  199  296  387  409  475  477  492  547
## [16]  686  826  862  896  916  925 1015 1019 1028 1083 1134 1209 1213 1267 1292
## [31] 1305 1362 1452 1458 1469 1478 1572 1586 1597 1607 1643 1686 1699 1759 1765
## [46] 1771 1825 1846 1865 1920 1963 1992 2007 2038 2050 2098 2126 2132 2228 2310
## [61] 2315 2319 2346 2430 2501 2528 2566 2567 2581

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

## [1] 2585

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

## [1] 2516

length(dnb)

## [1] 2585

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

## [1] 2516    2

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

## Neighbour list object:
## Number of regions: 2516 
## Number of nonzero links: 19158 
## Percentage nonzero weights: 0.3026418 
## Average number of links: 7.614467

length(dnb)

## [1] 2516

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,]  856  338 2112  884
## [2,]  563 1423  593 1740
## [3,] 1300  443 2200 1130
## [4,] 1001  931 1860  239
## [5,] 2160 1781 1351  295
## [6,]  148 2472 1558 1649

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.082, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1827057397     -0.0003976143      0.0001690665

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.127, p-value <
## 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1827057397     -0.0003976143      0.0001679827

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.18271, 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.493, p-value <
## 0.00000000000000022
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.8304448954      1.0000000000      0.0001841849

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.6924, p-value <
## 0.00000000000000022
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      0.8304448954      1.0000000000      0.0003060239

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.83044, 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 = 19.09, p-value < 0.00000000000000022
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##  0.004221275272300  0.003027621236018  0.000000003909745

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.09477571 -1.25380196  0.19837993  0.50535124  1.25007936  3.78966193
##  [7]  0.80693424 -0.06014283 -0.67080227  0.81556262  1.43978981 -0.89843790
## [13]  2.15651982 -1.49972334 -0.52231779 -1.07775044 -0.26400162  1.73014189
## [19] -0.48924125  0.90388678 -0.93486991  0.55476438 -0.30702859 -0.56807679
## [25]  0.61933939 -0.79582253  1.46418157  3.56010112 -0.35115130  1.08013998
## [31]  7.47831397 -1.40122076  0.23932971 -1.71522969 -1.42246295  0.55060654
## [37] -0.40678875 -0.70575363  1.44617992  2.19611017  1.52423306  0.11886609
## [43]  2.59971618  4.43797091 -0.78477079  0.72105292 -0.98717866 -0.37019483
## [49]  3.03100710 -0.33987113

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.76574385 -0.00336747335 0.493546528  1.09477571   0.2736149100
## 2   0.33676286 -0.00028834075 0.072265975  1.25380196   0.2099139952
## 3  -0.05439107 -0.00035554355 0.074193068 -0.19837993   0.8427478220
## 4  -0.15637082 -0.00022817117 0.095467694 -0.50535124   0.6133121581
## 5  -0.10609446 -0.00002295739 0.007199831 -1.25007936   0.2112705585
## 6   0.31585368 -0.00003049868 0.006947920  3.78966193   0.0001508525
## 7  -0.09887764 -0.00008986179 0.014987532 -0.80693424   0.4197043686
## 8  -0.01440817 -0.00022197656 0.055637012 -0.06014283   0.9520418784
## 9  -0.33274090 -0.00068384645 0.245039649 -0.67080227   0.5023465009
## 10 -0.19306649 -0.00020072384 0.055923674 -0.81556262   0.4147503363

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 
##  555 1028  329  604

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] 157

length(aplssub)

## [1] 157

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] 157

length(aplssub)

## [1] 157

length(APLs)

## [1] 157

#################################### resumo = resumo %>% filter(VarCPT > -59.61 | is.na(VarCPT))
#################################### APLs =subset(APLs, box_id %in% resumo$Cluster)

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] 134

APLs =subset(APLs, box_id %in% resumo$Cluster)
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] 134

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

## [1] 134

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")

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

## [1] 121

nrow(Nestacionarios)

## [1] 13

Foram encontrados entre os APLs 121 locais com média estacioária e 13 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 
##                            25                            96

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

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

4.4.2.2 Mann Whitney

table(estacionarios$nmw)

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

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

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

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 
##                            11                             2

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

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

4.4.3.2 Mann Whitney

table(Nestacionarios$nmw)

## 
##     Há Diferença entre médianas Não há Diferença entre médianas 
##                              11                               2

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

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

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 96 25
##       1  0 13
##                                          
##                Accuracy : 0.8134         
##                  95% CI : (0.737, 0.8755)
##     No Information Rate : 0.7164         
##     P-Value [Acc > NIR] : 0.00662        
##                                          
##                   Kappa : 0.427          
##                                          
##  Mcnemar's Test P-Value : 0.000001587    
##                                          
##             Sensitivity : 1.0000         
##             Specificity : 0.3421         
##          Pos Pred Value : 0.7934         
##          Neg Pred Value : 1.0000         
##              Prevalence : 0.7164         
##          Detection Rate : 0.7164         
##    Detection Prevalence : 0.9030         
##       Balanced Accuracy : 0.6711         
##                                          
##        '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 121   0
##       1   2  11
##                                           
##                Accuracy : 0.9851          
##                  95% CI : (0.9471, 0.9982)
##     No Information Rate : 0.9179          
##     P-Value [Acc > NIR] : 0.0008726       
##                                           
##                   Kappa : 0.9085          
##                                           
##  Mcnemar's Test P-Value : 0.4795001       
##                                           
##             Sensitivity : 0.9837          
##             Specificity : 1.0000          
##          Pos Pred Value : 1.0000          
##          Neg Pred Value : 0.8462          
##              Prevalence : 0.9179          
##          Detection Rate : 0.9030          
##    Detection Prevalence : 0.9030          
##       Balanced Accuracy : 0.9919          
##                                           
##        '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] 134

row.names(ccp) <- row.names(APLs@data)
require(maptools)    
APLs <- spCbind(APLs, ccp)
APLs$VarCPT = ifelse(is.na(APLs$VarCPT), 0, APLs$VarCPT)
length(APLs)

## [1] 134

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

## [1] 134

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 
##         10          3        121