Clusters de Acidentes por Affinity Propagation

Clusters de Acidentes por Affinity Propagation

library(leaflet.extras)

## Loading required package: leaflet

library(apcluster)

## 
## Attaching package: 'apcluster'

## The following object is masked from 'package:stats':
## 
##     heatmap

library(magrittr)
library(dplyr)   

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(leaflet)
library(rgdal)

## Loading required package: sp

## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-27, (SVN revision 1148)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Users/fagne/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/fagne/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was C:/Users/fagne/Documents/R/win-library/4.1/rgdal/proj

library(rgeos)

## rgeos version: 0.5-8, (SVN revision 679)
##  GEOS runtime version: 3.9.1-CAPI-1.14.2 
##  Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
##  GEOS using OverlayNG
##  Linking to sp version: 1.4-5 
##  Polygon checking: TRUE

library(geojsonio)

## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson

## 
## Attaching package: 'geojsonio'

## The following object is masked from 'package:base':
## 
##     pretty

library(mapview)
library(contoureR)

## Loading required package: geometry

library(geosphere)

## 
## Attaching package: 'geosphere'

## The following object is masked from 'package:geojsonio':
## 
##     centroid

library(cluster)
library(sp)
library(rgdal)
library(RColorBrewer)
library(sp)
library(foreign)

Carga de Dados

#dados = read.dbf("~/OneDrive/r-files/CIET/acidentes2020/_Base/acidentes_2014a2020_WGS84.dbf")
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

Preparação

x1 <- x2
#x2 <- x2[sample(nrow(x2), round(nrow(dados)*0.25, 0)), ]
x2 = as.data.frame(x2)
load("data/x2-20000.rda")
load("data/apres2-20000.rda")
names(x2) = c("LONGITUDE", "LATITUDE" )
head(x2)

#save(x2, file = "data/x2-20000-975.rda")
dim(x1)

## [1] 75102     2

dim(x2)

## [1] 18776     2

Treino

#apres <- apcluster(negDistMat(r=2), x2, q=0.975)
plot(apres, x2)

A caption

summary(apres)

##   Length    Class     Mode 
##     3610 APResult       S4

#save(apres, file = "data/apres2-20000-975.rda")

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] 3610    3

exemplars = poly
#save(exemplars, file = "data/exemplars-20000-975.rda")

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)

#save(dados, file = "data/acidentes-20000-975.rda")

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("Classificados em meio a ", summary(apres)[1], "Clusters"), title="Acidentes em Porto Alegre") %>% 
  addLayersControl(overlayGroups = c("Mapa", "Acidentes", "Legenda"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  addProviderTiles(providers$CartoDB.DarkMatter)

A caption

Enriquecimento Informacional

rm(apres)
clusters_encontrados = sort(unique(dados$cluster))
#clusters_encontrados
parq = dados
poly = data.frame()
for (i in clusters_encontrados){
  temp = parq[parq$"cluster" == i,  ]
  ch1 = convexHullAM_Indexes(temp[,2],temp[,3], includeColinear=FALSE,zeroBased = FALSE)
  #print(i)
  #print(ch1)
  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
  if(nrow(temp) >= anos) {
    for (ii in ch1) {
    polying = temp[ii,]
    polying$area = area
    polying$acidentes = acidentes
    polying$centroide_lat = centr_lat
    polying$centroide_lon = centr_lon
    poly = rbind(poly, polying)
    }  
  }
}
head(poly)

tail(poly)

mean(poly$area)

## [1] 11318.93

median(poly$area)

## [1] 4801.598

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% 
##   118.1853  1516.2621  4801.5985  9730.6296 73097.2575

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

## [1] 16329    48

class(poly)

## [1] "data.frame"

pol = poly

Combinando Informações

#save(pol, file = "data/pol.Rds")
#load("data/acidentes.rda")
#load("data/apres2.rda")
#load("data/exemplars.rda")
#load("data/x2.rda")
dados = poly[,c(1:9, 11:13,38,41,44:48)]
head(dados)

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)

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

## class       : SpatialPointsDataFrame 
## features    : 16329 
## extent      : -51.25813, -51.08553, -30.22011, -29.96328  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 19
## names       :     ID,                          log1,                                    Log2,  Pred,      Local,         Tipo,                          Via,       Data,         Dia,  Hora, Fx_horaria, UPS, box_id,             Area, Acidentes, ... 
## min values  : 601003, AC B CONJ RES ALTO PETROPOLIS, AC DOIS VOLUNTARIOS-PRES CASTELO BRANCO,     0, Cruzamento, ABALROAMENTO, 0 AV PRESIDENTE JOAO GOULART, 01/01/2015,     DOMINGO, 00:00,          0,   1,      1,  119.44694154337,         6, ... 
## max values  : 683133,            VDT LEONEL BRIZOLA,                             VDT OBIRICI, 24108, Logradouro,   TOMBAMENTO,           VDT LEONEL BRIZOLA, 31/12/2020, TERCA-FEIRA, 23:59,         23,  13,   3610, 18909.1945320629,       308, ...

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

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

Criação de Polígonos

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

#Criamos o SpatialPolygonsDataFrame
data$Log2 = NULL
spDF <- points2polygons(df,data)
spDF

## class       : SpatialPolygonsDataFrame 
## features    : 2620 
## extent      : -51.25813, -51.08553, -30.22011, -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 DIVINA PROVIDENCIA,     0, Cruzamento, ABALROAMENTO,              1 AV EDVALDO PEREIRA PAIVA, 01/01/2015,     DOMINGO, 00:01,          0,   1,        1,  119.44694154337,         6, ... 
## max values  :   3610, 683107,       VDT DOM PEDRO I, 15555, Logradouro,   TOMBAMENTO, VDT DOM PEDRO I & AV BORGES DE MEDEIROS, 31/12/2017, TERCA-FEIRA, 23:58,         23,  13,     3610, 18909.1945320629,       308, ...

class(spDF)

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

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

## [1] 2620   18

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

## [1] 2620   24

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

A caption

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

Acidentes por Cluster

#carregamos os dados SpatialPolygonsDataFrame
parqs <- geojsonio::geojson_read("data/myParq.json",
                                 what = "sp")
#Verificamos o objeto
parqs

## class       : SpatialPolygonsDataFrame 
## features    : 2620 
## extent      : -51.25813, -51.08553, -30.22011, -29.96328  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 24
## names       : box_id,     ID,                  log1,  Pred,      Local,         Tipo,                                     Via,       Data,         Dia,     Hora, Fx_horaria, UPS,           Area.x, Acidentes.x,    CentLon.x, ... 
## min values  :      1, 601034, AC DIVINA PROVIDENCIA,     0, Cruzamento, ABALROAMENTO,              1 AV EDVALDO PEREIRA PAIVA, 01/01/2015,     DOMINGO, 00:01:00,          0,   1,  119.44694154337,           6, -51.25697886, ... 
## max values  :   3610, 683107,       VDT DOM PEDRO I, 15555, Logradouro,   TOMBAMENTO, VDT DOM PEDRO I & AV BORGES DE MEDEIROS, 31/12/2017, TERCA-FEIRA, 23:58:00,         23,  13, 18909.1945320629,         308, -51.08580303, ...

dim(parqs)

## [1] 2620   24

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:magrittr':
## 
##     extract

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

parqs@data$densidade = (parqs@data$Acidentes.x/parqs@data$Area.x)*1000000
mapview(parqs, zcol = "densidade", col.regions=brewer.pal(9, "YlOrRd"))

Acidentes por poligono

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

A caption

Acidentes por KM2

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

A caption

Locais mais densos

quantile(parqs$densidade, probs = 0.99)

##      99% 
## 104685.9

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

Calculando a matriz de vizinhanças

ccods = coordinates(parqs)
points = cbind(ccods[,1],ccods[,2])
head(points)

##        [,1]      [,2]
## 1 -51.18244 -30.03713
## 2 -51.19929 -30.04108
## 3 -51.09693 -30.04702
## 4 -51.19974 -30.07228
## 5 -51.20602 -30.06893
## 6 -51.12612 -29.98625

library(spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

dnb = dnearneigh(points,0,0.1)
class(dnb)

## [1] "nb"

Matriz de Vizinhanca

Matriz Binária

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

Matriz Normalizada

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

KNN

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

## [1] "knn"

head(vizinhos_4$nn)

##      [,1] [,2] [,3] [,4]
## [1,] 1543 1798  403 2042
## [2,]  327  893 1507 2386
## [3,]  604 2085 2153 2004
## [4,]  820  377 1824 2014
## [5,]  397 1498  438 1769
## [6,] 1086 1975  822 1357

vizinhanca_4 <- knn2nb(vizinhos_4)
class(vizinhanca_4)

## [1] "nb"

Preparação para Análisis Global e Local

mv_simpl = st_as_sf(parqs)
plot(mv_simpl)

A caption

class(mv_simpl)

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

library(dplyr)
mv_simpl =  mv_simpl %>% dplyr::select(Acidentes.y)
#mv_simpl <- st_simplify(mv_simpl, preserveTopology = FALSE,                  dTolerance = 1)
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)

## although coordinates are longitude/latitude, st_intersects assumes that they are planar

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

Calculando o Índice de Moran Global

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

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 = 30.898, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.060209e-02     -3.818251e-04      1.263708e-07

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 = 31.031, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.060209e-02     -3.818251e-04      1.252931e-07

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.010602, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

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)

