…
## Cases and deaths of poisonning and case fatality rate
## http://sis.statistics.gov.lk/statHtml/statHtml.do?orgId=144&tblId=DT_HEA_ANN_112&conn_path=I2
poison <- read.table(file = "./data/poison.csv",
header = TRUE,
sep = ",", skip = 0, stringsAsFactors = FALSE)
str(poison)
## 'data.frame': 50 obs. of 5 variables:
## $ Period : int 2012 2012 2012 2012 2012 2012 2012 2012 2012 2012 ...
## $ Cases.or.deaths : chr "Cases" "Cases" "Cases" "Cases" ...
## $ District : chr "Colombo" "Gampaha" "Kalutara" "Kandy" ...
## $ Poisoning.by.drugs..medicaments.and.biological.substances : chr "2,670" "3,500" "1,395" "2,813" ...
## $ Toxic.effects.of.Pesticides_Organophosphate.and.carbamate.insecticides: chr "430" "617" "254" "916" ...
poison <- poison[,2:5]
poison.names <- names(poison)
poison.names
## [1] "Cases.or.deaths"
## [2] "District"
## [3] "Poisoning.by.drugs..medicaments.and.biological.substances"
## [4] "Toxic.effects.of.Pesticides_Organophosphate.and.carbamate.insecticides"
names(poison) <- c("cd","DISTRICT","P","T")
poison[,2] <- toupper(poison[,2])
poison.ma <- apply(X = poison[,3:4], MARGIN = 2, FUN = function(x){as.numeric(gsub(",","",x))})
poison.ma
## P T
## [1,] 2670 430
## [2,] 3500 617
## [3,] 1395 254
## [4,] 2813 916
## [5,] 652 682
## [6,] 774 978
## [7,] 1678 192
## [8,] 1053 113
## [9,] 1217 718
## [10,] 718 761
## [11,] 139 158
## [12,] 133 125
## [13,] 212 266
## [14,] 70 125
## [15,] 963 375
## [16,] 1095 514
## [17,] 862 224
## [18,] 3438 2577
## [19,] 1373 508
## [20,] 2188 1560
## [21,] 967 824
## [22,] 1183 1108
## [23,] 784 830
## [24,] 1611 955
## [25,] 1039 279
## [26,] 9 18
## [27,] 12 32
## [28,] 1 10
## [29,] 10 35
## [30,] 5 16
## [31,] 0 13
## [32,] 9 18
## [33,] 6 6
## [34,] 5 14
## [35,] 3 14
## [36,] 0 5
## [37,] 1 2
## [38,] 0 12
## [39,] 0 5
## [40,] 5 3
## [41,] 2 8
## [42,] 3 2
## [43,] 10 67
## [44,] 4 12
## [45,] 1 36
## [46,] 2 41
## [47,] 5 36
## [48,] 4 25
## [49,] 10 24
## [50,] 2 5
poison$P <- poison.ma[,1];poison$T <- poison.ma[,2]
head(poison)
## cd DISTRICT P T
## 1 Cases COLOMBO 2670 430
## 2 Cases GAMPAHA 3500 617
## 3 Cases KALUTARA 1395 254
## 4 Cases KANDY 2813 916
## 5 Cases MATALE 652 682
## 6 Cases NUWERA ELIYA 774 978
poison.c <- poison[poison$cd=="Cases",2:4]
names(poison.c)[2:3] <- c("CP","CT")
head(poison.c)
## DISTRICT CP CT
## 1 COLOMBO 2670 430
## 2 GAMPAHA 3500 617
## 3 KALUTARA 1395 254
## 4 KANDY 2813 916
## 5 MATALE 652 682
## 6 NUWERA ELIYA 774 978
poison.d <- poison[poison$cd=="Deaths",2:4]
names(poison.d)[2:3] <- c("DP","DT")
head(poison.d)
## DISTRICT DP DT
## 26 COLOMBO 9 18
## 27 GAMPAHA 12 32
## 28 KALUTARA 1 10
## 29 KANDY 10 35
## 30 MATALE 5 16
## 31 NUWERA ELIYA 0 13
poison.c <- poison.c[with(poison.c, order(DISTRICT)),]
poison.d <- poison.d[with(poison.d, order(DISTRICT)),]
head(poison.c)
## DISTRICT CP CT
## 16 AMPARA1 1095 514
## 20 ANURADHAPURA 2188 1560
## 22 BADULLA 1183 1108
## 15 BATTICOLOA 963 375
## 1 COLOMBO 2670 430
## 7 GALLE 1678 192
library(rgdal)
districts <- readOGR("./results/","districts")
## OGR data source with driver: ESRI Shapefile
## Source: "./results/", layer: "districts"
## with 25 features
## It has 2 fields
districts$DISTRICT %in% poison.c$DISTRICT
## [1] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [12] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE
## [23] TRUE TRUE TRUE
poison.c$DISTRICT[(districts$DISTRICT %in% poison.c$DISTRICT)==FALSE]
## [1] "AMPARA1" "BATTICOLOA" "MONERAGALA" "NUWERA ELIYA"
districts$DISTRICT[(districts$DISTRICT %in% poison.c$DISTRICT)==FALSE]
## [1] AMPARA BATTICALOA MONARAGALA NUWARA ELIYA
## 25 Levels: AMPARA ANURADHAPURA BADULLA BATTICALOA COLOMBO ... VAVUNIYA
poison.c$DISTRICT[(districts$DISTRICT %in% poison.c$DISTRICT)==FALSE] <- as.character(districts$DISTRICT[(districts$DISTRICT %in% poison.c$DISTRICT)==FALSE])
districts$DISTRICT %in% poison.c$DISTRICT
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
districts$DISTRICT %in% poison.d$DISTRICT
## [1] FALSE TRUE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [12] TRUE TRUE TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE
## [23] TRUE TRUE TRUE
poison.d$DISTRICT[(districts$DISTRICT %in% poison.d$DISTRICT)==FALSE]
## [1] "AMPARA1" "BATTICOLOA" "MONERAGALA" "NUWERA ELIYA"
districts$DISTRICT[(districts$DISTRICT %in% poison.d$DISTRICT)==FALSE]
## [1] AMPARA BATTICALOA MONARAGALA NUWARA ELIYA
## 25 Levels: AMPARA ANURADHAPURA BADULLA BATTICALOA COLOMBO ... VAVUNIYA
poison.d$DISTRICT[(districts$DISTRICT %in% poison.d$DISTRICT)==FALSE] <- as.character(districts$DISTRICT[(districts$DISTRICT %in% poison.d$DISTRICT)==FALSE])
districts$DISTRICT %in% poison.d$DISTRICT
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
library(dplyr)
##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:raster':
##
## intersect, select, union
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
districts@data <- districts@data %>%
left_join(poison.c, by = "DISTRICT")
## Warning in left_join_impl(x, y, by$x, by$y): joining character vector and
## factor, coercing into character vector
districts@data <- districts@data %>%
left_join(poison.d, by = "DISTRICT")
str(districts@data)
## 'data.frame': 25 obs. of 6 variables:
## $ DISTRICT: chr "AMPARA" "ANURADHAPURA" "BADULLA" "BATTICALOA" ...
## $ PROVINCE: Factor w/ 9 levels "CP","EP","NC",..: 2 3 8 2 9 7 9 7 4 9 ...
## $ CP : num 1095 2188 1183 963 2670 ...
## $ CT : num 514 1560 1108 375 430 ...
## $ DP : num 2 1 5 5 9 9 12 5 3 1 ...
## $ DT : num 8 36 36 3 18 18 32 14 14 10 ...
p1 <- spplot(districts,zcol=names(districts@data[3]),main = paste("Cases (2012) of \n",poison.names[3],sep = ""))
p2 <- spplot(districts,zcol=names(districts@data[4]),main = paste("Cases (2012) of \n",poison.names[4],sep = ""))
p3 <- spplot(districts,zcol=names(districts@data[5]),main = paste("Deaths (2012) of \n",poison.names[3],sep = ""))
p4 <- spplot(districts,zcol=names(districts@data[6]),main = paste("Deaths (2012) of \n",poison.names[4],sep = ""))
print(p1, position = c(0,0,.5,1), more = TRUE) # position xmin, ymin, xmax, ymax
print(p3, position = c(.5,0,1,1))
print(p2, position = c(0,0,.5,1), more = TRUE)
print(p4, position = c(.5,0,1,1))
##install.packages("spdep",dependencies = TRUE)
library(spdep)
## Loading required package: Matrix
d.co <- coordinates(districts)
IDs <- row.names(as(districts, "data.frame"))
d.del <- tri2nb(d.co, row.names = IDs)
##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").
d.soi <- graph2nb(soi.graph(d.del, d.co), row.names = IDs)
d.gab <- graph2nb(gabrielneigh(d.co), row.names = IDs)
d.rel <- graph2nb(relativeneigh(d.co), row.names = IDs)
str(d.rel)
## List of 25
## $ : int [1:2] 4 18
## $ : int [1:3] 22 24 25
## $ : int [1:2] 18 20
## $ : int 21
## $ : int [1:2] 7 10
## $ : int [1:2] 10 17
## $ : int 12
## $ : int [1:2] 17 18
## $ : int 13
## $ : int 23
## $ : int [1:2] 16 20
## $ : int [1:2] 14 20
## $ : int 19
## $ : int [1:2] 16 22
## $ : int 25
## $ : int 21
## $ : int 23
## $ : int 0
## $ : int 25
## $ : int 23
## $ : int 24
## $ : int 0
## $ : int 0
## $ : int 0
## $ : int 0
## - attr(*, "region.id")= chr [1:25] "1" "2" "3" "4" ...
## - attr(*, "call")= language relativeneigh(coords = d.co)
## - attr(*, "class")= chr "nb"
## - attr(*, "sym")= logi FALSE
par(mfrow = c(2,2))
plot(districts, main = "Delaunay graph")
plot.nb(d.del, coords = d.co,add=TRUE)
plot(districts, main = "Sphere of influence graph")
plot.nb(d.soi, coords = d.co,add=TRUE)
plot(districts, main = "Gabriel graph")
plot.nb(d.gab, coords = d.co,add=TRUE)
plot(districts, main = "Relative neighbour graph")
plot.nb(d.rel, coords = d.co,add=TRUE)
d.k2 <- knn2nb(knearneigh(d.co, k = 2), row.names = IDs)
d.k5 <- knn2nb(knearneigh(d.co, k = 5), row.names = IDs)
d.d50k <- dnearneigh(d.co, d1=0, d2=50000, row.names = IDs)
d.d200k <- dnearneigh(d.co, d1=0, d2=200000, row.names = IDs)
par(mfrow = c(2,2))
plot(districts, main = "K-nearest neighbour (k = 2)")
plot.nb(d.k2, coords = d.co,add=TRUE)
plot(districts, main = "K-nearest neighbour (k = 5)")
plot.nb(d.k5, coords = d.co,add=TRUE)
plot(districts, main = "Distance based neighbour (0-50km)")
plot.nb(d.d50k, coords = d.co,add=TRUE)
plot(districts, main = "Distance based neighbour (0-200km)")
plot.nb(d.d200k, coords = d.co,add=TRUE)
d.del.w <- nb2listw(d.del, style = "B",zero.policy = TRUE)
d.soi.w <- nb2listw(d.soi, style = "B",zero.policy = TRUE)
d.gab.w <- nb2listw(d.gab, style = "B",zero.policy = TRUE)
d.rel.w <- nb2listw(d.rel, style = "B",zero.policy = TRUE)
d.k2.w <- nb2listw(d.k2, style = "B",zero.policy = TRUE)
d.k5.w <- nb2listw(d.k5, style = "B",zero.policy = TRUE)
d.d50k.w <- nb2listw(d.d50k, style = "B",zero.policy = TRUE)
d.d200k.w <- nb2listw(d.d200k, style = "B",zero.policy = TRUE)
class(d.k5.w)
## [1] "listw" "nb"
summary(d.k5.w)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 25
## Number of nonzero links: 125
## Percentage nonzero weights: 20
## Average number of links: 5
## Non-symmetric neighbours list
## Link number distribution:
##
## 5
## 25
## 25 least connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 with 5 links
## 25 most connected regions:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 with 5 links
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 25 625 125 225 2588
CP.mor.del <- moran.test(x=districts$CP, listw = d.del.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.soi <- moran.test(x=districts$CP, listw = d.soi.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.gab <- moran.test(x=districts$CP, listw = d.gab.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.rel <- moran.test(x=districts$CP, listw = d.rel.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.k2 <- moran.test(x=districts$CP, listw = d.k2.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.k5 <- moran.test(x=districts$CP, listw = d.k5.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.d50k <- moran.test(x=districts$CP, listw = d.d50k.w ,zero.policy=TRUE, alternative = "two.sided")
CP.mor.d200k <- moran.test(x=districts$CP, listw = d.d200k.w ,zero.policy=TRUE, alternative = "two.sided")
library(knitr)
kable(x = t(data.frame(Delaunay=CP.mor.del$estimate,
SOI=CP.mor.soi$estimate,
Gabriel=CP.mor.gab$estimate,
Relative=CP.mor.rel$estimate,
K2=CP.mor.k2$estimate,
K5=CP.mor.k5$estimate,
D50km=CP.mor.d50k$estimate,
D200km=CP.mor.d200k$estimate)),
caption = paste("Moran's I for cases (2012) of ",poison.names[3],"; rows show different methods to determine weights/neighbors",sep = "")
)
Moran’s I for cases (2012) of Poisoning.by.drugs..medicaments.and.biological.substances; rows show different methods to determine weights/neighbors
| Delaunay |
0.3032680 |
-0.0416667 |
0.0119770 |
| SOI |
0.2997814 |
-0.0416667 |
0.0138244 |
| Gabriel |
0.2750768 |
-0.0500000 |
0.0174336 |
| Relative |
0.0928872 |
-0.0526316 |
0.0285042 |
| K2 |
0.2794636 |
-0.0416667 |
0.0292473 |
| K5 |
0.3140812 |
-0.0416667 |
0.0106351 |
| D50km |
0.2014661 |
-0.0526316 |
0.0369753 |
| D200km |
0.0123761 |
-0.0416667 |
0.0009117 |
par(mfrow = c(2,2))
moran.plot(x=districts$CP, listw = d.del.w, main = "Delaunay", spChk = FALSE, zero.policy = TRUE)
moran.plot(x=districts$CP, listw = d.gab.w, main = "Gabriel", spChk = FALSE, zero.policy = TRUE)
moran.plot(x=districts$CP, listw = d.k5.w, main = "K5")
moran.plot(x=districts$CP, listw = d.d200k.w, main = "200km", spChk = FALSE, zero.policy = TRUE)
CP.mor.k5.mc <- moran.mc(x=districts$CP, listw = d.k5.w,zero.policy=TRUE,nsim = 99)
CP.mor.k5.mc
##
## Monte-Carlo simulation of Moran's I
##
## data: districts$CP
## weights: d.k5.w
## number of simulations + 1: 100
##
## statistic = 0.31408, observed rank = 99, p-value = 0.01
## alternative hypothesis: greater
par(mfrow = c(1,1))
tmp <- CP.mor.k5.mc$res[1:length(CP.mor.k5.mc$res)-1]
zz <- density(tmp)
plot(zz,main="Moran’s I Permutation Test",xlab="Reference Distribution", col = 2,ylim=c(0,5),lwd=2)
hist(tmp,freq=F,add=T)
abline(v=CP.mor.k5.mc$statistic,lwd=2,col="green")
CP.mor.del.co <- sp.correlogram(neighbours = d.del, var = districts$CP, order = 5, method = "I",style = "B", zero.policy = TRUE)
CP.mor.soi.co <- sp.correlogram(neighbours = d.soi, var = districts$CP, order = 5, method = "I",style = "B", zero.policy = TRUE)
CP.mor.k5.co <- sp.correlogram(neighbours = d.k5, var = districts$CP, order = 5, method = "I",style = "B", zero.policy = TRUE)
CP.mor.d50k.co <- sp.correlogram(neighbours = d.d50k, var = districts$CP, order = 5, method = "I",style = "B", zero.policy = TRUE)
par(mfrow = c(2,2))
plot(CP.mor.del.co,main = "Delaunay")
plot(CP.mor.soi.co,main = "Sphere of Influence")
plot(CP.mor.k5.co,main = "K5")
plot(CP.mor.d50k.co,main = "50km distance")
