Abstract

This file contains all the code for reproducing the figures and tables for assessing childhood vaccination coverage in Kenya as discussed in Section 11.2.

library(bmstdr)
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
library(ggplot2)
library(rgdal)
library(spdep)
library(broom)
library(tidyr)
library(RColorBrewer)
library(ggpubr)
library(doBy)
library(sf)

# Download mapfiles and data files from https://github.com/sujit-sahu/bookbmstdr
## Note the start time
start.time<-proc.time()[3]
figurepath <- "figures/chap11figures/kenya/"
filepath <- "txttables/"

if (!file.exists(figurepath)) {
  dir.create(figurepath)
}  

if (!file.exists(filepath)) {
  dir.create(filepath)
}  

mappath <- "mapfiles"
dpath <- "datafiles/" # Folder containing the data files 
longrun <- TRUE  # 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  



Kmap <- readOGR(dsn=path.expand(mappath), layer="sdr_subnational_boundaries2")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/sks/Dropbox/sks/website/bookbmstdr/mapfiles", layer: "sdr_subnational_boundaries2"
## with 47 features
## It has 27 fields
kdat <- read.csv(file=paste0(dpath, "C11Kenya_vaccine.csv"))
kdat$vacprop <- kdat$yvac/kdat$n
WKenya <- read.csv(file=paste0(dpath,"C11WKenya.csv"))

1 Code to reproduce Figures in Section 11.2


head(kdat)
##   id      countyname   n yvac propprimary traveltime nightlight wted_n
## 1  1         baringo  92   75       0.245  42.290471 0.02899275     92
## 2  2           bomet  83   76       0.257  30.461168 0.07167645     83
## 3  3         bungoma 120  103       0.198  45.273968 0.08740881    120
## 4  4           busia  76   69       0.161   6.795021 0.14083326     76
## 5  5 elgeyo marakwet  77   67       0.286  27.275333 0.05260377     77
## 6  6            embu  54   50       0.292  22.485116 0.08914773     54
##   wted_yvac   vacprop
## 1        76 0.8152174
## 2        76 0.9156627
## 3       101 0.8583333
## 4        70 0.9078947
## 5        67 0.8701299
## 6        50 0.9259259
adf <- tidy(Kmap)
## Regions defined for each Polygons
adf$id <- as.numeric(adf$id)
kdat$id <- kdat$id-1
udf <- merge(kdat, adf)
head(udf)
##   id countyname  n yvac propprimary traveltime nightlight wted_n wted_yvac
## 1  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 2  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 3  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 4  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 5  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 6  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
##     vacprop     long      lat order  hole piece group
## 1 0.8152174 35.78413 1.652484     1 FALSE     1   0.1
## 2 0.8152174 35.78468 1.652248     2 FALSE     1   0.1
## 3 0.8152174 35.78502 1.652270     3 FALSE     1   0.1
## 4 0.8152174 35.78550 1.652377     4 FALSE     1   0.1
## 5 0.8152174 35.78586 1.652591     5 FALSE     1   0.1
## 6 0.8152174 35.78628 1.652847     6 FALSE     1   0.1

spatframe <- ggplot(data=adf, aes(x=long, y=lat, group = group, fill="White")) +
  # scale_fill_gradientn(colours=com, na.value="black",limits=range(bdf$no2))  +
  geom_path(colour='black',size=0.5)  +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "A map of 47 Counties in Kenya",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())

spatframe


head(adf)
## # A tibble: 6 × 7
##    long   lat order hole  piece group    id
##   <dbl> <dbl> <int> <lgl> <fct> <fct> <dbl>
## 1  35.8  1.65     1 FALSE 1     0.1       0
## 2  35.8  1.65     2 FALSE 1     0.1       0
## 3  35.8  1.65     3 FALSE 1     0.1       0
## 4  35.8  1.65     4 FALSE 1     0.1       0
## 5  35.8  1.65     5 FALSE 1     0.1       0
## 6  35.8  1.65     6 FALSE 1     0.1       0
udf <- merge(kdat, adf)
head(udf)
##   id countyname  n yvac propprimary traveltime nightlight wted_n wted_yvac
## 1  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 2  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 3  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 4  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 5  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
## 6  0    baringo 92   75       0.245   42.29047 0.02899275     92        76
##     vacprop     long      lat order  hole piece group
## 1 0.8152174 35.78413 1.652484     1 FALSE     1   0.1
## 2 0.8152174 35.78468 1.652248     2 FALSE     1   0.1
## 3 0.8152174 35.78502 1.652270     3 FALSE     1   0.1
## 4 0.8152174 35.78550 1.652377     4 FALSE     1   0.1
## 5 0.8152174 35.78586 1.652591     5 FALSE     1   0.1
## 6 0.8152174 35.78628 1.652847     6 FALSE     1   0.1

