# Download mapfiles and data files from https://github.com/sujit-sahu/bookbmstdr
## Note the start time
start.time<-proc.time()[3]
figurepath <- "figures/chap11figures/england_respiratory/"
filepath <- "txttables/"
mappath <- "mapfiles" # Folder containing the required map files
dpath <- "datafiles/" # Folder containing the data files
longrun <- FALSE # Should we run the model to check the results.
# If FALSE we use the results, e.g. tables, figures and model fits from previous runs # Load the required data set
load(file=paste0(dpath, "C11England_LA_map_and_Respiratory.RData"))library(bmstdr)
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
library(ggplot2)
library(RColorBrewer)
library(ggpubr)
library(doBy)1 Code to reproduce Table 11.3
The code is shown below. But longrun must be set to TRUE to run the model. Otherwise, it will load the fitted model output. Please see the load statement in this code block.
if (longrun) {
# Takes 1 hour 44 minutes
N <- 120000
burn.in <- 20000
thin <- 100
f0 <- observed ~ offset(logExpected) + no2 + price + jsa + temp.mean
m30 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="localised", G=2, package="CARBayesST",
scol="spaceid", tcol="timeid", N=N,
burn.in = burn.in, thin=thin, verbose=TRUE)
summary(m30)
table11.3 <- m30$params[, 1:3]
dput(table11.3, file=paste0(filepath, "table11.3.txt"))
save(m30, file="Eng_hosp_fits/m30fits.RData")
} else {
table11.3 <- dget(file=paste0(filepath, "table11.3.txt"))
load(file="Eng_hosp_fits/m30fits.RData")
}2 Reproduce Figure 11.7
u <- m30$fit
names(u)
## [1] "summary.results" "samples" "fitted.values"
## [4] "residuals" "modelfit" "accept"
## [7] "localised.structure" "formula" "model"
## [10] "X"
length(u$localised.structure)
## [1] 19380
table(u$localised.structure)
##
## 1 2
## 12877 6503
a <- u$samples
names(a)
## [1] "beta" "lambda" "Z" "delta" "phi" "tau2" "rho.T" "fitted"
dim(a$Z)
## [1] 1000 19380
table(a$Z[1, ])
##
## 1 2
## 12627 6753
v <- a$Z
dim(v)
## [1] 1000 19380
table(v[1, ])
##
## 1 2
## 12627 6753
head(v[1:5, 1:10])
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 2 1 1 1 1 1 1 1 1 1
## [2,] 2 1 1 1 1 1 1 1 1 1
## [3,] 2 1 1 1 1 1 1 1 1 1
## [4,] 2 1 1 1 1 1 1 1 1 1
## [5,] 2 1 1 1 1 1 1 1 1 1
dim(a$lambda)
## [1] 1000 2
# ?ST.CARlocalised
Engrespriratory$local <- u$localised.structure
head(Engrespriratory)
## LA Areacode spaceid timeid year month no2 price jsa
## 4081 00EB E06000001 1 1 2007 1 27.93062 11.71738 0.2242934
## 4085 00EB E06000001 1 2 2007 2 40.28149 11.67487 0.2346765
## 4086 00EB E06000001 1 3 2007 3 32.56677 11.71278 0.2570652
## 4087 00EB E06000001 1 4 2007 4 29.35142 11.43111 0.2516304
## 4088 00EB E06000001 1 5 2007 5 22.14810 11.43111 0.1999175
## 4089 00EB E06000001 1 6 2007 6 21.62227 11.71738 0.2242934
## ethnicity temp.mean temp.min temp.max observed expected logExpected
## 4081 98.2 6.084610 0.2182913 9.621449 173 165.4569 5.108711
## 4085 98.2 4.434712 -0.8721246 8.824551 116 175.3086 5.166548
## 4086 98.2 6.106849 1.0370813 10.222726 134 178.6722 5.185553
## 4087 98.2 8.484028 4.6467160 14.319693 108 184.9413 5.220038
## 4088 98.2 11.776571 7.5765405 15.700682 97 185.8099 5.224724
## 4089 98.2 12.940722 9.6615300 16.755153 77 165.4569 5.108711
## smr local
## 4081 1.0455897 1
## 4085 0.6616903 1
## 4086 0.7499768 1
## 4087 0.5839692 1
## 4088 0.5220388 1
## 4089 0.4653781 1
getmode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}
# Create the vector with numbers.
v <- c(2,1,2,3,1,2,3,4,1,5,5,3,2,3)
v <- c(2, 2, 1, 1)
table(v)
## v
## 1 2
## 2 2
# Calculate the mode using the user function.
result <- getmode(v)
result
## [1] 2
colnames(Engrespriratory)
## [1] "LA" "Areacode" "spaceid" "timeid" "year"
## [6] "month" "no2" "price" "jsa" "ethnicity"
## [11] "temp.mean" "temp.min" "temp.max" "observed" "expected"
## [16] "logExpected" "smr" "local"
a1 <- summaryBy(local ~ Areacode, data =Engrespriratory, FUN =getmode) # Weekly mean
colnames(a1)[2] <- "local"
head(a1)
## Areacode local
## 1 E06000001 1
## 2 E06000002 2
## 3 E06000003 2
## 4 E06000004 2
## 5 E06000005 1
## 6 E06000006 1
a2 <- summaryBy(local ~ Areacode + year, data =Engrespriratory, FUN =getmode) # Weekly mean
colnames(a2)[3] <- "local"
head(a2) # Reshape
## Areacode year local
## 1 E06000001 2007 1
## 2 E06000001 2008 1
## 3 E06000001 2009 1
## 4 E06000001 2010 1
## 5 E06000001 2011 1
## 6 E06000002 2007 2
a3 <- matrix(a2$local, byrow=T, ncol=5)
head(a3)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 1 1 1 1
## [2,] 2 2 2 2 2
## [3,] 2 2 2 2 2
## [4,] 2 2 2 2 2
## [5,] 1 1 1 1 1
## [6,] 2 2 1 1 1
a3 <- data.frame(a3)
colnames(a3) <- paste0("y", 2007:2011)
a3$Areacode <- unique(a2$Areacode)
head(a3)
## y2007 y2008 y2009 y2010 y2011 Areacode
## 1 1 1 1 1 1 E06000001
## 2 2 2 2 2 2 E06000002
## 3 2 2 2 2 2 E06000003
## 4 2 2 2 2 2 E06000004
## 5 1 1 1 1 1 E06000005
## 6 2 2 1 1 1 E06000006
a4 <- merge(a3, a1)
head(a4)
## Areacode y2007 y2008 y2009 y2010 y2011 local
## 1 E06000001 1 1 1 1 1 1
## 2 E06000002 2 2 2 2 2 2
## 3 E06000003 2 2 2 2 2 2
## 4 E06000004 2 2 2 2 2 2
## 5 E06000005 1 1 1 1 1 1
## 6 E06000006 2 2 1 1 1 1
head(eng323_la_map)
## # A tibble: 6 × 7
## long lat order hole piece group id
## <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>
## 1 548881. 191088. 1 FALSE 1 E09000002.1 E09000002
## 2 548852. 190846. 2 FALSE 1 E09000002.1 E09000002
## 3 548886. 190844. 3 FALSE 1 E09000002.1 E09000002
## 4 548881. 190797. 4 FALSE 1 E09000002.1 E09000002
## 5 549000. 190881. 5 FALSE 1 E09000002.1 E09000002
## 6 549097. 190702. 6 FALSE 1 E09000002.1 E09000002
adf <- merge(eng323_la_map, a4, by.x="id", by.y="Areacode", all.y=T, all.x=F)
head(adf)
## id long lat order hole piece group y2007 y2008 y2009
## 1 E06000001 449068.2 536802.2 1 FALSE 1 E06000001.1 1 1 1
## 2 E06000001 449422.7 536544.7 2 FALSE 1 E06000001.1 1 1 1
## 3 E06000001 449514.0 536430.8 3 FALSE 1 E06000001.1 1 1 1
## 4 E06000001 449662.1 536318.3 4 FALSE 1 E06000001.1 1 1 1
## 5 E06000001 449714.7 536227.0 5 FALSE 1 E06000001.1 1 1 1
## 6 E06000001 450174.1 535934.7 6 FALSE 1 E06000001.1 1 1 1
## y2010 y2011 local
## 1 1 1 1
## 2 1 1 1
## 3 1 1 1
## 4 1 1 1
## 5 1 1 1
## 6 1 1 1
dim(adf)
## [1] 187809 13
adf$y2007 <- factor(adf$y2007)
adf$y2008 <- factor(adf$y2008)
adf$y2009 <- factor(adf$y2009)
adf$y2010 <- factor(adf$y2010)
adf$y2011 <- factor(adf$y2011)
adf$local <- factor(adf$local)
cols <- c("1" = "grey", "2" = "yellow")
localmap <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill=local)) +
scale_fill_manual(values=cols, guides("Local")) +
geom_polygon(colour='black',size=0.25) + #longhurst outlines with half width size
coord_equal() +
theme_bw()+theme(text=element_text(family="Times")) +
labs(title= "Overall, 2007-11", x="", y = "", size=2.5) +
# theme(legend.position=c(0.15, 0.5)) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
plot(localmap)
i <- 2007
lmap2007 <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill=y2007)) +
# scale_fill_gradientn(colours=com,na.value="black",limits=mr) +
scale_fill_manual(values=cols, guides("Local")) +
geom_polygon(colour='black',size=0.25) + #longhurst outlines with half width size
coord_equal() +
theme_bw()+theme(text=element_text(family="Times")) +
labs(title= paste0("Year ", i), x="", y = "", size=2.5) +
theme(legend.position="none") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
lmap2007
i <- 2008
lmap2008 <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill=y2008)) +
# scale_fill_gradientn(colours=com,na.value="black",limits=mr) +
scale_fill_manual(values=cols, guides("Local")) +
geom_polygon(colour='black',size=0.25) + #longhurst outlines with half width size
coord_equal() +
theme_bw()+theme(text=element_text(family="Times")) +
labs(title= paste0("Year ", i), x="", y = "", size=2.5) +
theme(legend.position="none") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
lmap2008
i <- 2009
lmap2009 <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill=y2009)) +
# scale_fill_gradientn(colours=com,na.value="black",limits=mr) +
scale_fill_manual(values=cols, guides("Local")) +
geom_polygon(colour='black',size=0.25) + #longhurst outlines with half width size
coord_equal() +
theme_bw()+theme(text=element_text(family="Times")) +
labs(title= paste0("Year ", i), x="", y = "", size=2.5) +
theme(legend.position="none") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
lmap2009
i <- 2010
lmap2010 <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill=y2010)) +
# scale_fill_gradientn(colours=com,na.value="black",limits=mr) +
scale_fill_manual(values=cols, guides("Local")) +
geom_polygon(colour='black',size=0.25) + #longhurst outlines with half width size
coord_equal() +
theme_bw()+theme(text=element_text(family="Times")) +
labs(title= paste0("Year ", i), x="", y = "", size=2.5) +
theme(legend.position="none") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
lmap2010
i <- 2011
lmap2011 <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill=y2011)) +
# scale_fill_gradientn(colours=com,na.value="black",limits=mr) +
scale_fill_manual(values=cols, guides("Local")) +
geom_polygon(colour='black',size=0.25) + #longhurst outlines with half width size
coord_equal() +
theme_bw()+theme(text=element_text(family="Times")) +
labs(title= paste0("Year ", i), x="", y = "", size=2.5) +
theme(legend.position="none") +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
lmap2011
ggarrange(
lmap2007, lmap2008, lmap2009, lmap2010, lmap2011, localmap,
# common.legend = TRUE, legend = "none",
nrow = 2, ncol = 3
)
ggsave(filename=paste0(figurepath, "local_structures.png"), device = "png")
## Saving 7 x 5 in image3 Model comparison not used in the book
N <- 2000
burn.in <- 1000
thin <- 10
f0 <- observed ~ offset(logExpected) + no2 + price + jsa + temp.mean
m00 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="glm", package="CARBayesST", interaction=T,
N=N, burn.in = burn.in, thin=thin)
summary(m00)
m01 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="anova", package="CARBayesST", interaction=F,
scol="spaceid", tcol="timeid", N=N,
burn.in = burn.in, thin=thin)
summary(m01)
m10 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="anova", package="CARBayesST",
scol="spaceid", tcol="timeid", N=N,
burn.in = burn.in, thin=thin)
summary(m10)
m20 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="ar", package="CARBayesST",
scol="spaceid", tcol="timeid", N=N,
burn.in = burn.in, thin=thin)
summary(m20)
m30 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="localised", G=2, package="CARBayesST",
scol="spaceid", tcol="timeid", N=N,
burn.in = burn.in, thin=thin)
summary(m30)
m40 <- Bcartime(formula = f0, data =Engrespriratory, W = W323, family ="poisson",
model="sepspatial", package="CARBayesST",
scol="spaceid", tcol="timeid", N=N,
burn.in = burn.in, thin=thin)
summary(m40)
m00$mchoice
m01$mchoice
m10$mchoice
m20$mchoice
m30$mchoice
m40$mchoice
# All done
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
## elapsed
## 0.32345