## Note the start time
start.time<-proc.time()[3]
dpath <- "datafiles/" # Folder containing the data files
figurepath <- "figures/chap11figures/london/"
filepath <- "txttables/"
dpath <- "datafiles/" # Folder containing the data files
load(file=paste0(dpath, "C11childpovertylondon.RData"))
# Following MCMC output file is required to reproduce
load(file="mcmc_output_files/London/m13.RData")
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
# if (!longrun) load(file="London/m13.RData")library(bmstdr)
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
library(ggplot2)
require(ggsn)
library(RColorBrewer)
library(doBy)
library(GGally)
library(CARBayes)1 Code to reproduce Figure11.8
The code below can be used to reproduce Figure 11.8
## Get the average fix data
head(lainfo)
## id Areacode Areaname hectares Abbrevname centerx centery
## 1 0 E09000021 Kingston upon Thames 3726.117 Kingston 519297.6 166820.0
## 2 1 E09000008 Croydon 8649.441 Croydon 533290.2 163541.2
## 3 2 E09000006 Bromley 15013.487 Bromley 542895.5 165655.5
## 4 3 E09000018 Hounslow 5658.541 Hounslow 513515.5 175643.2
## 5 4 E09000009 Ealing 5554.428 Ealing 515887.9 181715.5
## 6 5 E09000016 Havering 11445.735 Havering 554049.0 187392.0
head(childpoverty)
## id Areacode LA Areaname Areahectares landareaSqKm
## 201 0 E09000021 00AX Kingston upon Thames 3726 37.25
## 202 0 E09000021 00AX Kingston upon Thames 3726 37.25
## 203 0 E09000021 00AX Kingston upon Thames 3726 37.25
## 204 0 E09000021 00AX Kingston upon Thames 3726 37.25
## 205 0 E09000021 00AX Kingston upon Thames 3726 37.25
## 206 0 E09000021 00AX Kingston upon Thames 3726 37.25
## imdavscore07 year population medianincome econinactive houseprice
## 201 13.1 2006 153700 22400 22.2 249999
## 202 13.1 2007 154500 22900 25.4 290000
## 203 13.1 2008 156000 24050 23.4 285000
## 204 13.1 2009 157300 25200 21.6 270000
## 205 13.1 2010 158600 25600 22.4 300000
## 206 13.1 2011 160400 24800 25.9 300000
## nGPregistration punder16 pall mincome log10price logitpall
## 201 3622 16.1 15.2 -0.330229493 5.397938 -1.719000
## 202 3721 16.8 16.1 -0.259822636 5.462398 -1.650806
## 203 3700 16.2 15.7 -0.097886866 5.454845 -1.680721
## 204 3576 16.1 15.8 0.064048904 5.431364 -1.673185
## 205 3632 15.0 14.9 0.120374389 5.477121 -1.742466
## 206 3950 13.6 13.8 0.007723419 5.477121 -1.832002
## logitinactive
## 201 -1.254049
## 202 -1.077391
## 203 -1.185861
## 204 -1.289131
## 205 -1.242506
## 206 -1.051173
colnames(childpoverty)
## [1] "id" "Areacode" "LA" "Areaname"
## [5] "Areahectares" "landareaSqKm" "imdavscore07" "year"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall" "mincome"
## [17] "log10price" "logitpall" "logitinactive"
u <- aggregate.data.frame(childpoverty[, c(7, 9:15)], by=list(id=childpoverty$id), FUN = mean)
head(lainfo)
## id Areacode Areaname hectares Abbrevname centerx centery
## 1 0 E09000021 Kingston upon Thames 3726.117 Kingston 519297.6 166820.0
## 2 1 E09000008 Croydon 8649.441 Croydon 533290.2 163541.2
## 3 2 E09000006 Bromley 15013.487 Bromley 542895.5 165655.5
## 4 3 E09000018 Hounslow 5658.541 Hounslow 513515.5 175643.2
## 5 4 E09000009 Ealing 5554.428 Ealing 515887.9 181715.5
## 6 5 E09000016 Havering 11445.735 Havering 554049.0 187392.0
lonfixdata <- merge(lainfo, u)
head(lonfixdata)
## id Areacode Areaname hectares Abbrevname centerx centery
## 1 0 E09000021 Kingston upon Thames 3726.117 Kingston 519297.6 166820.0
## 2 1 E09000008 Croydon 8649.441 Croydon 533290.2 163541.2
## 3 10 E09000022 Lambeth 2724.940 Lambeth 530844.7 174356.4
## 4 11 E09000028 Southwark 2991.340 Southwark 533820.6 176665.8
## 5 12 E09000023 Lewisham 3531.706 Lewisham 537670.4 173980.7
## 6 13 E09000011 Greenwich 5044.190 Greenwich 542908.1 176875.7
## imdavscore07 population medianincome econinactive houseprice nGPregistration
## 1 13.10 161472.5 25265 22.95 314599.9 3792.0
## 2 21.31 360603.1 22310 21.44 233495.0 6211.6
## 3 34.94 301793.1 23910 20.27 323925.0 7454.5
## 4 33.33 287420.1 23375 25.18 329775.0 7886.7
## 5 31.04 276252.5 22315 23.27 253029.8 5513.4
## 6 33.94 252620.3 21295 25.64 256749.0 6006.2
## punder16 pall
## 1 14.25 14.04
## 2 24.73 23.95
## 3 31.77 31.63
## 4 31.25 31.11
## 5 30.51 30.02
## 6 29.97 29.48
head(lainfo)
## id Areacode Areaname hectares Abbrevname centerx centery
## 1 0 E09000021 Kingston upon Thames 3726.117 Kingston 519297.6 166820.0
## 2 1 E09000008 Croydon 8649.441 Croydon 533290.2 163541.2
## 3 2 E09000006 Bromley 15013.487 Bromley 542895.5 165655.5
## 4 3 E09000018 Hounslow 5658.541 Hounslow 513515.5 175643.2
## 5 4 E09000009 Ealing 5554.428 Ealing 515887.9 181715.5
## 6 5 E09000016 Havering 11445.735 Havering 554049.0 187392.0
lonfixdata <- lonfixdata [order(as.numeric(lonfixdata$id)), ]
head(lonfixdata)
## id Areacode Areaname hectares Abbrevname centerx centery
## 1 0 E09000021 Kingston upon Thames 3726.117 Kingston 519297.6 166820.0
## 2 1 E09000008 Croydon 8649.441 Croydon 533290.2 163541.2
## 13 2 E09000006 Bromley 15013.487 Bromley 542895.5 165655.5
## 24 3 E09000018 Hounslow 5658.541 Hounslow 513515.5 175643.2
## 28 4 E09000009 Ealing 5554.428 Ealing 515887.9 181715.5
## 29 5 E09000016 Havering 11445.735 Havering 554049.0 187392.0
## imdavscore07 population medianincome econinactive houseprice nGPregistration
## 1 13.10 161472.5 25265 22.95 314599.9 3792.0
## 2 21.31 360603.1 22310 21.44 233495.0 6211.6
## 13 14.36 311275.7 24970 20.33 283400.0 1914.2
## 24 23.20 250287.0 21425 22.40 269389.5 8796.2
## 28 25.10 332975.9 21355 25.19 307380.0 11921.5
## 29 16.07 237488.5 22230 21.60 228175.0 1157.8
## punder16 pall
## 1 14.25 14.04
## 2 24.73 23.95
## 13 16.65 15.98
## 24 24.53 24.33
## 28 25.01 25.12
## 29 19.05 18.08
colnames(lonfixdata)
## [1] "id" "Areacode" "Areaname" "hectares"
## [5] "Abbrevname" "centerx" "centery" "imdavscore07"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall"
b <- c(0, 20, 25, 30, 40, 100)
labss <- c("< 20", "20-25", "25-30", "30-40", "> 40")
lonfixdata$cutpall <- factor(cut(lonfixdata$pall, breaks=b, labels=labss))
lonfixdata$cutpunder16 <- factor(cut(lonfixdata$punder16, breaks=b, labels=labss))
bdf <- merge(lonadf, lonfixdata, by="id")
head(bdf)
## id long lat order hole piece group Areacode Areaname
## 1 0 516401.6 160201.8 1 FALSE 1 0.1 E09000021 Kingston upon Thames
## 2 0 516407.3 160210.5 2 FALSE 1 0.1 E09000021 Kingston upon Thames
## 3 0 516413.3 160217.4 3 FALSE 1 0.1 E09000021 Kingston upon Thames
## 4 0 516419.9 160224.2 4 FALSE 1 0.1 E09000021 Kingston upon Thames
## 5 0 516427.9 160231.3 5 FALSE 1 0.1 E09000021 Kingston upon Thames
## 6 0 516453.8 160253.0 6 FALSE 1 0.1 E09000021 Kingston upon Thames
## hectares Abbrevname centerx centery imdavscore07 population medianincome
## 1 3726.117 Kingston 519297.6 166820 13.1 161472.5 25265
## 2 3726.117 Kingston 519297.6 166820 13.1 161472.5 25265
## 3 3726.117 Kingston 519297.6 166820 13.1 161472.5 25265
## 4 3726.117 Kingston 519297.6 166820 13.1 161472.5 25265
## 5 3726.117 Kingston 519297.6 166820 13.1 161472.5 25265
## 6 3726.117 Kingston 519297.6 166820 13.1 161472.5 25265
## econinactive houseprice nGPregistration punder16 pall cutpall cutpunder16
## 1 22.95 314599.9 3792 14.25 14.04 < 20 < 20
## 2 22.95 314599.9 3792 14.25 14.04 < 20 < 20
## 3 22.95 314599.9 3792 14.25 14.04 < 20 < 20
## 4 22.95 314599.9 3792 14.25 14.04 < 20 < 20
## 5 22.95 314599.9 3792 14.25 14.04 < 20 < 20
## 6 22.95 314599.9 3792 14.25 14.04 < 20 < 20
# 1 = City of London, 2 = Tower Hamlets,
# 3 = Islington, 4 = Camden, 5 = Westminster, 6 = Kensington and Chelsea,
# 7 = Hammersmith and Fulham.
londonframe <- ggplot(data=lonadf, aes(x=long, y=lat, group = group, fill="White")) +
# scale_fill_gradientn(colours=com, na.value="black",limits=range(bdf$no2)) +
geom_path(colour='red',size=0.5) +
geom_text(data=lainfo, inherit.aes=F, fontface="bold", aes(label =Abbrevname, x=centerx, y=centery), size = 3) +
coord_equal() +
theme_bw()+ # theme(text=element_text(family="Times")) +
labs(title= "A map of 32 boroughs and the City of London", x="", y = "", size=3) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank()) +
ggsn::scalebar(data=lonadf, dist =5, location = "topleft", transform=F, dist_unit = "km",
st.dist = .05, st.size =4, height = .06, st.bottom=T)
londonframe
ggsave(filename=paste0(figurepath, "london_frame.png"))
colnames(childpoverty)
## [1] "id" "Areacode" "LA" "Areaname"
## [5] "Areahectares" "landareaSqKm" "imdavscore07" "year"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall" "mincome"
## [17] "log10price" "logitpall" "logitinactive"
dp1 <- childpoverty[, c(16, 17, 19, 18)]
p3 <- dp1 %>% ggpairs(., # mapping = ggplot2::aes(colour="blue"),
lower = list(continuous = wrap("smooth", alpha = 0.3, size=0.1),
discrete = "blank", combo = wrap("box_no_facet", alpha=0.5)),
diag = list(discrete="barDiag",
continuous = wrap("densityDiag", alpha=0.5 )),
upper = list(combo = "blank", discrete = "blank", continuous = wrap("cor", family="sans", size=4, align_percent=0.5))) +
theme(panel.grid.major = element_blank())p3
ggsave(filename=paste0(figurepath, "pairs_logit.png"))
## Saving 7 x 5 in image2 Code for modeling
colnames(childpoverty)
## [1] "id" "Areacode" "LA" "Areaname"
## [5] "Areahectares" "landareaSqKm" "imdavscore07" "year"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall" "mincome"
## [17] "log10price" "logitpall" "logitinactive"
# data transform
childpoverty$mincome <- (childpoverty$medianincome - mean(childpoverty$medianincome))/sd(childpoverty$medianincome)
childpoverty$logitpall <- log(childpoverty$pall/(100-childpoverty$pall))
childpoverty$logitinactive <- log(childpoverty$econinactive/(100-childpoverty$econinactive))
childpoverty$log10price <- log10(childpoverty$houseprice)
childpoverty$id <- as.numeric(childpoverty$id)
childpoverty <- childpoverty[order(childpoverty$id, childpoverty$year), ]
childpoverty$yt <- (childpoverty$year - mean(2006:2015))/sd(2006:2015)
f1 <- logitpall ~ mincome + logitinactive + log10(houseprice) + yt
if (longrun) {
N <- 120000
burn.in <- 20000
thin <- 10
m01 <- Bcartime(formula = f1, data =childpoverty, family ="gaussian",
model="glm", package="CARBayes",
N=N, burn.in = burn.in, thin=thin, verbose=TRUE)
# summary(m01)
m12 <- Bcartime(formula = f1, data =childpoverty, W = Wlon, family ="gaussian",
model="anova", package="CARBayesST", interaction=T,
scol="id", tcol="year", N=N,
burn.in = burn.in, thin=thin, verbose=TRUE)
# summary(m12)
m13 <- Bcartime(formula = f1, data =childpoverty, W = Wlon, family ="gaussian",
model="adaptive", package="CARBayesST",
scol="id", tcol="year", N=N,
burn.in = burn.in, thin=thin, verbose=TRUE)
# summary(m13)
m14 <- Bcartime(formula = f1, data =childpoverty, W = Wlon, family ="gaussian",
model="ar", package="CARBayesST", interaction=T,
scol="id", tcol="year", N=N,
burn.in = burn.in, thin=thin, verbose=TRUE)
# summary(m14)
## Tables 11.4 and 11.5
calerrss <- function(modfit) {
# fits <- modfit$fitted.values
# errorss <- sum((childpoverty$logitpall - fits)^2)
# errorss ## Logit scale
x <- modfit$samples
y <- 100*exp(x$fitted)/(1+exp(x$fitted))
# dim(y)
pfits <- apply(y, 2, mean)
# length(pfits)
u <- childpoverty$pall - pfits
sqrt(errss <- sum(u^2)/330)
}
sd(childpoverty$pall)
ss01 <- calerrss(m01$fit)
ss01
ss12 <- calerrss(m12$fit)
ss12
ss13 <- calerrss(m13$fit)
ss13
ss14 <- calerrss(m14$fit)
ss14
# m12 ANOVA
# m14 AR
# m13 Adaptive
u <- rbind(m01$mchoice, m12$mchoice, m14$mchoice, m13$mchoice)
u <- u[, 1:4]
# u
# class(u)
a <- c(ss01, ss12, ss14, ss13)
table11.4 <- data.frame(u[, c(2, 1, 4, 3)], RMSE=a)
rownames(table11.4) <- c("Independent", "Anova", "AR (1)", "Adaptive")
dput(table11.4, file=paste0(filepath, "table11.4.txt"))
dput(table11.4, file=paste0(filepath, "table11.4.txt"))
table11.5 <- m13$params
dput(table11.5, file=paste0(filepath, "table11.5.txt"))
save(m13, file="London/m13.RData")
} else {
table11.4 <- dget(file=paste0(filepath, "table11.4.txt"))
table11.5 <- dget(file=paste0(filepath, "table11.5.txt"))
}
table11.4
## p.d DIC p.w WAIC RMSE
## Independent 5.909947 235.97242 5.370669 235.58108 6.319141
## Anova 39.374537 51.94284 34.316313 50.52369 4.260647
## AR (1) 118.756900 -14.25834 82.832166 -29.91024 2.767106
## Adaptive 125.252927 -19.93633 85.797060 -37.86067 2.651545
table11.5
## Median 2.5% 97.5% n.sample % accept n.effective
## (Intercept) -1.0021 -3.5321 1.6079 10000 100.0 5066.0
## mincome -0.2691 -0.3255 -0.2104 10000 100.0 6665.3
## logitinactive 0.6175 0.4279 0.8042 10000 100.0 5903.5
## log10(houseprice) 0.1073 -0.3635 0.5610 10000 100.0 5053.4
## yt -0.1035 -0.1478 -0.0590 10000 100.0 7163.9
## nu2 0.0430 0.0340 0.0541 10000 100.0 5290.4
## tau2 0.0631 0.0445 0.0892 10000 100.0 4845.5
## rho.S 0.9677 0.9270 0.9884 10000 45.8 6555.6
## rho.T 0.0537 0.0022 0.2151 10000 100.0 9433.1
## tau2.w 71.2979 31.3615 141.7916 10000 100.0 685.7
## Geweke.diag
## (Intercept) 0.1
## mincome 0.5
## logitinactive -2.6
## log10(houseprice) -0.3
## yt -0.9
## nu2 1.6
## tau2 -0.6
## rho.S 0.4
## rho.T 0.5
## tau2.w -0.7# load(file="London/m13.RData")
modfit <- m13$fit
names(modfit)
## [1] "summary.results" "samples" "fitted.values"
## [4] "residuals" "modelfit" "accept"
## [7] "localised.structure" "formula" "model"
## [10] "X"
colnames(childpoverty)
## [1] "id" "Areacode" "LA" "Areaname"
## [5] "Areahectares" "landareaSqKm" "imdavscore07" "year"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall" "mincome"
## [17] "log10price" "logitpall" "logitinactive" "yt"
fits <- modfit$fitted.values
errorss <- sum((childpoverty$logitpall - fits)^2)
errorss ## Logit scale
## [1] 6.89703
names(modfit)
## [1] "summary.results" "samples" "fitted.values"
## [4] "residuals" "modelfit" "accept"
## [7] "localised.structure" "formula" "model"
## [10] "X"
x <- modfit$samples
names(x)
## [1] "beta" "phi" "rho" "tau2" "nu2" "w" "fitted"
y <- 100*exp(x$fitted)/(1+exp(x$fitted))
dim(y)
## [1] 10000 330
pfits <- apply(y, 2, mean)
length(pfits)
## [1] 330
fullres <- data.frame(childpoverty[, c("id", "Areaname", "year", "pall")], pfits)
head(fullres)
## id Areaname year pall pfits
## 201 0 Kingston upon Thames 2006 15.2 19.67534
## 202 0 Kingston upon Thames 2007 16.1 21.01528
## 203 0 Kingston upon Thames 2008 15.7 20.83288
## 204 0 Kingston upon Thames 2009 15.8 16.21089
## 205 0 Kingston upon Thames 2010 14.9 15.51006
## 206 0 Kingston upon Thames 2011 13.8 17.75272
mean(fullres$pfits[fullres$id==26])
## [1] 43.16635
lafits <- tapply(pfits, childpoverty$id, FUN=mean)
lafits
## 0 1 2 3 4 5 6 7
## 16.50981 21.27845 16.63611 22.49399 24.70637 19.91651 21.77460 20.60261
## 8 9 10 11 12 13 14 15
## 26.34477 22.07402 29.67100 30.92352 30.84192 28.70378 21.10988 29.55818
## 16 17 18 19 20 21 22 23
## 27.26924 22.66004 16.29485 12.01931 20.76152 20.71481 27.63619 28.82647
## 24 25 26 27 28 29 30 31
## 34.55966 34.39498 43.16635 37.46931 38.24909 34.16544 36.12360 32.50646
## 32
## 15.28906
res <- data.frame(lonfixdata[, c("id", "Areaname", "pall")], fits=lafits)
colnames(res)
## [1] "id" "Areaname" "pall" "fits"
# seeres <- res[, c("id", "Areaname", "pall", "fits", "sd", "low", "up")]
# seeres
## find la means
acodes <- childpoverty$id
callameans <- function(x, lapattern=acodes) {
u <- data.frame(x=x, lacodes=lapattern)
v <- summaryBy(x~lacodes, data=u)
as.vector(v$x.mean)
}
u <- callameans(y[1,])
u
## [1] 17.76851 20.44064 15.24474 22.92457 24.73422 20.35880 21.94929 20.02658
## [9] 26.44217 22.55434 28.68263 28.56578 29.22301 27.84435 20.92088 28.76000
## [17] 27.92137 22.72080 17.41364 12.37421 23.09401 19.54338 27.47551 28.42044
## [25] 34.88597 33.73336 41.15856 36.89097 38.30072 34.32172 37.22469 31.70244
## [33] 15.63689
a <- apply(y, 1, callameans)
dim(a)
## [1] 33 10000
lafits <- apply(a, 1, mean)
summary(lafits)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.02 20.76 26.34 26.22 30.92 43.17
lalims <- apply(a, 1, FUN=quantile, probs=c(0.025, 0.975))
lalims
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## 2.5% 15.10546 19.49810 15.27667 20.90036 22.82161 18.34003 20.17642 18.97490
## 97.5% 17.97891 23.08081 18.04579 24.11269 26.65340 21.55020 23.44040 22.31038
## [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16]
## 2.5% 24.42805 20.29135 27.52264 28.79202 28.69879 26.83493 19.35525 27.44566
## 97.5% 28.33113 23.97715 31.86903 33.16502 33.07536 30.60119 22.95610 31.67773
## [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24]
## 2.5% 25.46661 20.83895 14.95480 10.89256 19.01706 19.02246 25.46569 26.69940
## 97.5% 29.14454 24.54151 17.75005 13.24768 22.63646 22.48171 29.85201 31.01365
## [,25] [,26] [,27] [,28] [,29] [,30] [,31] [,32]
## 2.5% 32.23353 32.16612 40.85088 34.99281 35.95388 31.89448 33.69742 30.27964
## 97.5% 36.94225 36.66839 45.47938 39.95876 40.59077 36.50162 38.67552 34.72702
## [,33]
## 2.5% 13.73232
## 97.5% 17.02120
sdfit <- apply(a, 1, sd)
#head(lainfo)
fits <- data.frame(id=as.character(0:32), fits=lafits, sd=sdfit, low=lalims[1, ], up=lalims[2, ] )
head(fits)
## id fits sd low up
## 1 0 16.50981 0.7392087 15.10546 17.97891
## 2 1 21.27845 0.9218175 19.49810 23.08081
## 3 2 16.63611 0.7085619 15.27667 18.04579
## 4 3 22.49399 0.8171824 20.90036 24.11269
## 5 4 24.70637 0.9804032 22.82161 26.65340
## 6 5 19.91651 0.8152407 18.34003 21.55020
colnames(lonfixdata)
## [1] "id" "Areacode" "Areaname" "hectares"
## [5] "Abbrevname" "centerx" "centery" "imdavscore07"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall" "cutpall"
## [17] "cutpunder16"
res <- merge(lonfixdata, fits)
class(lonfixdata$id)
## [1] "character"
class(fits$id)
## [1] "character"
# head(res)
dim(res)
## [1] 33 21
res <- res[order(as.numeric(res$id)), ]
colnames(res)
## [1] "id" "Areacode" "Areaname" "hectares"
## [5] "Abbrevname" "centerx" "centery" "imdavscore07"
## [9] "population" "medianincome" "econinactive" "houseprice"
## [13] "nGPregistration" "punder16" "pall" "cutpall"
## [17] "cutpunder16" "fits" "sd" "low"
## [21] "up"
seeres <- res[, c("id", "Areaname", "pall", "fits", "sd", "low", "up")]
b <- c(0, 20, 25, 30, 40, 100)
labss <- c("< 20", "20-25", "25-30", "30-40", "> 40")
res$cutpfits <- factor(cut(res$fits, breaks=b, labels=labss))
table(res$cutpfits)
##
## < 20 20-25 25-30 30-40 > 40
## 6 10 7 9 1
table(res$cutpall)
##
## < 20 20-25 25-30 30-40 > 40
## 9 6 5 12 1
dim(x$w)
## [1] 10000 68
head(x$w[, 1:10])
## Markov Chain Monte Carlo (MCMC) output:
## Start = 1
## End = 7
## Thinning interval = 1
## 2.3 4.5 4.7 5.7 5.8 7.8 5.9
## [1,] 0.9997544 0.9994071 0.9999927 0.9861772 0.9999838 0.9963097 0.9994696
## [2,] 0.9997623 0.9975050 0.9983554 0.9999996 0.9991456 0.9999946 0.9249218
## [3,] 0.9999972 0.9999913 0.9996040 0.9999738 0.9999993 0.3437588 0.8179523
## [4,] 0.9999972 0.9999913 0.9996040 0.9999738 0.9999993 0.3437588 0.8179523
## [5,] 0.9993117 0.9998677 0.9755170 0.8952121 0.9999996 0.9950778 0.2481994
## [6,] 0.9998937 0.9954617 0.9999927 0.6012351 0.9999930 0.9999995 0.9619096
## [7,] 0.9857703 0.9965268 0.9999856 0.8285541 0.9999953 0.9999729 0.9999996
## 8.9 8.10 9.10
## [1,] 0.9995922 0.9985290 0.9999871
## [2,] 0.9988934 0.9999868 0.9999973
## [3,] 0.9999983 0.9999966 0.9999627
## [4,] 0.9999983 0.9999966 0.9999627
## [5,] 0.9999692 0.9999977 0.9999876
## [6,] 0.9901540 0.9977261 0.9999846
## [7,] 0.9987821 0.9999471 0.9969344
colnames(x$w)
## [1] "2.3" "4.5" "4.7" "5.7" "5.8" "7.8" "5.9" "8.9" "8.10"
## [10] "9.10" "2.11" "3.11" "3.12" "11.12" "3.13" "12.13" "3.14" "13.14"
## [19] "3.15" "14.15" "10.16" "16.17" "6.18" "17.18" "1.19" "2.19" "1.20"
## [28] "4.20" "1.21" "2.21" "11.21" "19.21" "1.22" "11.22" "20.22" "21.22"
## [37] "4.23" "5.23" "9.23" "9.24" "23.24" "9.25" "24.25" "9.26" "10.26"
## [46] "25.26" "26.28" "17.29" "27.29" "28.29" "10.30" "16.30" "17.30" "26.30"
## [55] "28.30" "29.30" "17.31" "18.31" "27.31" "29.31" "6.32" "18.32" "31.32"
## [64] "25.33" "26.33" "27.33" "28.33" "29.33"
y <- modfit$localised.structure
names(y)
## [1] "Wmedian" "W99"
Wmed <- modfit$localised.structure$Wmedian
W99 <- modfit$localised.structure$W99
border.locations <- W99
boundary.final <- highlight.borders(border.locations=W99,
spdata=london)
boundary.final
## class : SpatialPoints
## features : 13052
## extent : 507187.7, 549856.1, 159658.1, 198211.5 (xmin, xmax, ymin, ymax)
## crs : NA
hbs <- as.data.frame(boundary.final)
head(hbs)
## X Y
## 1 539596.3 160796.3
## 2 539595.5 160800.1
## 3 539592.8 160816.1
## 4 539590.5 160832.6
## 5 539588.5 160866.6
## 6 539597.4 160883.8
dim(hbs)
## [1] 13052 2
justres <- res[, c("Areacode", "fits", "cutpfits" )]
head(justres)
## Areacode fits cutpfits
## 1 E09000021 16.50981 < 20
## 2 E09000008 21.27845 20-25
## 13 E09000006 16.63611 < 20
## 24 E09000018 22.49399 20-25
## 28 E09000009 24.70637 20-25
## 29 E09000016 19.91651 < 20
bdf <- merge(lonadf, lonfixdata, by="id")
# head(bdf)
ddf <- merge(bdf, justres, by="Areacode")
head(ddf)
## Areacode id long lat order hole piece group Areaname
## 1 E09000001 32 531145.1 180782.1 1 FALSE 1 32.1 City of London
## 2 E09000001 32 531143.8 180799.3 2 FALSE 1 32.1 City of London
## 3 E09000001 32 531143.5 180844.6 3 FALSE 1 32.1 City of London
## 4 E09000001 32 531144.3 180858.7 4 FALSE 1 32.1 City of London
## 5 E09000001 32 531145.6 180880.9 5 FALSE 1 32.1 City of London
## 6 E09000001 32 531142.8 180890.6 6 FALSE 1 32.1 City of London
## hectares Abbrevname centerx centery imdavscore07 population medianincome
## 1 314.942 1 532479.6 181271.8 12.84 7656 57005
## 2 314.942 1 532479.6 181271.8 12.84 7656 57005
## 3 314.942 1 532479.6 181271.8 12.84 7656 57005
## 4 314.942 1 532479.6 181271.8 12.84 7656 57005
## 5 314.942 1 532479.6 181271.8 12.84 7656 57005
## 6 314.942 1 532479.6 181271.8 12.84 7656 57005
## econinactive houseprice nGPregistration punder16 pall cutpall cutpunder16
## 1 38.33 510285 218.6 15.22 15.9 < 20 < 20
## 2 38.33 510285 218.6 15.22 15.9 < 20 < 20
## 3 38.33 510285 218.6 15.22 15.9 < 20 < 20
## 4 38.33 510285 218.6 15.22 15.9 < 20 < 20
## 5 38.33 510285 218.6 15.22 15.9 < 20 < 20
## 6 38.33 510285 218.6 15.22 15.9 < 20 < 20
## fits cutpfits
## 1 15.28906 < 20
## 2 15.28906 < 20
## 3 15.28906 < 20
## 4 15.28906 < 20
## 5 15.28906 < 20
## 6 15.28906 < 20
com <- c("lightyellow2", "lightyellow3" , "yellow2", "firebrick2", "purple")
lonfitmap <- ggplot(data=ddf, aes(x=long, y=lat, group = group, fill=cutpfits)) +
geom_polygon(colour='black',size=0.25) +
geom_point(data=hbs, aes(x=X, y=Y), inherit.aes = FALSE, size = 0.2, colour = "blue") +
geom_text(data=lainfo, inherit.aes=F, fontface="bold", aes(label =Abbrevname, x=centerx, y=centery), size = 3) +
scale_fill_manual(values =com, guides("Percentage"), guide = guide_legend(reverse=TRUE)) +
theme_bw()+theme(text=element_text(family="Times")) +
coord_equal() +
theme(legend.position = c(0.9, 0.25)) +
labs(title= "Fitted percentage of children living in poverty, 2006-15", x="", y = "", size=2.5) +
theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
lonfitmap
ggsave(filename=paste0(figurepath, "fitted_london_poverty.png"))
## Saving 7 x 5 in image
# 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.3795