## The default for subtract_mean_in_numerator set TRUE from February 2016

## 
##  Monte-Carlo simulation of Empirical Bayes Index (mean subtracted)
## 
## data:  cases: parqs$Acidentes.y, risk population: parqs$Area.y
## weights: nb2listw(dnb, style = "B", zero.policy = TRUE)
## number of simulations + 1: 1000
## 
## statistic = -0.00091013, observed rank = 3, p-value = 0.997
## alternative hypothesis: greater

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)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  shapeCG.p 
## weights: nb2listw(dnb, style = "B", zero.policy = TRUE)  
## number of simulations + 1: 1000 
## 
## statistic = -0.00086881, observed rank = 5, p-value = 0.995
## alternative hypothesis: greater

Calculando a Estatística C de Geary Global

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 = -1.3298, p-value = 0.9082
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      1.003540e+00      1.000000e+00      7.087918e-06

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 = -0.37695, p-value = 0.6469
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##      1.003540e+00      1.000000e+00      8.821701e-05

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 = 1.0035, observed rank = 605, p-value = 0.605
## alternative hypothesis: greater

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 = 8.3096, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       8.308567e-01       7.792060e-01       3.863626e-05

Getis e Ord Local

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

##    [1]  6.01685436  5.67440690 -3.37928997  5.26582295  5.03563680  1.70771957
##    [7]  7.63086602  6.19643705  4.97635473  6.71253736  6.86579195  4.93681362
##   [13]  6.19109007  7.16704671  6.18958494  2.82888719  5.43473558  6.22438030
##   [19]  6.27687466  3.65954636  4.20355797  6.47510024  6.78446317  2.05472838
##   [25]  5.10115600  7.64283696  1.47170290 -0.01906698  6.14239760  5.67989616
##   [31]  0.66930228  6.27308513  5.49404050  7.39219611  5.85986329  6.53629027
##   [37]  4.53965593  0.09633320  6.40132150  6.19996535  2.02726645 -0.53569896
##   [43]  6.62742340 -0.14942660 -6.26690921  2.98056225  7.03085766  7.23541086
##   [49]  6.34314511  5.40965884  6.99557879  6.61237075  6.99233859  6.87712349
##   [55]  6.18328555  5.51849889  5.80261509  6.54346540  4.58423153  8.10255101
##   [61]  5.07347978  6.43836004  5.96980470  4.61425504  6.73410690  6.67007186
##   [67]  6.66513145  5.82824274  5.05664741  6.05519402  5.83397595  5.96162267
##   [73]  6.85916433  2.04342743 -4.05310267 -5.89504245  6.43994722  6.23091517
##   [79]  5.08059641  6.15565243  1.62011232 -1.33302280  5.66485024  2.46488103
##   [85] -1.23891613  7.04792101 -0.75080368 -1.99895763  3.63593495  6.18737267
##   [91]  2.55152746  6.80344117  6.42880780  6.90904358  5.23501531  6.94131804
##   [97]  7.59836496  6.33919023 -2.56625119  2.39361194  6.76882591  7.62829789
##  [103]  6.47187986  0.43207972  4.18959718  7.34171373  6.14820882  7.82930627
##  [109] -0.62616072  6.07907778  5.86659113  6.75496079  5.37702981  2.41173411
##  [115]  7.63039328  2.65394726  7.27869374  5.82632262 -0.45729687 -0.72059473
##  [121]  6.87649960  4.04123623  2.59503446  7.66678752 -0.92707885  5.48085038
##  [127]  4.99648795  6.58034641 -5.30344766  5.99452315  6.51025063  5.67838926
##  [133]  6.03215108  7.01676053  5.93380502  7.46569587  6.07160847 -6.23917962
##  [139]  6.11079034 -4.90147309  5.90629180  5.83317673 -4.71445825  2.69183262
##  [145]  6.05784103  5.88902591  5.87268316  6.01604011  6.59030064 -5.90754997
##  [151]  7.66169775  7.44586684  4.11822981  6.92013801  5.80673195  2.81476616
##  [157]  6.23152630  5.04635376  6.43710556  5.59296923 -1.64853467  6.47363612
##  [163]  6.36423168  5.98558813  6.12923088  6.40777025  6.60298264  4.26004683
##  [169]  6.28379761  1.55370276 -1.02975744  6.05936059  6.35980699  6.02454208
##  [175]  6.71928478  6.62983652  0.93991525 -0.46809353  5.78723335  4.61279456
##  [181]  5.78518634  7.52661465  4.45842859  6.18187577  5.66729097  3.58252827
##  [187]  6.77688750  7.28519236 -7.25521999  5.17627047  5.30729171  5.68861695
##  [193]  6.71072625  1.00866656  6.03112913  5.47451837 -0.77437531  5.52454589
##  [199] -0.33305507 -6.53735010  4.84008297  4.18704083  5.90350025  7.99249975
##  [205]  5.24666556  7.29836935  6.82234298  7.28558863  7.41792541  6.46008243
##  [211]  6.45398709  6.38798449 -3.64447659  6.33586105  5.76823825  6.29947202
##  [217]  7.87646084  5.14775000  4.61006448  5.74061573  5.74950330  6.54803772
##  [223]  1.90210644  6.57514547  6.55943546  7.60842181 -4.61919055  7.56389793
##  [229]  7.15225139  6.13393294  5.50622431  5.82314138  4.30405876  3.94602518
##  [235]  6.16822005 -5.10009224  6.46353197  5.07950199  7.98527775  4.25392272
##  [241]  6.82130294  4.28756263  5.43144622  4.29294649  3.41126757  6.66111100
##  [247] -0.14406064  7.68531677  6.35157487  6.16443931 -0.61417553  5.94166228
##  [253]  6.36223669  5.58854698 -5.22143415  6.25188912  6.35594673  3.98899846
##  [259]  4.79504436  6.59630601 -5.38585982  5.77074339  6.52924370  7.54044778
##  [265]  6.17303093  6.02900720  5.99962310  5.70263516  5.30238647  8.01607902
##  [271]  6.09316431  5.80762188  5.42209420  6.34035307  2.68344460  5.56825900
##  [277]  2.38564588  6.72789190  6.45933624  3.92046952  7.11390292  6.19178308
##  [283]  8.11787612  6.19943812  2.20282868  6.39275058  7.60120386  7.07549664
##  [289]  6.45520224  5.63587730  6.30996626  6.08854949  6.25552349  5.34052143
##  [295]  6.29000725  5.55016927  6.27830330  6.34224391  6.77745531  5.37521399
##  [301]  6.45039068 -0.51234236  5.92583847 -6.46884244  0.82485475  6.75730253
##  [307]  7.20218562  6.16340660  4.75691620  7.21782453  5.59318452 -1.60841101
##  [313]  4.24815171  6.11845015  1.04781651  7.09398730  6.92804166  6.32031093
##  [319]  2.47351379  2.06358855  4.43598607 -0.12363659  6.14360245  4.67919259
##  [325] -3.83643297  5.07147695  5.55197445  7.06470881  1.69566108  2.72169919
##  [331]  7.91013214  4.70237186  5.71615600  7.52703194  6.52390186  5.26607040
##  [337] -4.62839257  2.57063655  6.21919798  6.29097696  6.01371974  3.78464018
##  [343]  6.73921653  4.86412466  6.21846336  6.33367792  6.37756561  6.94509208
##  [349]  5.41858843  2.38905626  5.58415276 -6.45810295  5.46195933  6.74577784
##  [355]  3.14434618 -4.42922332  7.17631357  5.49774364  5.30743085  5.68970671
##  [361]  0.46515766  6.28140220  6.45196982 -0.69931279  4.68554612  7.06296496
##  [367]  4.94122493  2.50557359 -4.26874842  5.03924367  7.99479176  7.69792929
##  [373]  6.67751580  0.66206354  1.97521202  6.38097620  5.14377400  6.38628436
##  [379]  5.89446027  7.91286764  5.42950030  6.76446306  7.79930902  6.62840274
##  [385]  6.79241727  6.44153662  6.22115167  6.19926502  2.49070919  5.63551644
##  [391] -6.20204290  6.00174150  7.54714487 -5.81556733  6.24417494  6.04566539
##  [397]  5.13309131  6.37177524  5.82052049  6.42690321  6.69024112  5.94777143
##  [403]  5.95744791 -6.03385959  2.55731681  5.43148430  6.42823099 -5.55702926
##  [409]  2.44752091 -3.66559640  5.02682832  5.55616668  6.73198020  6.52950075
##  [415] -4.71757821  5.15975388  4.48811885  6.85674141  6.86179551  6.04249043
##  [421]  2.58257669  3.43017134 -0.69349950  6.82607899 -6.13406565  5.19859283
##  [427]  6.34004032  7.53208586 -1.16528120  5.76347257  5.68602381  6.73310845
##  [433] -0.42370977  5.33836324 -5.29880361  6.63020969  7.32948441  4.92132202
##  [439]  2.35870011  5.78809860  2.48054882  5.71657487  2.85683035  7.67085582
##  [445]  6.58784742  5.39352840  4.29470100  6.41978495  6.01982473  6.76929016
##  [451]  6.91314886  6.14155161  6.08394189  7.06164509  5.74009012  6.89893213
##  [457]  5.61380292  4.42781018  7.68422594  5.21866037  6.37921206  7.04949217
##  [463] -5.92174773  4.48425271  2.77346649  7.98043237  6.12360866  5.99914907
##  [469]  6.85084958  6.81691140  5.54727885  6.42317182  5.47604278  6.36778397
##  [475]  6.67820916  5.47738807  6.03752576  6.86516712  2.94810270  6.66558245
##  [481]  2.82215517  6.39544831  6.17657688  8.06376804  5.75265378  5.69814446
##  [487]  5.54828331  7.83717409 -5.61475229  7.85997024  6.03405640  7.47916719
##  [493]  6.72040115  6.25087315  7.03958056  6.96376646 -0.42744834  5.01600141
##  [499]  7.56484541  5.96346946  6.45192227  6.63074489  5.55475566  6.34870317
##  [505]  6.54370281  2.53352052  4.86600067  5.21840721  7.05433457  4.12356486
##  [511]  6.85816734  7.73949308 -0.81612692  6.67068395  5.11285186  6.73075106
##  [517] -6.29892843  6.84098251  5.91198935  5.55558373  5.12525970  6.54710418
##  [523]  2.01563818  7.87873218  6.17490077 -6.30750997 -5.88894000  6.71085026
##  [529]  5.32109216  4.63633853  5.60245633  1.90185426  6.32736981  2.93202121
##  [535]  6.98905534  2.78226020  6.28065324  5.86809629  6.60187708 -4.81317195
##  [541]  2.61872212 -1.51926916  6.49393322  7.39442608  1.53496874  8.00825413
##  [547] -0.38760349  3.01495236  0.04041500  5.12370011  6.06787807  5.71954572
##  [553]  7.58680688  5.75347385  3.45807253  6.04975334  5.95832968  5.93577484
##  [559]  6.78116996  3.04547480  3.86543911  6.27024066  7.20563330  5.69870894
##  [565]  5.57612237  5.13624411  6.16242713 -6.99875610  1.92111589  6.00236125
##  [571]  6.10367776  5.42394475  6.80176719 -1.37487450  6.11891835  4.50460707
##  [577]  6.78916585  3.81682408  6.22490514  6.95372087 -5.08543490 -0.24543175
##  [583]  7.41856193 -6.08017513  6.66148343  5.91949904  7.16598211 -0.45295624
##  [589]  0.18965202  5.42710893  6.96945811  6.76187868  4.21059797  5.27667519
##  [595]  6.31044695  7.62182013  5.28030850  5.93256338  5.60231346  5.51058745
##  [601] -6.59954274  5.46056232  5.81605427 -3.59420346  6.48870140  1.62149020
##  [607]  5.71095314  6.51632119  6.01148808  4.91232369  2.52466432  7.01540640
##  [613]  5.45089043 -6.57075891  4.95580831  5.04369620  6.48184919  4.29552991
##  [619]  6.63785829  6.53040879  5.21409448  3.41824164  6.69538414 -0.51307503
##  [625]  4.96484987  6.51841166  2.48326884  6.11805163  5.59763644  2.75887750
##  [631]  7.23700540  2.91364698  3.21930691  5.85174130  6.49785822  4.10621108
##  [637]  6.38949465  5.43552799  5.44810973  5.96137462  7.35597728  6.25650642
##  [643] -6.25579411  6.89787131  0.09109833  7.78893212  7.05337049  5.12424942
##  [649]  0.86779686  5.54306325  5.94918094 -1.22143182  2.85635539  6.18067925
##  [655]  3.96658020  7.88532396  5.64522584  7.60248696 -0.59140144  0.31349758
##  [661]  6.85827059  6.99638487  6.31282082 -1.04336438  4.00480426  6.53013517
##  [667]  6.35040055 -6.20796377  6.56329241  5.66892917  6.03131900  6.10057573
##  [673]  5.72522935  6.85916494  6.62978298  6.61276026  5.21239724  6.09417231
##  [679] -5.30749237  6.69510718  5.42530686  6.35995441  5.79826333  5.31972398
##  [685]  3.55252004  7.07092464 -6.77246597  6.67956556  1.67918630  5.23398488
##  [691]  6.82016990 -6.89529223  4.21813124  4.91463930  2.85320716  6.43636860
##  [697]  5.29791101  4.16712773  7.22895886 -0.05625871  6.51753073  1.11456371
##  [703]  5.25982133  2.83407940  6.29639562  3.18260680  5.74733450 -6.05797061
##  [709]  7.63937780  6.41478331  6.49685263 -6.74819187 -5.58163866  6.36735684
##  [715]  2.92498765  5.70028114  3.21962682  6.88468892 -4.47071234  4.67391612
##  [721]  6.14314440 -1.26457263  3.89472203  6.43342978  7.45783226  6.22605078
##  [727]  1.47231301  5.60652911  3.62598022 -6.59715883 -2.80904921  6.24499717
##  [733]  6.14004177  6.43943574  4.34897794  6.92671840  6.06703802  2.68351597
##  [739]  6.07093102  6.87687022  7.70901526 -6.42383203  6.49567062  6.10192718
##  [745]  7.53534147  5.81617682  6.57316085  5.47077596  3.53139553  6.22211583
##  [751] -0.44477046  5.42457640  5.19999196  4.77104741  5.43910931  5.55878340
##  [757]  7.07355481  4.36546590  7.18880467  1.89606245  7.70227330  6.96478696
##  [763]  5.92354474  6.12966830  6.41163789  5.28569484  5.79108147  5.97726910
##  [769]  2.08225139  6.15388705 -6.49940605  5.59691351  5.84251628  6.71882327
##  [775]  5.21048977  7.64644683  7.34627684  4.97703040  1.06566198  5.91824705
##  [781]  6.38130653  7.88221566  7.99628288 -0.07806224 -0.58043889 -5.15099010
##  [787]  6.58334378  6.59227011  7.44509703 -0.80111062  6.32267561  6.84660939
##  [793]  5.27187670  6.07883689  5.88260605  5.86408540  4.90659871  5.63174893
##  [799] -1.81984681  6.04485046  8.05825566  6.98263107  6.65611138  6.49918482
##  [805]  6.70537774  6.78992569  2.55301813 -4.09623619  5.68306858  5.18581516
##  [811]  6.43145348  7.49014602 -1.91594863  7.45333303 -5.31792088  6.19445756
##  [817]  6.15856687  7.46543501  7.16637591  4.96097718 -1.43008382  1.43555046
##  [823]  4.29912556  5.98892450  2.45983987  6.28718783  6.87636436  6.16020832
##  [829]  4.93670566  0.69794422  7.00383948  4.09096289  6.52743914  7.08925856
##  [835] -7.16897363  5.88758951  7.41871500 -1.08916048  5.81994567  4.06025759
##  [841]  1.10173674  5.47963921  6.80904985  7.16795372  7.75782311  6.65061881
##  [847]  6.49332515  4.45270840  5.32028637  5.18701596  7.47645527  6.18618964
##  [853]  6.51723789  6.82579074  6.94089680  6.78861892  5.69229130  5.09519028
##  [859] -0.94244706  6.33007536  5.66681580  6.55955360  6.33874817  7.56856680
##  [865]  6.93804346  6.92603987  6.27772494 -0.93642465 -5.00243443  7.44636936
##  [871]  0.25695102  4.76090470  6.70674887  6.14086079  6.74822370 -4.74059482
##  [877]  6.48489882  1.75310806 -0.17534617  7.03288289  7.79138587  6.11542081
##  [883]  7.71463487  7.57369205  5.30859677  3.66148044  2.21201097  5.55075740
##  [889] -6.71338616  5.23121789  7.12938458  6.89231408  5.49167207  7.49240016
##  [895]  7.09917477 -0.76162579  4.43807104 -6.05206525 -1.08394224  5.34384714
##  [901]  6.49643956 -6.15439689  2.42056880  7.24131804 -0.15956640 -6.57364784
##  [907]  6.43014356  6.47432994  6.25555950  6.91759490  6.50938897  3.72930930
##  [913]  5.92724186  6.67671357  2.27728586  5.08329606  6.48219986  5.98789188
##  [919]  6.46170143  5.98740674  5.80352099  6.45994779  6.39064751  4.19581515
##  [925] -5.41437691  5.60601149  6.35067994  2.32219928  7.15950201  5.71279485
##  [931]  7.59358226  6.98596094  6.25566563  4.23222449  6.15153589 -0.97834525
##  [937] -3.45232933  7.63047428  5.54598490  6.31564231  8.27233542  7.13739825
##  [943]  5.18343280  5.89194676 -0.68188204  6.78913831  6.37235365  7.18056525
##  [949]  2.16621743  5.39819558  5.74689684  6.66748322  5.38986610  2.83023056
##  [955]  6.79003567 -5.20965145  5.95819830  6.55212233  1.90123022  7.60379280
##  [961]  7.68952343  6.93090697  7.14750678  6.11954323  4.74818962  2.63137294
##  [967]  4.42536191 -5.48315401  5.67838187  6.87872172  2.27174939  6.31910114
##  [973]  6.05829040  7.01110968 -5.79160684  4.39775713  7.82612163  4.91013222
##  [979]  2.41602355  5.36773533  1.67284157 -0.86716511  7.55633907  4.14218656
##  [985]  5.24350573  2.33777544  6.09181989  5.86874297  3.54429492  6.42860567
##  [991]  1.40971413 -1.21005755  6.54519149 -0.57036536  5.55219987  5.98141416
##  [997] -0.34755611  2.26862725 -5.38231834  7.93610089  4.73021123  6.82528361
## [1003]  6.83424871  5.68924138  6.87158079  6.16476843  6.88969115 -0.83709251
## [1009]  5.79082080  7.31857003  6.23104451  6.54747750  6.28375260  5.71113249
## [1015]  2.11077818  6.32398323  6.98557279  5.07248443  6.37966635  5.70198560
## [1021]  5.65340725  6.54212402 -4.69634951  6.24449283 -0.58753404 -1.53976919
## [1027]  1.87048014  7.37559892  5.31835471  3.22997998  6.31526178  7.34194411
## [1033]  4.51796523  6.77435135  2.81608430  7.00472317  2.44485125 -7.49511479
## [1039]  5.91423082 -1.16781149 -4.95101344  6.16594275  5.77836430  4.35993874
## [1045]  6.26199042  5.96444193 -0.49951980  6.13624786  4.91863036  3.86580842
## [1051]  5.21009057  5.92845777  5.38539664  7.63476900  7.19334045 -6.30055826
## [1057]  5.53928551 -5.90917153  5.77269127  6.16146772  6.97790361  7.32620112
## [1063]  5.87989636  5.98342861  2.79331040  6.69524155  5.83306560  6.84589254
## [1069]  6.54579981  7.45713793  7.17624567  2.17877998  6.03480008  5.47881299
## [1075] -4.51906882  7.05274461  5.60314586  5.76441699  2.53943116  4.29082893
## [1081]  3.47152537  6.19049708  5.35617427  5.74759082  6.23695884  1.53143410
## [1087]  6.91742796 -4.33674987  5.45554679  5.12376775  5.66932269  6.86308744
## [1093]  5.34029893 -1.70508009 -5.85286803  5.72394958  7.51998908  6.64522391
## [1099] -6.02209792  6.70200116 -4.67640840  6.25298729  5.54846783  6.32674270
## [1105]  6.30955387  6.39465706  6.73854644  7.02529680  6.63058196  7.47255707
## [1111]  6.83334636  5.10912788  8.10309567  5.65351501  5.09903085  4.71116009
## [1117]  5.49447488  6.11025171  5.23693052 -0.57860706  7.15069488  7.04556449
## [1123]  4.93875623  5.84167253  6.32095162 -0.68409184  5.18182432  4.41297438
## [1129]  6.86672979  6.51298999  0.40770985  5.84786684  1.82222659  5.64419521
## [1135]  6.77547258  7.12087930 -0.75105153  6.40062549  6.70424037  4.62868701
## [1141]  6.32368752  6.16797714 -7.18416935  5.57120575  4.87314923 -2.06805031
## [1147]  7.54423012  3.77589657  5.69289764  5.88780650  6.06333124  6.32578546
## [1153]  3.10961116  6.42216479  4.07208880 -0.82949888 -6.55283522 -0.42827484
## [1159]  5.46243191  2.76905183  5.47977963  6.35564559  7.22193290 -0.91854681
## [1165]  4.78341213  0.72131540  4.32275148  4.85012602  0.33533928  7.66079814
## [1171]  2.38930675  5.68776131  2.34186877  6.22456049  7.09228038  5.87324287
## [1177]  5.99566007  5.70965986  6.81540872  6.00268840 -5.44674572  6.41964494
## [1183]  6.66033102  7.20749491  6.25271356  7.24051615  7.79088813  7.82563911
## [1189]  4.63543114  6.91312280  5.70983487  6.84551038  6.45401092  2.05020647
## [1195]  6.81993462  5.89187085  7.06691291  6.00322609  5.91888004  4.12284377
## [1201]  6.33573084  7.15734637  6.77196493  4.11406909 -6.85855713  5.73322201
## [1207] -4.67910940  0.93224468  7.71229273  5.65406517  4.58678990  2.11395179
## [1213]  6.95225290  6.16051500 -6.46971844  6.19115760  7.64229405  6.89015849
## [1219]  6.73670534  2.60359775  7.52966262  6.89354158  7.79248656  5.60691235
## [1225] -7.22239936  7.06131319  6.36740658  7.13503119  5.91349349  6.10652719
## [1231]  4.44149813  6.93595907  6.08497768  6.32158785  3.81028322 -3.05738209
## [1237] -0.93024476  6.89018889  6.08681672  6.53322755  4.67938320  4.76745209
## [1243]  5.67239371 -1.81073740  5.46558755  6.83386457  2.66881332  6.50105216
## [1249]  5.86614382  4.23500013  6.39830336  4.63558512  5.89473603 -0.79426035
## [1255]  4.90903356 -0.86344814  8.24648465  5.15362235 -4.25364289  6.19299496
## [1261]  7.50594895 -3.54620926  7.67498737  4.36908908  6.92785056 -5.81247876
## [1267]  6.19063758  6.50487279  6.02457869  6.05849007 -0.45656081  5.52144547
## [1273]  5.71919099  5.92333043  6.28459093 -4.65907158  5.43093211  5.40616396
## [1279]  6.66393437  7.58761381  5.95518680  7.72002497  2.80922642  4.88482249
## [1285]  6.24120188  6.04729596  5.36368938  1.98826572  6.81202751  7.46505291
## [1291] -2.84873299  6.59189186  6.41805162  1.72589244  6.73181329 -0.36711104
## [1297]  2.19712657  6.70058368  6.44085801  6.52652969  1.01189803  3.85094095
## [1303]  6.91210759  5.85907864  5.99152930  6.90777812  2.67679880  5.78386913
## [1309]  8.21753541  4.22641005  6.79044649  6.19535145  3.55582018 -0.13681662
## [1315]  5.35512706  3.73798138  7.07661377  6.00503601 -0.41715721  6.48075636
## [1321]  6.18115696  6.12129920  1.51936362  4.75933261  6.96330823  6.09781869
## [1327]  2.64825617  7.80058139  7.65248045  4.85817682  5.29604479  6.15362389
## [1333]  3.65286323  7.69278535  6.93826943  6.33040469  4.14263640  6.19878743
## [1339]  3.35536113  3.01559015  6.33272137  1.94924701 -5.04216989  5.57942898
## [1345]  5.85955794  2.36928861  0.77097928 -1.22761739  2.48576961  6.20322749
## [1351]  7.73573119  4.69096010  5.85126938 -0.41794977 -0.29568621  5.56941007
## [1357]  1.65969596  6.82222384 -2.98643937  2.43201987  7.76425464 -1.14808532
## [1363]  5.89298122  4.38548903  7.29604254  2.63896131  5.14111012  6.12701320
## [1369] -6.21542300  6.23325638  6.80993094  7.63518879  6.77059067  2.16653335
## [1375] -0.36412787  6.65905580  7.85969210  7.26724922  5.95103486  6.06204281
## [1381]  3.88632603  5.91108965  5.88269874  6.32518497  6.96620305  6.69518416
## [1387]  2.89182663  6.39704393  6.74089752  6.03866805  2.34435943  6.38797382
## [1393] -3.72597276  7.93410156  5.83964048  6.34852855  4.26230237  7.50800687
## [1399]  4.24535828  6.78459223 -0.65717838  7.50030064  5.92551258  7.36554852
## [1405]  6.14270318  6.83156557 -5.21605942 -1.33155366  4.17314241  6.09268872
## [1411]  5.46401485 -7.32646243  6.10413085  4.23900704  2.64940250  5.49277673
## [1417]  5.79692585  6.82315245  5.15216013  6.38263741 -1.48464740  7.87144269
## [1423]  7.13848959  6.13668654  8.09164070  4.91267621  2.32736650 -6.48988183
## [1429]  5.74280210  5.56013825  1.36163181 -5.64422055  5.79587441  5.68718468
## [1435]  2.32182611  6.42142216  5.30579995  4.86211416  5.70255408  2.19977740
## [1441]  6.50107317  5.93030211  5.46292899  4.08053601  6.24717778  6.12870307
## [1447]  5.26517213  7.11775053  5.97287959 -6.54588701  6.93980870  5.78593574
## [1453]  6.90428677  5.24589755  6.37163916  6.86932452  0.39941586  7.50118782
## [1459]  6.86970578  5.72070611  2.28472474 -7.29079198  7.28752416  6.18127070
## [1465]  5.63697193 -4.77116239  6.11356648  6.60245112  5.93284083  1.55945646
## [1471]  6.35977271 -0.29509405  6.39720544 -2.49121455 -1.81530265  6.65895012
## [1477]  2.11696321 -6.15351747  5.09362832  4.88574379  5.81340088  2.11743763
## [1483]  6.65576310  2.07471618  5.63063413  6.25989964  2.40333312  1.56496587
## [1489]  6.53500670  6.31036169 -1.03527476  5.68908620  5.62479208  6.19300632
## [1495]  6.99249231  7.05546191  6.13186757  5.18573283  5.24642717  4.99086661
## [1501]  4.21853725  4.84085732  5.80382147 -4.15586516  5.64883029  4.71643482
## [1507]  5.87907121  3.16695253 -0.88343945  6.01035112  6.29707544  7.00892154
## [1513]  5.62943864  6.06918906  7.71580840  5.53770638  6.80508517 -1.02843492
## [1519]  6.83225952  6.60036080  7.28266943  6.81504259  7.12405613  7.85661658
## [1525]  5.46559106  7.95843765  6.66972715  7.39310940  7.00224058  2.55465695
## [1531]  6.35770365  6.05282364  5.34601709 -3.01498712  6.18094847  5.71028829
## [1537]  7.56638872  6.87916704  5.56257063 -5.31885846  6.31755042  5.77573472
## [1543]  6.02417171  6.37180405  7.03024438  1.44487504  5.19776654  8.02463830
## [1549]  7.52424801  7.55545625  2.11176475  5.57601180  5.45365613  5.96628350
## [1555]  4.23967934  6.32650234  4.69266693  2.17816832  5.33863524  5.42810892
## [1561]  5.36473803  6.04154607  5.94966364  7.31562237  6.96936661  5.91728550
## [1567]  2.75286015  6.15971252 -6.27615839  7.38642812  6.82118071 -7.19837958
## [1573]  2.72509606 -0.14688024  6.74977864 -0.18654060 -7.51817101  6.93673564
## [1579]  6.24732317  5.05722405  6.53213806  6.05496163  5.92935703  5.96338709
## [1585]  6.34644896  7.46538404  6.51460961  5.50715740  5.92916702  4.80527048
## [1591]  6.19623236  6.14564599  7.10806060  5.84602763  6.90611486 -6.36708140
## [1597]  5.35391323  6.75697601  2.88004903  5.85740214  4.85552051 -6.04522020
## [1603]  6.14771724  5.64998090  1.14104851  5.49261910  5.94302009  6.37566949
## [1609]  6.31672062  5.14043010  5.92692793  6.94100682  6.46830799  6.93914050
## [1615] -0.12595256  6.19388505 -0.25916713  6.81987221  6.76184111  6.10131086
## [1621]  5.66930557  6.82828910  7.73876895  5.44638993  6.91680825 -3.87937797
## [1627]  5.45494875  5.72914173  3.01230497  3.87762223  6.63800372  5.17323746
## [1633]  6.02491357 -5.07999227  6.13361637  6.35938603 -0.56283377  6.29639805
## [1639] -1.06097103  2.30277663  6.81775124 -0.74365292  0.63290588  2.96912303
## [1645]  3.12450855  0.80976140  6.03988943 -7.06064823  3.08597047  6.34008825
## [1651]  6.37723063  6.36330812 -6.25710938  6.17145618  5.71714298  5.92884030
## [1657]  2.28427969  6.32049712  5.99885618  6.06454705  6.00174807  2.07425525
## [1663]  3.30730745 -6.86544765  5.39606893  5.71973336  5.78273165  5.03854613
## [1669]  6.14848811  6.95531969  6.12842040  5.18531864  6.08483684  5.72237067
## [1675]  5.14451575 -3.38628464  6.16559272  6.17420490  6.35012978  0.17121582
## [1681] -0.50597218  6.82925859  6.97000959  3.22772307  5.44955354  6.01482623
## [1687]  6.50007910  5.37675209  7.24335749 -5.17709744 -1.22445447  7.07299050
## [1693]  5.45292622  5.59240692  7.66730421  2.37058472  7.04468015  6.61264688
## [1699]  3.63172021 -2.21282297  6.76279349  6.79842496  5.93244001  6.32307558
## [1705]  5.55573786  0.38540959  5.83357279  6.43236572  6.55140557  4.35868977
## [1711]  5.41278789 -0.36058050  1.20567506  6.56980923  6.13158471  1.32370919
## [1717]  6.38701455 -5.86083610  6.05285233  5.75971847  4.44663232  7.00986602
## [1723]  6.16748655  6.11214590  5.25327801 -0.29372799  6.35515711  5.47862913
## [1729] -4.26095707  6.55083251  7.10522668  5.33068574 -4.23085507  6.90260015
## [1735]  2.70576833  4.85338637  2.79385147  2.07779525  6.19242437  5.58693338
## [1741]  5.24207812  5.35042512  6.21760641  6.93574085  5.84205429 -5.43613734
## [1747] -6.33816190 -0.47479089  6.12201295  6.77421616  4.69755310  6.42676293
## [1753]  6.47243964  7.71871173  4.25530495  6.86304502  5.95871350 -4.66032748
## [1759]  6.20505387  0.98270293  5.66541476  2.01767442  5.74505374  3.44194480
## [1765]  5.81080343  4.59588221  5.55590149  6.38695353  4.89339182  2.03328237
## [1771]  5.94662201  1.81871358  6.49821041  6.34053645  6.05231669  7.59385975
## [1777] -1.50843218  6.94751134  4.62821039  0.94135050 -4.65716633  6.81582894
## [1783]  6.39272189  1.86438558  5.63814668  6.23603726  6.28594431 -0.20466061
## [1789]  6.66882858  6.32761312  6.19923042  6.79907430  5.24263470  5.81265156
## [1795]  2.57162190 -4.42271585  5.73065580  6.12960128  6.03151900  5.76535381
## [1801] -0.70490296  5.65717036 -5.66955848  6.13599259  7.02927766  5.00007031
## [1807]  5.72792074  6.66644491  5.66705591  6.32412978  6.55592635 -1.64767849
## [1813]  5.54832843  3.90792867  5.25418046  5.64864824  6.28113753 -0.50259909
## [1819] -1.03639538  7.89860145  6.83196581  7.35370526  6.14952007  4.47240617
## [1825]  5.79982260  5.38261953  6.16293226  6.44632210  5.83680620  6.47569596
## [1831]  5.17257584  3.96351599  5.48572646  2.76694649  5.23495910  6.86530189
## [1837]  6.70376683 -1.81902753  4.53514589  6.32605406  5.28151668  3.08817615
## [1843]  5.52567444  5.83602311  6.76502712  3.49177379  5.39844583  6.43596105
## [1849]  3.98012098  5.57593093  4.23678734  5.96659887  7.41425200  7.06701752
## [1855]  6.86549692  5.12653890  6.39584807  5.63095272 -0.95834798  4.92074879
## [1861]  7.58116878  6.19063961  6.66901989  7.48787059  5.86378924  6.48336443
## [1867]  6.25085778  7.95210979  6.70043505 -2.99932556  6.52297681  4.46898271
## [1873]  5.53689409  4.86752064 -4.26344584  4.16453578  6.15173716  6.78017667
## [1879]  3.02450855  6.45348244 -5.73287419 -1.45148318 -6.57635845  6.64988202
## [1885]  6.07945333  6.49071573  6.30449723  5.71461396  6.22203819  7.19039615
## [1891]  6.31865311  7.18610867  6.81187534  1.65569854  6.09580503  5.74267792
## [1897]  3.91841101  5.84074766  2.94257456  6.12485178  3.72985924 -5.36432747
## [1903]  4.57515142 -6.40193380  5.23927207  2.46282536 -1.38969680  5.63727389
## [1909]  2.45812731  6.55725598 -7.26462263  6.98485489  7.62445652  5.63591735
## [1915]  6.79623874  5.67301564  5.12555838  4.85693409  7.30715824  5.40222096
## [1921]  6.06771394  6.30424866  2.44181327 -6.85107048 -6.40454499  7.03856722
## [1927]  5.14149277 -2.35148374  6.61612796  2.14490962  3.39990805  7.08744479
## [1933]  5.93004517  6.20968731 -5.24303279  5.47490505  6.12372285  6.04786473
## [1939]  5.08888993 -1.47323063  6.80182866  6.07275533  6.68433456 -4.46823848
## [1945]  7.76112251  2.57076976  6.53539016  6.96349925  5.92440455  4.84153054
## [1951]  5.13668840  6.23596635  3.88519389  6.97012354  4.01551238  2.28731292
## [1957]  6.21978711  5.29010603  5.98484012  6.81057167  6.17102561  5.84142404
## [1963]  5.41797936  2.45241577  4.15316731  5.63272472  4.04173692  6.64009955
## [1969]  6.68304845  7.52436430  0.62027308 -2.24527697  6.17553648 -5.10787971
## [1975]  1.74196196 -6.67264039  7.53351928  6.30903343  4.22843946  4.49886848
## [1981]  2.36459955  6.02601197  7.05986772 -6.67670514 -0.14986494  6.14317517
## [1987]  7.53947017 -0.47203907  7.39326182  6.14623506  6.65785406  5.68633838
## [1993]  6.14195187  6.25846608  7.07680445  2.18224950  5.99134732  6.01973414
## [1999] -6.62947182  6.69731694  5.78938398  7.30942648  4.28318918 -2.19402859
## [2005]  6.16601490  5.34541763  6.53483537  6.28629655  5.36581469 -4.59400722
## [2011]  5.42953089  6.86702355  5.76006731  5.32330120  6.94644580  7.53364284
## [2017]  4.34065372 -3.92611125  6.56502508  6.90690198  6.97050004  6.32373311
## [2023]  6.58028018 -6.41291663  6.11876994  3.50607761  0.18663229  5.25291563
## [2029]  6.36794222  5.87320753  3.02235477 -1.50783385  2.02109508  6.56358883
## [2035]  6.23941825  4.94102469  6.43582076  5.29418937 -0.35742164 -6.50565599
## [2041]  5.78711922  6.14182517  7.29783904  0.66850782  5.42623571  7.89213184
## [2047]  7.31040752  5.42829260  6.44727824  3.57258802  6.14756100  5.34936768
## [2053]  4.68027790  7.28077002  4.71936644  5.00105803 -0.03734291  5.49276727
## [2059]  6.25637046  4.36224490  6.14964621  3.56655477  5.89503186  7.49595142
## [2065]  3.35458634  6.18491959  6.83137742  6.74045417  6.10807664  4.84035029
## [2071]  6.58134025  5.41890553  2.67021610  6.96143104  7.11452063  6.20945552
## [2077]  6.60787348  2.37400726 -0.47201296  0.15985355  6.51942948  4.80539892
## [2083] -4.54180516  5.39246365 -2.23679557  5.55176093 -1.48871920  6.48987906
## [2089]  5.81665042 -6.87209532  1.59261835  4.32643413 -1.44003008  7.94361701
## [2095]  4.90559825  6.96585436  5.78081880  6.02622202  6.30070058  5.62754800
## [2101]  5.29226724  5.67355472  6.17647499  7.12309846 -1.11608684  6.61234062
## [2107]  6.65900650 -4.31645277  7.63563547 -0.62993749  6.21280463  3.08219207
## [2113]  5.51743136  5.60814987  6.35675487  6.34788289  6.58008149  5.49843228
## [2119]  2.68547067  6.69736989  7.33194112  6.64466995  3.18068628  5.64964863
## [2125]  6.41032488  7.61941150  4.99959550 -2.34612657  0.31668462  2.18899647
## [2131]  6.80455078 -2.05400145  6.59173721  5.76730411  6.69938415  6.03718060
## [2137]  7.82116011  1.05288014  6.42201398  4.81912171 -5.25693555  3.36398098
## [2143]  6.48823549  5.76100658  5.57705658 -0.60047065 -5.42197606  2.34496320
## [2149]  6.10023481  2.62282448  3.11884046  7.03825338 -3.13022635 -1.56222168
## [2155]  6.26331892  6.00464508  6.12676635  7.00251788  6.63573996  4.79288306
## [2161]  5.85453250  2.71184272  6.64328541  4.84388752  5.99922913 -6.16221270
## [2167]  2.52612625  2.72739819  5.79285845  7.23519210  7.18439335  5.57342608
## [2173]  5.01174491  6.87681206  6.94574729  6.28738022  4.76499242  6.86357934
## [2179] -2.20439484  5.16954956  5.72574622  2.43072501  6.63348087  5.67161777
## [2185]  2.24721354  4.96606240  5.85474618  7.24610541  2.85015589  5.76469650
## [2191]  6.43069459  7.81019126  1.96164104  6.25556208 -5.26916604  6.71144278
## [2197]  6.45562390  6.45014096  5.79027945 -6.31843735  5.85321814  6.71928442
## [2203]  5.63521305  6.41527064  6.75046090  6.51704567  4.77833224  6.08523949
## [2209]  5.99047510 -5.52648477  5.51630189  4.54749980  7.54654745  6.09510177
## [2215]  4.30342323  6.84366488  6.43625409  6.09493821  6.62409632 -4.86868257
## [2221]  1.73008302  5.66246355  5.36303957  4.62531478  6.77491450  6.38836963
## [2227] -1.05055990  5.84343921  6.84548314  6.40641119  5.76211791  5.30201096
## [2233]  4.99115536  2.08761923  6.62514144  5.34576145  7.45467894  4.30085049
## [2239]  5.53715852  6.66966573  2.40315220  5.36243146  6.79426889  5.08449828
## [2245] -1.24214140  6.63119316  6.12878390  6.19839574 -4.15655195  7.05586159
## [2251]  5.41440832  7.41580892  6.94741450 -0.51025118  6.81254930  4.68205974
## [2257]  1.36245623  6.34482836  6.30937577  5.26014248  4.26652124  1.22689988
## [2263]  6.20335449  5.54240980  5.06659025  4.73335192  1.88749309  6.24075228
## [2269]  6.63929252  5.28472285  5.08413148 -0.34043148  8.06624563  6.37921853
## [2275]  5.79221682  5.82949250  5.91789973  4.76375209  4.33607430  6.92013110
## [2281]  6.47317569  7.41396151  2.33002552  2.42171303  6.51674400  6.25873495
## [2287]  5.09147056  3.20468985  5.44261638  6.65799195 -5.92823530  5.64960299
## [2293]  6.70129118 -1.46428840  2.07200301  6.56928686  5.33960940  4.78877321
## [2299]  6.07167606  5.44091686  6.86042187  6.93020859  5.38765963  6.00714320
## [2305]  5.65374841 -6.36002502 -6.22533069  0.53410071  5.56881143  5.18949821
## [2311] -7.21446911 -0.58068884  7.47199819  7.39529939 -0.86438146  1.97006319
## [2317]  5.75868256  7.90654354  5.99846802  7.13916860  2.47732907  5.92938689
## [2323]  5.89355531 -5.74347531  6.53623321  1.79460547  7.20728126  6.57634012
## [2329]  5.99073060  4.63418703  4.97505899  5.86335387  1.86241123  5.96064122
## [2335]  5.76530425  4.03694187  6.16823234  4.94276527 -1.82884488  4.78803055
## [2341]  5.10274109  6.93685572 -5.94131444  6.61560084  5.76576332  7.24963528
## [2347]  5.14610452  5.11365249  6.58303067  5.42760981  1.71315277  1.84710058
## [2353]  6.87445519  6.60226564  2.55604529  3.85619604 -0.47129366  6.26323015
## [2359]  5.47771712  7.86897105  5.86684483  5.57795679  6.26378401  6.01850638
## [2365]  6.06529524  6.09268829  6.75066738  6.58644812  5.14781873  6.16218554
## [2371]  6.11871481 -2.27094751  7.69939705 -6.64013467  5.99713901  5.58034883
## [2377]  6.62067831  5.30548717 -4.94969094  7.41903706  6.87163877  5.72235503
## [2383]  5.90773745 -1.17534070  7.82468888  5.48151849  7.12794621  5.27230509
## [2389] -6.06723500 -3.91247851  0.90838362  6.01012332  7.10457153  6.50789736
## [2395] -0.46283086  5.81776661  7.47170305  6.70501847 -3.62484937 -1.77873585
## [2401] -0.45577661  4.62381669  0.23543375  6.87054976  5.97518891  4.53486125
## [2407]  6.38124803 -6.69175388 -1.03155327  5.82973033  5.96299923  2.73319500
## [2413]  7.72871980  8.08797934  5.43224844  6.05949863  2.92602748  6.99452811
## [2419]  1.71494617  5.82245746 -3.45502914  6.26342601  5.93890456  6.46693495
## [2425]  5.47236588  2.49379831  7.25952053  5.86013419  2.63041334  5.45477416
## [2431]  3.07931274  6.96710288  6.43022117 -7.40836772 -6.00764725  0.46910090
## [2437]  4.61920703  5.75161849  5.96408233 -1.17828792  6.70712749  6.40130627
## [2443]  5.90028866  6.77314447  7.46279944  2.11916798  4.23498573  3.89249490
## [2449]  5.41208559  5.35803591 -4.10721685  6.24211173  5.68487907  6.42437649
## [2455]  6.30140422  6.88592144  6.48698981  3.96975734 -6.04194962  6.69342474
## [2461]  6.19494256  8.02497351  6.66206742  2.70464550 -5.17026931  6.76802060
## [2467]  6.02397301  6.33677792  0.51111259  0.68868721  5.48030067  6.72552590
## [2473]  7.80237389  7.05506980  6.85677971  6.97922924  2.08703587  7.07409779
## [2479]  6.57917525  4.08526922  6.68027928  1.50187980  5.96287365  6.65730723
## [2485]  6.20828202  7.15833792 -4.86681724  6.68905993  6.61278209 -1.09563587
## [2491]  6.32060097  6.87265907  6.21731468  7.15545238  6.41589572  6.46750636
## [2497]  5.50856860  1.87777895 -5.83908724 -0.44717972  6.15521850  8.16654514
## [2503]  7.44367768  5.86439382  6.27471129  6.37243452  6.41396934  2.68208843
## [2509]  1.42103837 -7.30406852  7.28432506  5.84799431  6.19509917  4.67167139
## [2515]  6.18224667  1.91596951  7.89395676  3.18755211  6.61303755  6.11494182
## [2521]  6.40535801  4.64533017  6.69434257  5.84159609  8.03796459  4.07175313
## [2527]  6.70521451 -2.80665924  5.54217668  5.44180998  6.26030227 -0.47111109
## [2533]  5.18941027 -4.63901394  5.19246802  6.11859229  6.43793604  6.58709405
## [2539]  7.78239360 -1.26785470  6.74142595 -4.35106074  5.20729423  6.64424814
## [2545]  6.00323388  7.20094651  5.75256624 -3.10091034  5.29540840  6.66182010
## [2551]  5.07727580  5.43261067  6.31302422 -4.64395229 -4.37852838  6.90061220
## [2557]  6.10726780  5.18009955  2.84852178 -6.41045939  5.65187964 -6.74024865
## [2563]  2.44784270  5.20790202  7.25749593  6.50226354  5.56239457 -1.02370088
## [2569]  6.55354380  5.15030978  6.79477782 -0.36573190  6.08227778  6.69008570
## [2575]  6.28920033  6.20252123 -0.78157865  4.45361479 -1.33858184  1.20571637
## [2581] -0.02781695  6.13913580  5.08107117  6.11583537  4.65553675  4.43779128
## [2587]  5.93235233  0.74675712  8.04653460  6.98807620  4.97424753  5.63089139
## [2593] -4.31065989 -6.01524233  6.41915335  6.78972767  5.71998169  6.92675630
## [2599]  6.86297847  7.93713400  7.19660516  5.85743041  5.98537887  6.89554448
## [2605] -0.61084129  5.59380162  0.66117026  7.82742615  6.76200061  6.34201139
## [2611]  6.89368860  5.17938473  6.61426583  6.13956568  6.38544287  0.15321392
## [2617]  6.26725004  2.23266637  5.14302085  5.57880857
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = parqs$Acidentes.y, listw = nb2listw(dnb, style = "B"), 
##     zero.policy = NULL, spChk = NULL, return_internals = FALSE)
## attr(,"class")
## [1] "localG"

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)