a <- range(udf$vacprop)
vmap <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=vacprop)) +
   scale_fill_gradientn(colours=colpalette, na.value="black",limits=a)  +
  # geom_path(colour='black',size=0.5)  +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Observed vaccination map of 47 Counties in Kenya",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
vmap


ggsave(filename=paste0(figurepath, "Observed_vaccination_rate.png"))
## Saving 7 x 5 in image


a <- range(udf$nightlight)
nightmap <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=nightlight)) +
  scale_fill_gradientn(colours=colpalette, na.value="black",limits=a)  +
  # geom_path(colour='black',size=0.5)  +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Nightlight map of 47 Counties in Kenya",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
nightmap

ggsave(filename=paste0(figurepath, "nightlightmap.png"))
## Saving 7 x 5 in image

a <- range(udf$propprimary)
primarymap <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=propprimary)) +
  scale_fill_gradientn(colours=colpalette, na.value="black",limits=a)  +
  # geom_path(colour='black',size=0.5)  +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Proportion of primary educated women in Kenya",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
primarymap

ggsave(filename=paste0(figurepath, "primarymap.png"))
## Saving 7 x 5 in image

a <- range(udf$traveltime)
traveltimemap <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=traveltime)) +
  scale_fill_gradientn(colours=colpalette, na.value="black",limits=a)  +
  # geom_path(colour='black',size=0.5)  +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Travel time map",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
traveltimemap


ggsave(filename=paste0(figurepath, "traveltimemap.png"))
## Saving 7 x 5 in image


# Standardize the covariates 

for (i in 5:7) {
  kdat[, i] <- (kdat[, i] - mean( kdat[, i] ))/sd(kdat[, i] )
}

# Construct the adjacency matrix

# head(Kmap@data)
sf::sf_use_s2(FALSE)
## Spherical geometry (s2) switched off
nbhood <- poly2nb(Kmap)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
# nbhood
W <- nb2mat(nbhood, style = "B")
dim(W)
## [1] 47 47
## Read W


dim(WKenya)
## [1] 47 48
countyids <- WKenya[, 1]
table(countyids)
## countyids
##  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
WKenya <- as.matrix(WKenya[, -1])
a <- apply(WKenya, 1, sum)
table(a)
## a
##  1  2  3  4  5  6  7  8 
##  1  2  8  9 13  3  6  5
dim(WKenya)
## [1] 47 47
class(WKenya)
## [1] "matrix" "array"

# Check 
summary(c(W-WKenya))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0       0       0       0       0


N <- 220000
burn.in <- 20000
thin <- 100
# Spatial Binomial GLM

head(kdat)
##   id      countyname   n yvac propprimary  traveltime nightlight wted_n
## 1  0         baringo  92   75   0.2416790 -0.14231623 -0.2566808     92
## 2  1           bomet  83   76   0.3771730 -0.38350771 -0.2273954     83
## 3  2         bungoma 120  103  -0.2890058 -0.08148476 -0.2166014    120
## 4  3           busia  76   69  -0.7067790 -0.86604443 -0.1799468     76
## 5  4 elgeyo marakwet  77   67   0.7046169 -0.44846474 -0.2404812     77
## 6  5            embu  54   50   0.7723639 -0.54613402 -0.2154083     54
##   wted_yvac   vacprop
## 1        76 0.8152174
## 2        76 0.9156627
## 3       101 0.8583333
## 4        70 0.9078947
## 5        67 0.8701299
## 6        50 0.9259259
colnames(kdat)
##  [1] "id"          "countyname"  "n"           "yvac"        "propprimary"
##  [6] "traveltime"  "nightlight"  "wted_n"      "wted_yvac"   "vacprop"
f1 <-  yvac ~ propprimary + traveltime + nightlight

## independent logistic regression
M1 <- Bcartime(formula=f1,   data=kdat, family="binomial", trials="n",
               N=N, burn.in = burn.in, thin=thin, verbose=FALSE)
