## Note the start time
start.time<-proc.time()[3]
# figurepath <- "YourPath" # Folder to save or read the figures from
# filepath <- "YourFilePath" # Folder to save or read the tables from
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)
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" "residuals" "modelfit"
## [6] "accept" "localised.structure" "formula" "model" "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 ethnicity temp.mean temp.min
## 4081 00EB E06000001 1 1 2007 1 27.93062 11.71738 0.2242934 98.2 6.084610 0.2182913
## 4085 00EB E06000001 1 2 2007 2 40.28149 11.67487 0.2346765 98.2 4.434712 -0.8721246
## 4086 00EB E06000001 1 3 2007 3 32.56677 11.71278 0.2570652 98.2 6.106849 1.0370813
## 4087 00EB E06000001 1 4 2007 4 29.35142 11.43111 0.2516304 98.2 8.484028 4.6467160
## 4088 00EB E06000001 1 5 2007 5 22.14810 11.43111 0.1999175 98.2 11.776571 7.5765405
## 4089 00EB E06000001 1 6 2007 6 21.62227 11.71738 0.2242934 98.2 12.940722 9.6615300
## temp.max observed expected logExpected smr local
## 4081 9.621449 173 165.4569 5.108711 1.0455897 1
## 4085 8.824551 116 175.3086 5.166548 0.6616903 1
## 4086 10.222726 134 178.6722 5.185553 0.7499768 1
## 4087 14.319693 108 184.9413 5.220038 0.5839692 1
## 4088 15.700682 97 185.8099 5.224724 0.5220388 1
## 4089 16.755153 77 165.4569 5.108711 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" "month" "no2"
## [8] "price" "jsa" "ethnicity" "temp.mean" "temp.min" "temp.max" "observed"
## [15] "expected" "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 y2010 y2011 local
## 1 E06000001 449068.2 536802.2 1 FALSE 1 E06000001.1 1 1 1 1 1 1
## 2 E06000001 449422.7 536544.7 2 FALSE 1 E06000001.1 1 1 1 1 1 1
## 3 E06000001 449514.0 536430.8 3 FALSE 1 E06000001.1 1 1 1 1 1 1
## 4 E06000001 449662.1 536318.3 4 FALSE 1 E06000001.1 1 1 1 1 1 1
## 5 E06000001 449714.7 536227.0 5 FALSE 1 E06000001.1 1 1 1 1 1 1
## 6 E06000001 450174.1 535934.7 6 FALSE 1 E06000001.1 1 1 1 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.2689667