First Example
#### Source: https://rpubs.com/bms63/346190
library(CARBayes)
data("housedata", package = "CARBayes")
data("shp", package = "CARBayes")
data("dbf", package = "CARBayes")
# loading pakcages
library("sp");
library("dplyr");
library("tidyr");
library(reshape2);
library(surveillance);
library(spdep);
library(maps);
library(maptools);
library(usmap);
library(CARBayes);
library(CARBayesdata);
library(CARBayesST);
library(shapefiles);
library(sp);
library(knitr);
library(RColorBrewer);
library(coda);
library(spam);
library(cdcfluview)
#data(statepov)#data(salesdata)#View(salesdata)
data(lipdata);data(lipdbf);data(lipshp)
kable(summary(lipdata), caption="Summary Statistics")
Summary Statistics
|
Min. : 0.000 |
Min. : 1.100 |
Min. : 0.000 |
Min. :54.94 |
Min. :1.430 |
|
1st Qu.: 4.750 |
1st Qu.: 4.050 |
1st Qu.: 1.000 |
1st Qu.:55.78 |
1st Qu.:3.288 |
|
Median : 8.000 |
Median : 6.300 |
Median : 7.000 |
Median :56.04 |
Median :4.090 |
|
Mean : 9.571 |
Mean : 9.575 |
Mean : 8.661 |
Mean :56.40 |
Mean :4.012 |
|
3rd Qu.:11.000 |
3rd Qu.:10.125 |
3rd Qu.:11.500 |
3rd Qu.:57.02 |
3rd Qu.:4.730 |
|
Max. :39.000 |
Max. :88.700 |
Max. :24.000 |
Max. :60.24 |
Max. :6.800 |
lipdbf$dbf <- lipdbf$dbf[ ,c(2,1)]
data.combined <- combine.data.shapefile(data=lipdata, shp=lipshp, dbf=lipdbf)
W.nb <-poly2nb(data.combined, row.names = rownames(lipdata))
W.mat <- nb2mat(W.nb, style="B")
#str(data.combined)
nb.bound <- poly2nb(data.combined) # shared boundaries
#summary(nb.bound)
coords <- coordinates(data.combined)
#plot(data.combined, border = "gray", main="Scottland")
#plot(nb.bound, coords, pch = 19, cex = 0.6, add = TRUE)
breakpoints <- seq(min(lipdata$observed)-1, max(lipdata$observed)+1, length.out=8)
my.palette <- brewer.pal(n = 7, name = "OrRd")
spplot(data.combined, c("observed", "expected"), main="Scottish Lip Cancer",at=breakpoints,col.regions=my.palette, col="grey")
spplot(data.combined, c("observed", "pcaff"), main="Scottish Lip Cancer", at=breakpoints,col.regions=my.palette, col="black")
glmmodel <- glm(observed~., family="poisson", data=data.combined@data)
#summary(glmmodel)
glmmodel <- glm(observed~., family="poisson", data=data.combined@data)
#summary(glmmodel)
resid.glmmodel <- residuals(glmmodel)
W.nb <- poly2nb(data.combined, row.names = rownames(data.combined@data))
W.list <- nb2listw(W.nb, style="B")
testglm <- moran.mc(x=resid.glmmodel, listw=W.list, nsim=1000)
W <- nb2mat(W.nb, style="B")
formula <- observed ~ expected+pcaff+latitude+longitude
model.spatial1 <- S.CARleroux(formula=formula, data=data.combined@data,family="gaussian", W=W, burnin=20000, n.sample=120000, thin=10, verbose=FALSE)
betas1 <- summarise.samples(model.spatial1$samples$beta, quantiles=c(0.5, 0.025, 0.975))
resultsMS1 <- betas1$quantiles
rownames(resultsMS1) <- c("Intercept", "Expected", "Pcaff", "Latitude", "Longitude")
#kable(resultsMS1, caption="95% Credible Intervals Model 1")
resid.glmmodel <- residuals(glmmodel)
W.nb <- poly2nb(data.combined, row.names = rownames(data.combined@data))
W.list <- nb2listw(W.nb, style="B")
testglm <- moran.mc(x=resid.glmmodel, listw=W.list, nsim=1000)
W <- nb2mat(W.nb, style="B")
formula <- observed ~ expected+pcaff+latitude+longitude
model.spatial1 <- S.CARleroux(formula=formula, data=data.combined@data,family="gaussian", W=W, burnin=20000, n.sample=120000, thin=10, verbose=FALSE)
betas1 <- summarise.samples(model.spatial1$samples$beta, quantiles=c(0.5, 0.025, 0.975))
resultsMS1 <- betas1$quantiles
rownames(resultsMS1) <- c("Intercept", "Expected", "Pcaff", "Latitude", "Longitude")
#kable(resultsMS1, caption="95% Credible Intervals Model 1")
model.spatial2 <- S.CARleroux(formula=formula, data=data.combined@data,family="poisson", W=W, burnin=20000, n.sample=120000, thin=10, verbose=FALSE)
betas2 <- summarise.samples(model.spatial2$samples$beta, quantiles=c(0.5, 0.025, 0.975))
resultsMS2 <- betas2$quantiles
rownames(resultsMS2) <- c("Intercept", "Expected", "Pcaff", "Latitude", "Longitude")
#kable(resultsMS2, caption="95% Credible Intervals Model 2")
model.spatial3 <- S.CARlocalised(formula=formula,G=5, data=data.combined@data,family="poisson", W=W, burnin=20000, n.sample=120000, thin=10, verbose=FALSE)
betas3 <- summarise.samples(model.spatial3$samples$beta, quantiles=c(0.5, 0.025, 0.975))
resultsMS3 <- betas3$quantiles
rownames(resultsMS3) <- c("Expected", "Pcaff", "Latitude", "Longitude")
#kable(resultsMS3, caption="95% Credible Model 3")
getfancy <- matrix(c(model.spatial1$modelfit[1],
model.spatial2$modelfit[1],
model.spatial3$modelfit[1]), nrow = 3, ncol = 1, byrow = TRUE,
dimnames = list(c("Model 1", "Model 2", "Model 3"),
c("DIC")))
getfancy
## DIC
## Model 1 322.3380
## Model 2 318.2261
## Model 3 184.1960
#kable(getfancy, caption="Summary of Deviance")
data("GGHB.IG")
data("pollutionhealthdata")
pollutionhealthdata$SMR <- pollutionhealthdata$observed/ pollutionhealthdata$expected
pollutionhealthdata$logSMR <- log(pollutionhealthdata$SMR)
#par(pty="s", cex.axis=1.5, cex.lab=1.5)
#pairs(pollutionhealthdata[ , c(9, 5:7)], pch = 19, cex = 0.5, lower.panel=NULL, panel=panel.smooth,labels = c("ln(SMR)", "PM10", "JSA", "Price (*100,000)"))
SMR.av <- summarise(group_by(pollutionhealthdata,IG), SMR.mean =mean(SMR))
GGHB.IG@data$SMR <- SMR.av$SMR.mean
price.av <- summarise(group_by(pollutionhealthdata,IG), price.mean =mean(price))
GGHB.IG@data$price <- price.av$price.mean
W.nb <- poly2nb(GGHB.IG, row.names = SMR.av$IG)
W.list <- nb2listw(W.nb, style = "B")
W <- nb2mat(W.nb, style = "B")
nb.bound <- poly2nb(GGHB.IG) # shared boundaries
#summary(nb.bound)
kable(summary(pollutionhealthdata[,c(2:8)]), caption="Summary Statistics")
Summary Statistics
|
Min. :2007 |
Min. : 10.0 |
Min. : 44.47 |
Min. : 7.838 |
Min. : 0.300 |
Min. :0.228 |
Min. :0.2091 |
|
1st Qu.:2008 |
1st Qu.: 54.0 |
1st Qu.: 72.71 |
1st Qu.:10.953 |
1st Qu.: 2.050 |
1st Qu.:0.880 |
1st Qu.:0.6020 |
|
Median :2009 |
Median : 74.0 |
Median : 89.00 |
Median :12.117 |
Median : 3.700 |
Median :1.150 |
Median :0.8458 |
|
Mean :2009 |
Mean : 79.2 |
Mean : 92.35 |
Mean :12.326 |
Mean : 4.191 |
Mean :1.277 |
Mean :0.8650 |
|
3rd Qu.:2010 |
3rd Qu.: 99.0 |
3rd Qu.:109.64 |
3rd Qu.:13.466 |
3rd Qu.: 5.925 |
3rd Qu.:1.541 |
3rd Qu.:1.0820 |
|
Max. :2011 |
Max. :213.0 |
Max. :180.54 |
Max. :19.612 |
Max. :13.750 |
Max. :4.300 |
Max. :2.1871 |
coords <- coordinates(GGHB.IG)
#plot(GGHB.IG, border = "gray", main="Glasgow")
#plot(nb.bound, coords, pch = 19, cex = 0.6, add = TRUE)
par(mfrow=c(1,2))
breakpoints <- seq(min(SMR.av$SMR.mean)-0.1, max(SMR.av$SMR.mean)+0.1,length.out = 11)
spplot(GGHB.IG, "SMR", main="SMR", xlab = "", ylab = "", scales = list(draw = FALSE), at = breakpoints, col.regions = terrain.colors(n = length(breakpoints)-1),par.settings=list(fontsize=list(text=20)))
breakpoints <- seq(min(price.av$price.mean)-0.1, max(price.av$price.mean)+0.1,length.out = 11)
spplot(GGHB.IG, "price", main="Price", xlab = "", ylab = "", scales = list(draw = FALSE), at = breakpoints, col.regions = terrain.colors(n = length(breakpoints)-1),par.settings=list(fontsize=list(text=20)))
formula <- observed ~ offset(log(expected)) + jsa + price + pm10
model1 <- glm(formula = formula, family = "quasipoisson", data = pollutionhealthdata)
resid.glm <- residuals(model1)
Second Example
library(CARBayesdata)
library(shapefiles)
library(sp)
data(lipdata)
data(lipdbf)
data(lipshp)
library(CARBayes)
lipdbf$dbf <- lipdbf$dbf[ ,c(2,1)]
data.combined <- combine.data.shapefile(data=lipdata, shp=lipshp, dbf=lipdbf)
data(GGHB.IG)
data(pricedata)
propertydata.spatial <- merge(x=GGHB.IG, y=pricedata, by="IG", all.x=FALSE)
library(leaflet)
library(rgdal)
propertydata.spatial <- spTransform(propertydata.spatial,
CRS("+proj=longlat +datum=WGS84 +no_defs"))
library(leaflet)
colours <- colorNumeric(palette = "BuPu", domain = propertydata.spatial@data$price)
map1 <- leaflet(data=propertydata.spatial) %>%
addTiles() %>%
addPolygons(fillColor = ~colours(price), color="red", weight=1,
fillOpacity = 0.7) %>%
addLegend(pal = colours, values = propertydata.spatial@data$price, opacity = 1,
title="Price") %>%
addScaleBar(position="bottomleft")
map1
propertydata.spatial@data$logprice <- log(propertydata.spatial@data$price)
propertydata.spatial@data$logdriveshop <- log(propertydata.spatial@data$driveshop)
library(splines)
form <- logprice~ns(crime,3)+rooms+sales+factor(type) + logdriveshop
model <- lm(formula=form, data=propertydata.spatial@data)
library(spdep)
W.nb <- poly2nb(propertydata.spatial, row.names = rownames(propertydata.spatial@data))
W.list <- nb2listw(W.nb, style="B")
resid.model <- residuals(model)
moran.mc(x=resid.model, listw=W.list, nsim=1000)
##
## Monte-Carlo simulation of Moran I
##
## data: resid.model
## weights: W.list
## number of simulations + 1: 1001
##
## statistic = 0.2733, observed rank = 1001, p-value = 0.000999
## alternative hypothesis: greater
library(CARBayes)
W <- nb2mat(W.nb, style="B")
model.spatial <- S.CARleroux(formula=form, data=propertydata.spatial@data,
family="gaussian", W=W, burnin=100000, n.sample=300000, thin=20)
## Setting up the model.
## Generating 10000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 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%
|
|================= | 26%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 29%
|
|==================== | 30%
|
|==================== | 31%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 34%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 37%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 40%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 43%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 46%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 51%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 54%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 57%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 60%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 66%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================== | 70%
|
|============================================== | 71%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 74%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 77%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 80%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 83%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 91%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 94%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 100%
## Summarising results.
## Finished in 80.1 seconds.
##
## #################
## #### Model fitted
## #################
## Likelihood model - Gaussian (identity link function)
## Random effects model - Leroux CAR
## Regression equation - logprice ~ ns(crime, 3) + rooms + sales + factor(type) + logdriveshop
## Number of missing observations - 0
##
## ############
## #### Results
## ############
## Posterior quantities and DIC
##
## Median 2.5% 97.5% n.sample % accept n.effective
## (Intercept) 4.2449 3.9698 4.5292 10000 100.0 10000.0
## ns(crime, 3)1 -0.2456 -0.4046 -0.0963 10000 100.0 9191.7
## ns(crime, 3)2 -0.4047 -0.7062 -0.1145 10000 100.0 8919.3
## ns(crime, 3)3 -0.2026 -0.4104 0.0021 10000 100.0 10875.7
## rooms 0.2195 0.1689 0.2700 10000 100.0 10000.0
## sales 0.0022 0.0016 0.0029 10000 100.0 10350.5
## factor(type)flat -0.2481 -0.3663 -0.1288 10000 100.0 9442.2
## factor(type)semi -0.1631 -0.2654 -0.0611 10000 100.0 10000.0
## factor(type)terrace -0.2939 -0.4245 -0.1680 10000 100.0 10000.0
## logdriveshop -0.0052 -0.0609 0.0515 10000 100.0 8526.0
## nu2 0.0249 0.0150 0.0343 10000 100.0 4303.9
## tau2 0.0416 0.0189 0.0804 10000 100.0 3667.9
## rho 0.9416 0.7640 0.9925 10000 45.3 6155.2
## Geweke.diag
## (Intercept) 0.5
## ns(crime, 3)1 -1.5
## ns(crime, 3)2 -0.8
## ns(crime, 3)3 0.0
## rooms 0.2
## sales -1.7
## factor(type)flat 0.4
## factor(type)semi 0.0
## factor(type)terrace 0.3
## logdriveshop -1.6
## nu2 -0.7
## tau2 -0.1
## rho -0.6
##
## DIC = -143.9527 p.d = 92.56169 LMPL = 56.87
## Length Class Mode
## summary.results 91 -none- numeric
## samples 7 -none- list
## fitted.values 270 -none- numeric
## residuals 2 data.frame list
## modelfit 6 -none- numeric
## accept 5 -none- numeric
## localised.structure 0 -none- NULL
## formula 3 formula call
## model 2 -none- character
## X 2700 -none- numeric
summarise.samples(model.spatial$samples$beta, quantiles=c(0.5, 0.025, 0.975))
## $quantiles
## 0.5 0.025 0.975
## [1,] 4.244881250 3.969815335 4.529218175
## [2,] -0.245578046 -0.404586574 -0.096347539
## [3,] -0.404703786 -0.706225869 -0.114481401
## [4,] -0.202614684 -0.410394134 0.002122644
## [5,] 0.219522004 0.168864815 0.270003683
## [6,] 0.002245438 0.001624014 0.002883847
## [7,] -0.248060580 -0.366279417 -0.128841956
## [8,] -0.163090901 -0.265428870 -0.061093012
## [9,] -0.293941679 -0.424472217 -0.167960876
## [10,] -0.005182311 -0.060875402 0.051505572
##
## $exceedences
## NULL
crime.effect <- summarise.lincomb(model=model.spatial, columns=c(2,3,4),
quantiles=c(0.5, 0.025, 0.975), distribution=FALSE)
plot(propertydata.spatial@data$crime, crime.effect$quantiles[ ,1], pch=19,
ylim=c(-0.55,0.05), xlab="Number of crimes", ylab="Effect of crime")
points(propertydata.spatial@data$crime, crime.effect$quantiles[ ,2], pch=19,
col="red")
points(propertydata.spatial@data$crime, crime.effect$quantiles[ ,3], pch=19,
col="red")
Third Example
data(respiratorydata)
respiratorydata.spatial <- merge(x=GGHB.IG, y=respiratorydata, by="IG", all.x=FALSE)
head(respiratorydata.spatial@data)
## IG name easting northing observed
## 1 S02000260 Auchinairn 261624.5 669657.4 107
## 2 S02000261 Woodhill East 262927.1 670027.8 23
## 3 S02000262 Woodhill West 262142.9 670428.0 53
## 4 S02000263 Westerton East 254570.5 670593.8 40
## 5 S02000264 Bishopbriggs West and Cadder 261248.4 670928.0 60
## 6 S02000265 Westerton West 253764.4 670982.6 25
## expected incomedep SMR
## 1 106.45661 22 1.0051044
## 2 50.97354 7 0.4512145
## 3 104.49236 6 0.5072141
## 4 90.35747 5 0.4426861
## 5 140.16546 7 0.4280655
## 6 63.93549 6 0.3910191
respiratorydata.spatial <- spTransform(respiratorydata.spatial,
CRS("+proj=longlat +datum=WGS84 +no_defs"))
colours <- colorNumeric(palette = "BuPu", domain = respiratorydata.spatial@data$SMR)
map2 <- leaflet(data=respiratorydata.spatial) %>%
addTiles() %>%
addPolygons(fillColor = ~colours(SMR), color="red", weight=1,
fillOpacity = 0.7) %>%
addLegend(pal = colours, values = respiratorydata.spatial@data$SMR, opacity = 1,
title="SMR") %>%
addScaleBar(position="bottomleft")
map2
W.nb <- poly2nb(respiratorydata.spatial, row.names =
rownames(respiratorydata.spatial@data))
W <- nb2mat(W.nb, style="B")
income <- respiratorydata.spatial@data$incomedep
Z.incomedep <- as.matrix(dist(income, diag=TRUE, upper=TRUE))
formula <- observed ~ offset(log(expected))
model.dissimilarity <- S.CARdissimilarity(formula=formula,
data=respiratorydata.spatial@data, family="poisson", W=W,
Z=list(Z.incomedep=Z.incomedep), W.binary=TRUE, burnin=100000,
n.sample=300000, thin=20)
## Setting up the model.
## Generating 10000 post burnin and thinned (if requested) samples.
##
|
| | 0%
|
|= | 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%
|
|================= | 26%
|
|================== | 27%
|
|================== | 28%
|
|=================== | 29%
|
|==================== | 30%
|
|==================== | 31%
|
|===================== | 32%
|
|===================== | 33%
|
|====================== | 34%
|
|======================= | 35%
|
|======================= | 36%
|
|======================== | 37%
|
|========================= | 38%
|
|========================= | 39%
|
|========================== | 40%
|
|=========================== | 41%
|
|=========================== | 42%
|
|============================ | 43%
|
|============================= | 44%
|
|============================= | 45%
|
|============================== | 46%
|
|=============================== | 47%
|
|=============================== | 48%
|
|================================ | 49%
|
|================================ | 50%
|
|================================= | 51%
|
|================================== | 52%
|
|================================== | 53%
|
|=================================== | 54%
|
|==================================== | 55%
|
|==================================== | 56%
|
|===================================== | 57%
|
|====================================== | 58%
|
|====================================== | 59%
|
|======================================= | 60%
|
|======================================== | 61%
|
|======================================== | 62%
|
|========================================= | 63%
|
|========================================== | 64%
|
|========================================== | 65%
|
|=========================================== | 66%
|
|============================================ | 67%
|
|============================================ | 68%
|
|============================================= | 69%
|
|============================================== | 70%
|
|============================================== | 71%
|
|=============================================== | 72%
|
|=============================================== | 73%
|
|================================================ | 74%
|
|================================================= | 75%
|
|================================================= | 76%
|
|================================================== | 77%
|
|=================================================== | 78%
|
|=================================================== | 79%
|
|==================================================== | 80%
|
|===================================================== | 81%
|
|===================================================== | 82%
|
|====================================================== | 83%
|
|======================================================= | 84%
|
|======================================================= | 85%
|
|======================================================== | 86%
|
|========================================================= | 87%
|
|========================================================= | 88%
|
|========================================================== | 89%
|
|========================================================== | 90%
|
|=========================================================== | 91%
|
|============================================================ | 92%
|
|============================================================ | 93%
|
|============================================================= | 94%
|
|============================================================== | 95%
|
|============================================================== | 96%
|
|=============================================================== | 97%
|
|================================================================ | 98%
|
|================================================================ | 99%
|
|=================================================================| 100%
## Summarising results.
## Finished in 697.1 seconds.
print(model.dissimilarity)
##
## #################
## #### Model fitted
## #################
## Likelihood model - Poisson (log link function)
## Random effects model - Binary dissimilarity CAR
## Dissimilarity metrics - Z.incomedep
## Regression equation - observed ~ offset(log(expected))
## Number of missing observations - 0
##
## ############
## #### Results
## ############
## Posterior quantities and DIC
##
## Median 2.5% 97.5% n.sample % accept n.effective
## (Intercept) -0.2197 -0.2416 -0.1982 10000 35.3 9632.4
## tau2 0.1366 0.0977 0.1911 10000 100.0 9476.0
## Z.incomedep 0.0502 0.0465 0.0513 10000 45.3 9626.2
## Geweke.diag alpha.min
## (Intercept) 1.5 NA
## tau2 -0.3 NA
## Z.incomedep -1.6 0.0139
##
## DIC = 1070.285 p.d = 105.2792 LMPL = -611.12
##
## The number of stepchanges identified in the random effect surface
## no stepchange stepchange
## [1,] 261 99
border.locations <- model.dissimilarity$localised.structure$W.posterior
respiratorydata.spatial@data$risk <- model.dissimilarity$fitted.values /
respiratorydata.spatial@data$expected
boundary.final <- highlight.borders(border.locations=border.locations,
spdata=respiratorydata.spatial)
colours <- colorNumeric(palette = "BuPu", domain = respiratorydata.spatial@data$risk)
map3 <- leaflet(data=respiratorydata.spatial) %>%
addTiles() %>%
addPolygons(fillColor = ~colours(risk), color="red", weight=1,
fillOpacity = 0.7) %>%
addLegend(pal = colours, values = respiratorydata.spatial@data$risk, opacity = 1,
title="Risk") %>%
addCircles(lng = ~boundary.final$X, lat = ~boundary.final$Y, weight = 1,
radius = 2) %>%
addScaleBar(position="bottomleft")
map3
#### SOURCE: https://cran.r-project.org/web/packages/CARBayes/vignettes/CARBayes.pdf