Section 11.5: Assessing trend in childpoverty in London

Sujit K. Sahu

University of Southampton

2022-03-30

Abstract

This file contains all the code for reproducing the results for the child poverty example in Section 11.5.

## 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 image

2 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