## ##
## # Total time taken:: 5.49  - Sec.
summary(M1)
## 
##  The  glm  model has been fitted using the  CARBayes  package.
## Call:
## Bcartime(formula = f1, data = kdat, family = "binomial", trials = "n", 
##     N = N, burn.in = burn.in, thin = thin, verbose = FALSE)
## ##
## # Total time taken:: 5.49  - Sec.
## Model formula
## yvac ~ propprimary + traveltime + nightlight
## 
## 
## Parameter Estimates:
##               Mean   2.5%  97.5% n.sample % accept n.effective Geweke.diag
## (Intercept)  1.904  1.809  2.002     2000     47.9        2000        -0.8
## propprimary  0.459  0.347  0.575     2000     47.9        2000         0.9
## traveltime  -0.121 -0.221 -0.019     2000     47.9        2000         0.3
## nightlight   0.102 -0.012  0.227     2000     47.9        2000        -0.6
## 
## Model Choice Statistics:
##           DIC           p.d          WAIC           p.w          LMPL 
##       326.618         3.972       332.201         8.693      -166.420 
## loglikelihood 
##      -159.337
# plot(M1)
# Leroux model
M2 <- Bcartime(formula=f1, data=kdat, scol="id", W=WKenya,
                      family="binomial", trials="n", model="leroux",
                      N=N, burn.in = burn.in, thin=thin, verbose=FALSE)
## ##
## # Total time taken:: 21.91  - Sec.
summary(M2)
## 
##  The  leroux  model has been fitted using the  CARBayes  package.
## Call:
## Bcartime(formula = f1, data = kdat, family = "binomial", trials = "n", 
##     scol = "id", model = "leroux", W = WKenya, N = N, burn.in = burn.in, 
##     thin = thin, verbose = FALSE)
## ##
## # Total time taken:: 21.91  - Sec.
## Model formula
## yvac ~ propprimary + traveltime + nightlight
## 
## 
## Parameter Estimates:
##               Mean   2.5%  97.5% n.sample % accept n.effective Geweke.diag
## (Intercept)  2.030  1.922  2.144     2000     48.5      2000.0         2.3
## propprimary  0.314  0.049  0.573     2000     48.5       856.1        -0.1
## traveltime  -0.317 -0.597 -0.063     2000     48.5       945.6        -1.2
## nightlight   0.017 -0.198  0.230     2000     48.5      1484.9         0.4
## tau2         0.533  0.226  1.041     2000    100.0      1797.9         0.1
## rho          0.489  0.067  0.927     2000     52.5      1809.6        -0.6
## 
## Model Choice Statistics:
##           DIC           p.d          WAIC           p.w          LMPL 
##       263.603        30.713       261.496        21.458      -140.445 
## loglikelihood 
##      -101.088

u <- M2$params[, 1:3]
v <- as.data.frame(M1$params[, 1:3])
dim(u)
## [1] 6 3
dim(v)
## [1] 4 3
k <- nrow(v)
u
##                Mean    2.5%   97.5%
## (Intercept)  2.0299  1.9218  2.1441
## propprimary  0.3139  0.0493  0.5726
## traveltime  -0.3174 -0.5973 -0.0633
## nightlight   0.0166 -0.1983  0.2295
## tau2         0.5331  0.2256  1.0408
## rho          0.4887  0.0666  0.9271
v
##                Mean    2.5%   97.5%
## (Intercept)  1.9042  1.8090  2.0020
## propprimary  0.4593  0.3473  0.5751
## traveltime  -0.1209 -0.2208 -0.0190
## nightlight   0.1025 -0.0122  0.2273
v[k+1, ] <- rep(NA, 3)
v[k+2, ] <- rep(NA, 3)
rownames(v) <- rownames(u)
w <- cbind.data.frame(v, u) # Now repair
dim(w)
## [1] 6 6
w
##                Mean    2.5%   97.5%    Mean    2.5%   97.5%
## (Intercept)  1.9042  1.8090  2.0020  2.0299  1.9218  2.1441
## propprimary  0.4593  0.3473  0.5751  0.3139  0.0493  0.5726
## traveltime  -0.1209 -0.2208 -0.0190 -0.3174 -0.5973 -0.0633
## nightlight   0.1025 -0.0122  0.2273  0.0166 -0.1983  0.2295
## tau2             NA      NA      NA  0.5331  0.2256  1.0408
## rho              NA      NA      NA  0.4887  0.0666  0.9271
dics <- round(c(M1$mchoice[c(2, 1)], NA , M2$mchoice[c(2, 1)], NA), 3)
waics <- round(c(M1$mchoice[c(4, 3)], NA , M2$mchoice[c(4, 3)], NA), 3)
k <- nrow(w)
w[k+1, ] <- dics
w[k+2, ] <- waics

