CARBayesST.R
In CARBayesST: Spatio-Temporal Generalised Linear Mixed Models for Areal Unit Data

### R code from vignette source 'CARBayesST.Rnw'

###################################################
### code chunk number 1: CARBayesST.Rnw:92-93
###################################################
options(prompt = "R>  ")


###################################################
### code chunk number 2: CARBayesST.Rnw:537-541
###################################################
library("CARBayesdata")
library("sf")
data("GGHB.IZ")
data("pollutionhealthdata")


###################################################
### code chunk number 3: CARBayesST.Rnw:547-551
###################################################
class(GGHB.IZ)
head(GGHB.IZ)
class(pollutionhealthdata)
head(pollutionhealthdata)


###################################################
### code chunk number 4: CARBayesST.Rnw:556-561
###################################################
library(dplyr)
pollutionhealthdata <- pollutionhealthdata %>% mutate( 
          SMR = pollutionhealthdata$observed / pollutionhealthdata$expected, 
          logSMR = log(pollutionhealthdata$observed / pollutionhealthdata$expected))
head(pollutionhealthdata)


###################################################
### code chunk number 5: CARBayesST.Rnw:567-569
###################################################
library(GGally)
ggpairs(pollutionhealthdata, columns=c(9, 5:7))


###################################################
### code chunk number 6: CARBayesST.Rnw:584-588
###################################################
group_IZ <- group_by(pollutionhealthdata, IZ)
SMR.av <- summarise(group_IZ, SMR.mean = mean(SMR))
GGHB.IZ$SMR <- SMR.av$SMR.mean
head(GGHB.IZ)


###################################################
### code chunk number 7: CARBayesST.Rnw:595-596
###################################################
GGHB.IZ <- st_transform(x=GGHB.IZ, crs='+proj=longlat +datum=WGS84 +no_defs')


###################################################
### code chunk number 8: CARBayesST.Rnw:601-611
###################################################
library(leaflet)
colours <- colorNumeric(palette = "YlOrRd", domain = GGHB.IZ$SMR)
leaflet(data=GGHB.IZ) %>% 
    addTiles() %>% 
    addPolygons(fillColor = ~colours(GGHB.IZ$SMR), 
                color="grey", weight=1, 
                fillOpacity = 0.7) %>%
    addLegend(pal = colours, values = GGHB.IZ$SMR, 
              opacity = 1, title="SMR") %>%
    addScaleBar(position="bottomleft")


###################################################
### code chunk number 9: CARBayesST.Rnw:626-630
###################################################
library("spdep")
W.nb <- poly2nb(GGHB.IZ, row.names = GGHB.IZ$IZ)
W.list <- nb2listw(W.nb, style = "B")
W <- nb2mat(W.nb, style = "B")


###################################################
### code chunk number 10: CARBayesST.Rnw:639-645
###################################################
formula <- observed ~ offset(log(expected)) + jsa + price + pm10
model1 <- glm(formula = formula, family = "quasipoisson", 
    data = pollutionhealthdata)
resid.glm <- residuals(model1)
summary(model1)$coefficients
summary(model1)$dispersion


###################################################
### code chunk number 11: CARBayesST.Rnw:650-651
###################################################
moran.mc(x = resid.glm[1:271], listw = W.list, nsim = 10000)


###################################################
### code chunk number 12: CARBayesST.Rnw:783-788
###################################################
library("CARBayesdata")
library("sf")
data("GGHB.IZ")
data("salesdata")
head(salesdata)


###################################################
### code chunk number 13: CARBayesST.Rnw:794-801
###################################################
salesdata <- salesdata %>% mutate(salesprop = salesdata$sales / salesdata$stock)
library(ggplot2)
ggplot(salesdata, aes(x = factor(year), y = salesprop)) +
    geom_boxplot(fill="red", alpha=0.7) + 
    scale_x_discrete(name = "Year") +
    scale_y_continuous(name = "Sales proportion") + 
    theme(text=element_text(size=16), plot.title=element_text(size=18, face="bold")) 


###################################################
### code chunk number 14: CARBayesST.Rnw:814-819
###################################################
library(dplyr)
group_IZ <- group_by(salesdata, IZ)
salesprop <- summarise(group_IZ, salesproprtion.mean = mean(salesprop))
GGHB.IZ$sales <- salesprop$salesproprtion.mean
head(GGHB.IZ)


###################################################
### code chunk number 15: CARBayesST.Rnw:824-836
###################################################
GGHB.IZ <- st_transform(x=GGHB.IZ, crs='+proj=longlat +datum=WGS84 +no_defs')

library(leaflet)
colours <- colorNumeric(palette = "YlOrRd", domain = GGHB.IZ$sales)
map1 <- leaflet(data=GGHB.IZ) %>% 
    addTiles() %>% 
    addPolygons(fillColor = ~colours(sales), color="grey", weight=1, 
                fillOpacity = 0.7) %>%
    addLegend(pal = colours, values = GGHB.IZ$sales, opacity = 1, 
              title="Sales") %>%
    addScaleBar(position="bottomleft")
map1


###################################################
### code chunk number 16: CARBayesST.Rnw:854-857
###################################################
library("spdep")
W.nb <- poly2nb(GGHB.IZ, row.names = GGHB.IZ$IZ)
W <- nb2mat(W.nb, style = "B")