##            Ii          E.Ii       Var.Ii      Z.Ii Pr(z != E(Ii))
## 1 -0.02435956 -1.960189e-04 1.612806e-05 -6.016854   1.778389e-09
## 2 -0.02880919 -2.783009e-04 2.528074e-05 -5.674407   1.391700e-08
## 3  0.01517515 -1.737105e-05 2.021197e-05  3.379290   7.267331e-04
## 4 -0.03024888 -3.089244e-04 3.232733e-05 -5.265823   1.395627e-07
## 5 -0.02803859 -2.492755e-04 3.045419e-05 -5.035637   4.762628e-07
## 6 -0.02338610 -2.783009e-04 1.830981e-04 -1.707720   8.768839e-02

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)

A caption

mapview(parqs, zcol = "p_valor")

temp = parqs[parqs$p_valor != "[0.05,0.1]", ]
mapview(temp, zcol = "p_valor")

Montando matrix W de vizinhança

ShapeCG.nb1.mat <- nb2mat(dnb)

Incidência de acidentes padronizada

Acidentes_SD <- scale(parqs$Acidentes.y)

Média das incidências de acedentes padronizada

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

Diagrama de espalhamento de Moran

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

A caption

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

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)

A caption

CORES.5[Q][1:5]

## [1] "yellow" "yellow" "green"  "yellow" "yellow"

head(CORES.5[Q])

## [1] "yellow" "yellow" "green"  "yellow" "yellow" "gray"

#save(parqs, file = "parqsFinal-975.Rds")

parqs$cores5 =  ifelse(parqs$cores5Q == "blue", "A-A", ifelse(parqs$cores5Q == "green", "B-B", 
            ifelse(parqs$cores5Q == "red", "A-B", ifelse(parqs$cores5Q == "yellow", "B-A", "NA"))))
mapview(parqs, zcol = "cores5")

temp = parqs[parqs$cores5 == "A-A", ]
mapview(temp, zcol = "cores5")

mapview(temp, zcol = "Acidentes.x", col.regions=brewer.pal(9, "YlOrRd"))