rownames(w)[c(k+1, k+2)] <- c("DIC", "WAIC")
table11.1 <- w
round(table11.1, 3)
##               Mean    2.5%  97.5%   Mean    2.5%  97.5%
## (Intercept)  1.904   1.809  2.002  2.030   1.922  2.144
## propprimary  0.459   0.347  0.575  0.314   0.049  0.573
## traveltime  -0.121  -0.221 -0.019 -0.317  -0.597 -0.063
## nightlight   0.102  -0.012  0.227  0.017  -0.198  0.230
## tau2            NA      NA     NA  0.533   0.226  1.041
## rho             NA      NA     NA  0.489   0.067  0.927
## DIC          3.972 326.618     NA 30.713 263.603     NA
## WAIC         8.693 332.201     NA 21.458 261.496     NA
dput(table11.1, file=paste0(filepath, "table11.1.txt"))

## For generating maps
v <- M2$fit$samples
names(v)
## [1] "beta"   "phi"    "tau2"   "rho"    "fitted" "Y"
dim(v$fitted)
## [1] 2000   47

dim(v$phi)
## [1] 2000   47
x <- get_parameter_estimates(v$phi)
dim(x)
## [1] 47  4
head(x)
##             mean        sd        low         up
## var1 -0.47136249 0.2318861 -0.9230020 -0.0285085
## var2 -0.01432582 0.2930798 -0.5649303  0.5747487
## var3 -0.13377626 0.2462137 -0.6028074  0.3716530
## var4  0.08816217 0.3440417 -0.5787545  0.7891661
## var5 -0.41836042 0.2809500 -0.9770375  0.1240398
## var6  0.25407600 0.3600059 -0.4519933  0.9758183
x$id <- 0:46
udf <- merge(x, adf)
head(udf)
##   id       mean        sd       low         up     long      lat order  hole
## 1  0 -0.4713625 0.2318861 -0.923002 -0.0285085 35.78413 1.652484     1 FALSE
## 2  0 -0.4713625 0.2318861 -0.923002 -0.0285085 35.78468 1.652248     2 FALSE
## 3  0 -0.4713625 0.2318861 -0.923002 -0.0285085 35.78502 1.652270     3 FALSE
## 4  0 -0.4713625 0.2318861 -0.923002 -0.0285085 35.78550 1.652377     4 FALSE
## 5  0 -0.4713625 0.2318861 -0.923002 -0.0285085 35.78586 1.652591     5 FALSE
## 6  0 -0.4713625 0.2318861 -0.923002 -0.0285085 35.78628 1.652847     6 FALSE
##   piece group
## 1     1   0.1
## 2     1   0.1
## 3     1   0.1
## 4     1   0.1
## 5     1   0.1
## 6     1   0.1

reffectsamps <- M2$fit$samples$phi
reffs <- get_parameter_estimates(reffectsamps)
reffs$id <- 0:46
udf <- merge(reffs, adf)
a <- range(udf$mean)

randomeffectmap <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=mean)) +
  scale_fill_gradientn(colours=colpalette, na.value="black",limits=a)  +
  # geom_path(colour='black',size=0.5)  +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Spatial random effects",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
randomeffectmap


ggsave(filename=paste0(figurepath, "randomeffectmap.png"))
## Saving 7 x 5 in image


a <- range(udf$sd)
sdrandomeffectmap <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=sd)) +
  scale_fill_gradientn(colours=colpalette, na.value="black",limits=a)  +
  # geom_path(colour='black',size=0.5)  +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "sd of spatial random effects",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
sdrandomeffectmap


ggsave(filename=paste0(figurepath, "sdrandomeffectmap.png"))
## Saving 7 x 5 in image

a <- as.matrix(M2$fit$samples$fitted)
dim(a)
## [1] 2000   47
ns <- matrix(rep(kdat$n, nrow(a)), byrow = T, nrow=nrow(a))
dim(ns)
## [1] 2000   47
a <- a/ns

propthresh <- function(b, thresh=0.95) {
  length(b[b>thresh])/length(b)
}

head(kdat)
##   id      countyname   n yvac propprimary  traveltime nightlight wted_n
## 1  0         baringo  92   75   0.2416790 -0.14231623 -0.2566808     92
## 2  1           bomet  83   76   0.3771730 -0.38350771 -0.2273954     83
## 3  2         bungoma 120  103  -0.2890058 -0.08148476 -0.2166014    120
## 4  3           busia  76   69  -0.7067790 -0.86604443 -0.1799468     76
## 5  4 elgeyo marakwet  77   67   0.7046169 -0.44846474 -0.2404812     77
## 6  5            embu  54   50   0.7723639 -0.54613402 -0.2154083     54
##   wted_yvac   vacprop
## 1        76 0.8152174
## 2        76 0.9156627
## 3       101 0.8583333
## 4        70 0.9078947
## 5        67 0.8701299
## 6        50 0.9259259
propthresh(a[1,])
## [1] 0.0212766
pcov95 <- apply(a, 2, propthresh)
pcov80 <- apply(a, 2, propthresh, thresh=0.80)

pcovdat <- data.frame(id=0:46, pcov95=pcov95, pcov80=pcov80)
head(pcovdat)
##      id pcov95 pcov80
## var1  0 0.0000 0.8975
## var2  1 0.0135 1.0000
## var3  2 0.0000 0.9730
## var4  3 0.0100 0.9955
## var5  4 0.0000 0.9805
## var6  5 0.2400 1.0000
summary(pcovdat)
##        id           pcov95           pcov80      
##  Min.   : 0.0   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:11.5   1st Qu.:0.0000   1st Qu.:0.7975  
##  Median :23.0   Median :0.0030   Median :0.9955  
##  Mean   :23.0   Mean   :0.1201   Mean   :0.7878  
##  3rd Qu.:34.5   3rd Qu.:0.1502   3rd Qu.:1.0000  
##  Max.   :46.0   Max.   :0.7240   Max.   :1.0000


udf <- merge(pcovdat, adf)
head(udf)
##   id pcov95 pcov80     long      lat order  hole piece group
## 1  0      0 0.8975 35.78413 1.652484     1 FALSE     1   0.1
## 2  0      0 0.8975 35.78468 1.652248     2 FALSE     1   0.1
## 3  0      0 0.8975 35.78502 1.652270     3 FALSE     1   0.1
## 4  0      0 0.8975 35.78550 1.652377     4 FALSE     1   0.1
## 5  0      0 0.8975 35.78586 1.652591     5 FALSE     1   0.1
## 6  0      0 0.8975 35.78628 1.652847     6 FALSE     1   0.1
b <- c(-Inf, 0.25, 0.5, 0.75, 1)
labss <- c("0-0.25", "0.25-0.50", "0.50-0.75", "0.75-1")
udf$col95 <- factor(cut(udf$pcov95, breaks=b,  labels=labss))
udf$col80 <- factor(cut(udf$pcov80, breaks=b,  labels=labss))
table(udf$col95)
## 
##    0-0.25 0.25-0.50 0.50-0.75 
##    285144      8335     53631
table(udf$col80)
## 
##    0-0.25 0.25-0.50 0.50-0.75    0.75-1 
##     43935     22557      6137    274481
com <- c("lightyellow2", "yellow2", "green1", "green4")

cov95map <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=col95)) +
  geom_polygon(colour='black',size=0.25) +
  scale_fill_manual(values =com, guides("Probability"), guide = guide_legend(reverse=TRUE)) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Probability of attaining 95% Coverage",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
cov95map

ggsave(filename =paste0(figurepath, "cov95map.png"))
## Saving 7 x 5 in image

cov80map <- ggplot(data=udf, aes(x=long, y=lat, group = group, fill=factor(col80))) +
  scale_fill_manual(values =com,  guides("Probability"),guide = guide_legend(reverse=TRUE)) +
  geom_polygon(colour='black',size=0.25) +
  theme_bw()+theme(text=element_text(family="Times")) +
  labs(title= "Probability of attaining 80% Coverage",  x="", y = "", size=2.5) +
  theme(axis.text.x = element_blank(), axis.text.y = element_blank(),axis.ticks = element_blank())
cov80map


ggsave(filename =paste0(figurepath, "cov80map.png"))
## Saving 7 x 5 in image

2 Other models not reported



set.seed(5)
vs <- sample(nrow(kdat), 0.1*nrow(kdat))


## BYM model
M1.bym <- Bcartime(formula=f1, data=kdat, scol="id", W=WKenya,
                   family="binomial", trials="n", model="bym",
                   N=N, burn.in = burn.in, thin=thin, verbose = FALSE)
summary(M1.bym)
M1.bym.0 <- Bcartime(formula=f1, data=kdat, scol="id", W=WKenya,
                     family="binomial", trials="n", model="bym",
                     N=N, burn.in = burn.in, thin=thin,
                     validrows = vs, verbose = FALSE)
summary(M1.bym.0)


# 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.96805

Section 11.2: Assessing childhood vaccination coverage in Kenya

Sujit K. Sahu

2022-07-19

1 Code to reproduce Figures in Section 11.2

2 Other models not reported