#Abstract

#Introduction

#Literature Review

Data Analysis

Load Data

library(magrittr)
library(utils)
library(ggplot2)
library(sf)
library(sp)
library(rgdal)
library(dplyr)
library(kableExtra)

### Prepare data
load("cases.RData")

### Observed/expected table (for knitr)
kable(
  cases %>%
    group_by(Year) %>%
    summarise(`Observed` = format(sum(Observed), big.mark=","), `Expected` = format(sum(Expected), big.mark=",")),
  booktabs = TRUE,
  caption = "Aggregated Admissions"
) %>% kable_styling(., full_width = TRUE)
Aggregated Admissions
Year Observed Expected
2003 NA NA
2004 NA NA
2005 NA NA
2006 NA NA
2007 NA NA
2008 NA NA
2009 NA NA
2010 NA NA
2011 NA NA
2012 NA NA
2013 42,919 44,918.5
2014 45,580 44,380.3
2015 50,401 44,339.2
2016 53,565 44,635.6
2017 51,405 43,954.6
2018 57,085 43,190
2019 57,330 42,190 
missing <- data.frame()
missing_LA <- list()
for (i in 2003:2019) {
  o <- cases[cases$Year == i, ]$Observed %>%
    is.na() %>%
    which() %>%
    length()
  e <- cases[cases$Year == i, ]$Expected %>%
    is.na() %>%
    which() %>%
    length()
  s <- cases[cases$Year == i, ]$SAR %>%
    is.na() %>%
    which() %>%
    length()
  p <- cases[cases$Year == i, ]$Population %>%
    is.na() %>%
    which() %>%
    length()
  row <- data.frame(Year = i, Observed = o, Expected = e, SAR = s, Population = p)
  missing <- rbind(missing, row)

  which <- cases[cases$Year == i, ]$SAR %>%
    is.na() %>%
    which()
  missing_LA <- c(missing_LA, list(cases[cases$Year == i, ][which, ]$Level.description))
}
kable(missing, caption = "Missing Observations") %>% kable_styling(., full_width = TRUE)
Missing Observations
Year Observed Expected SAR Population
2003 121 121 121 0
2004 99 99 99 0
2005 115 115 115 0
2006 127 127 127 0
2007 109 109 109 0
2008 77 77 77 0
2009 106 106 106 0
2010 130 130 130 0
2011 116 116 116 0
2012 73 73 73 0
2013 0 0 0 0
2014 0 0 0 0
2015 0 0 0 0
2016 0 0 0 0
2017 0 0 0 0
2018 0 0 0 0
2019 0 0 0 0

Maps

Including Plots

You can also embed plots, for example:

# Create colour pallette
library(RColorBrewer)
colors <- c("#D0F0C0", brewer.pal(4, "Reds"))

cases$SAR.cut <- cut(
  cases$SAR,
  breaks = c(min(cases$SAR[!is.na(cases$SAR)]), 1, 2.5, 5, max(cases$SAR[!is.na(cases$SAR)])),
)

# Read shapefile
shapefile_path <- "../Boundaries (2019)/England/Local_Authority_Districts_December_2019_England.shp"
map <- st_read(shapefile_path)
## Reading layer `Local_Authority_Districts_December_2019_England' from data source `/Users/cam/OneDrive - Imperial College London/LRTI project/Boundaries (2019)/England/Local_Authority_Districts_December_2019_England.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 317 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -714517.9 ymin: 6422874 xmax: 196334.8 ymax: 7520900
## Projected CRS: WGS 84 / Pseudo-Mercator
# Prep data
for (i in 2013:2019) {
  data <- cases[cases$Year == i, ]
  map_i <- left_join(map, data, by = join_by(lad19cd == Level))

  p <- ggplot() +
    geom_sf(data = map_i) +
    aes(fill = SAR.cut) +
    theme_bw() +
    scale_fill_manual(values = colors) +
    guides(fill = guide_legend(title = "SAR")) +
    labs(title = paste0("SAR in ", i))
  
  print(p)
}

##Analysis

library(spdep)

map_nb <- poly2nb(map)
summary(map_nb)
## Neighbour list object:
## Number of regions: 317 
## Number of nonzero links: 1610 
## Percentage nonzero weights: 1.602165 
## Average number of links: 5.078864 
## 2 regions with no links:
## 44 50
## Link number distribution:
## 
##  0  1  2  3  4  5  6  7  8  9 10 11 
##  2  8 25 39 54 58 46 49 21 10  3  2 
## 8 least connected regions:
## 10 26 61 67 88 89 115 263 with 1 link
## 2 most connected regions:
## 51 227 with 11 links
# Convert to INLA format
nb2INLA("map.graph", map_nb)
map.adj <- paste(getwd(), "/map.graph", sep = "")

# select only complete year (2013:2019)
cases <- cases[cases$Year >= 2013, ]

# Aggregate observed and expected cases by LA (independent of year)
cases_agg <- cases %>%
  group_by(Level) %>%
  summarise(Observed = sum(Observed), Expected = sum(Expected), Population = sum(Population)) %>%
  mutate(SAR = Observed / Expected) %>%
  mutate(log.Population = log(Population)) %>%
  dplyr::rename(O = Observed, E = Expected)

kable(cases_agg) %>% kable_styling(., full_width=TRUE)
Level O E Population SAR log.Population
E06000001 637 487.3 148300 1.3072030 11.90699
E06000002 1729 882.6 238100 1.9589848 12.38045
E06000003 1212 675.9 203300 1.7931647 12.22244
E06000004 1503 1067.1 317100 1.4084903 12.66697
E06000005 991 553.3 166600 1.7910718 12.02335
E06000006 1098 707.1 209200 1.5528214 12.25105
E06000007 1253 1083.3 328300 1.1566510 12.70168
E06000008 1750 983.0 283600 1.7802645 12.55532
E06000009 1486 765.9 213600 1.9402011 12.27186
E06000010 1763 1593.2 414100 1.1065780 12.93386
E06000011 1593 1395.3 467800 1.1416900 13.05580
E06000012 801 860.3 253600 0.9310706 12.44351
E06000013 919 832.1 261700 1.1044346 12.47495
E06000014 1415 910.1 274600 1.5547742 12.52307
E06000015 1132 1542.0 436000 0.7341115 12.98540
E06000016 1795 2298.8 609700 0.7808422 13.32072
E06000017 159 167.8 58300 0.9475566 10.97336
E06000018 2244 1895.6 500000 1.1837940 13.12236
E06000019 1153 846.5 266700 1.3620791 12.49388
E06000020 1154 982.1 293800 1.1750331 12.59065
E06000021 3022 1572.4 417100 1.9219028 12.94108
E06000022 1277 835.7 263800 1.5280603 12.48295
E06000023 2829 2760.2 686100 1.0249257 13.43878
E06000024 1176 1027.3 316900 1.1447484 12.66634
E06000025 1554 1450.7 428700 1.0712070 12.96851
E06000026 2088 1387.6 387100 1.5047564 12.86644
E06000027 982 634.1 187100 1.5486516 12.13940
E06000030 1307 1313.1 363300 0.9953545 12.80298
E06000031 1725 1406.9 358000 1.2260999 12.78829
E06000032 3289 1586.9 411800 2.0725944 12.92829
E06000033 593 1017.8 286000 0.5826292 12.56375
E06000034 709 1142.8 307500 0.6204060 12.63623
E06000035 1927 1651.2 469100 1.1670300 13.05857
E06000036 516 679.8 206700 0.7590468 12.23902
E06000037 889 822.8 263800 1.0804570 12.48295
E06000038 1022 1101.7 268200 0.9276573 12.49949
E06000039 1699 1163.2 301000 1.4606259 12.61487
E06000040 1042 790.9 250700 1.3174864 12.43201
E06000041 806 867.2 281300 0.9294280 12.54718
E06000042 2728 1730.1 488500 1.5767875 13.09909
E06000043 1888 1301.8 379500 1.4502996 12.84661
E06000044 1020 1178.5 324500 0.8655070 12.69004
E06000045 1908 1446.9 369900 1.3186813 12.82099
E06000046 737 587.8 187500 1.2538278 12.14153
E06000047 3804 2471.2 746400 1.5393331 13.52302
E06000049 2888 1814.0 559900 1.5920617 13.23551
E06000050 2806 1658.2 495900 1.6921964 13.11413
E06000051 1556 1351.0 443800 1.1517395 13.00313
E06000052 4286 2553.2 791100 1.6786777 13.58118
E06000053 4286 2553.2 791100 1.6786777 13.58118
E06000054 2715 2437.9 774600 1.1136634 13.56010
E06000055 960 1005.4 285200 0.9548438 12.56095
E06000056 2075 1548.3 445800 1.3401796 13.00763
E06000057 1884 1349.4 439700 1.3961761 12.99385
E06000058 2777 1907.5 548200 1.4558322 13.21440
E06000059 2438 1461.4 507300 1.6682633 13.13686
E07000004 1700 1078.1 326500 1.5768482 12.69619
E07000005 566 441.7 162100 1.2814127 11.99597
E07000006 485 350.1 110400 1.3853185 11.61187
E07000007 1662 992.4 298600 1.6747279 12.60686
E07000008 470 632.1 172100 0.7435532 12.05583
E07000009 344 468.7 142600 0.7339450 11.86780
E07000010 696 523.5 146800 1.3295129 11.89683
E07000011 1243 919.5 271300 1.3518216 12.51098
E07000012 689 814.9 256000 0.8455025 12.45293
E07000026 632 416.8 134700 1.5163148 11.81081
E07000027 674 334.4 99200 2.0155502 11.50489
E07000028 802 539.8 156600 1.4857355 11.96145
E07000029 563 318.7 96900 1.7665516 11.48143
E07000030 254 197.2 68500 1.2880325 11.13459
E07000031 477 379.3 132300 1.2575798 11.79283
E07000032 437 572.0 177700 0.7639860 12.08785
E07000033 661 394.1 116400 1.6772393 11.66479
E07000034 978 502.6 149000 1.9458814 11.91170
E07000035 296 253.3 94700 1.1685748 11.45847
E07000036 534 577.0 170600 0.9254766 12.04708
E07000037 616 417.7 133100 1.4747426 11.79886
E07000038 721 422.7 136700 1.7057014 11.82554
E07000039 640 519.8 161000 1.2312428 11.98916
E07000040 924 557.9 185900 1.6562108 12.13296
E07000041 1064 601.3 165900 1.7694994 12.01914
E07000042 588 374.8 123700 1.5688367 11.72561
E07000043 694 436.1 136400 1.5913781 11.82335
E07000044 425 317.9 113900 1.3368984 11.64308
E07000045 819 557.7 177300 1.4685315 12.08560
E07000046 441 283.2 91200 1.5572034 11.42081
E07000047 267 217.4 73700 1.2281509 11.20776
E07000061 662 491.6 147300 1.3466233 11.90023
E07000062 662 508.3 142500 1.3023805 11.86710
E07000063 762 429.7 147200 1.7733302 11.89955
E07000064 410 348.6 119900 1.1761331 11.69441
E07000065 712 653.1 227500 1.0901853 12.33491
E07000066 849 1127.2 309800 0.7531938 12.64368
E07000067 719 775.6 240200 0.9270242 12.38923
E07000068 271 391.7 117200 0.6918560 11.67164
E07000069 254 393.7 126600 0.6451613 11.74879
E07000070 778 902.8 271400 0.8617634 12.51135
E07000071 1808 1022.5 287500 1.7682152 12.56898
E07000072 508 727.1 200300 0.6986659 12.20757
E07000073 442 582.7 151800 0.7585378 11.93032
E07000074 227 260.7 88900 0.8707326 11.39527
E07000075 222 361.7 124800 0.6137683 11.73447
E07000076 1085 626.2 196000 1.7326733 12.18587
E07000077 314 440.4 143800 0.7129882 11.87618
E07000078 775 591.3 171700 1.3106714 12.05350
E07000079 414 358.4 118900 1.1551339 11.68604
E07000080 521 365.5 121600 1.4254446 11.70849
E07000081 1346 780.2 211300 1.7251987 12.26103
E07000082 690 525.2 175800 1.3137852 12.07710
E07000083 759 464.0 134000 1.6357759 11.80560
E07000084 1632 998.0 287900 1.6352705 12.57037
E07000085 670 529.9 182800 1.2643895 12.11615
E07000086 738 681.2 205400 1.0833823 12.23271
E07000087 354 488.2 164800 0.7251127 12.01249
E07000088 373 433.5 132800 0.8604383 11.79660
E07000089 375 474.0 156800 0.7911392 11.96273
E07000090 599 593.5 184300 1.0092671 12.12432
E07000091 903 717.4 240700 1.2587120 12.39131
E07000092 493 609.7 156100 0.8085944 11.95825
E07000093 566 616.8 189300 0.9176394 12.15109
E07000094 495 558.2 190900 0.8867789 12.15950
E07000095 445 569.1 160200 0.7819364 11.98418
E07000096 880 883.4 251600 0.9961512 12.43560
E07000098 461 607.0 178100 0.7594728 12.09010
E07000099 538 728.0 212800 0.7390110 12.26811
E07000102 446 491.0 154200 0.9083503 11.94601
E07000103 578 655.4 170100 0.8819042 12.04414
E07000105 841 708.3 214700 1.1873500 12.27700
E07000106 710 655.6 221300 1.0829774 12.30727
E07000107 966 708.8 184300 1.3628668 12.12432
E07000108 672 537.8 169700 1.2495351 12.04179
E07000109 721 642.3 182600 1.1225284 12.11505
E07000110 563 927.8 267100 0.6068118 12.49538
E07000111 607 609.7 194900 0.9955716 12.18024
E07000112 627 503.5 158600 1.2452830 11.97414
E07000113 832 816.3 240600 1.0192331 12.39089
E07000114 1061 739.4 220500 1.4349473 12.30365
E07000115 475 678.0 217800 0.7005900 12.29133
E07000116 323 583.7 197600 0.5533665 12.19400
E07000117 1088 533.1 147000 2.0408929 11.89819
E07000118 718 567.6 173000 1.2649753 12.06105
E07000119 432 301.3 103300 1.4337869 11.54539
E07000120 956 492.6 137000 1.9407227 11.82774
E07000121 1146 673.8 202600 1.7008014 12.21899
E07000122 1071 566.3 155100 1.8912237 11.95183
E07000123 1446 842.6 230600 1.7161168 12.34844
E07000124 395 226.4 88100 1.7446996 11.38623
E07000125 693 369.6 113100 1.8750000 11.63603
E07000126 915 537.5 166500 1.7023256 12.02275
E07000127 728 504.4 165700 1.4432990 12.01793
E07000128 757 447.2 147700 1.6927549 11.90294
E07000129 486 492.5 152200 0.9868020 11.93295
E07000130 736 843.8 255500 0.8722446 12.45098
E07000131 388 397.4 140300 0.9763463 11.85154
E07000132 463 528.6 161600 0.8758986 11.99288
E07000133 238 234.7 75500 1.0140605 11.23189
E07000134 491 477.4 150900 1.0284876 11.92437
E07000135 260 265.2 85400 0.9803922 11.35510
E07000136 345 379.7 103700 0.9086121 11.54926
E07000137 559 568.7 178000 0.9829436 12.08954
E07000138 651 531.6 136100 1.2246050 11.82115
E07000139 613 513.2 165900 1.1944661 12.01914
E07000140 466 449.3 131500 1.0371689 11.78676
E07000141 817 665.3 213800 1.2280174 12.27280
E07000142 540 417.7 135000 1.2927939 11.81303
E07000143 835 666.9 197200 1.2520618 12.19197
E07000144 616 527.0 176400 1.1688805 12.08051
E07000145 648 501.6 146700 1.2918660 11.89615
E07000146 979 739.6 215400 1.3236885 12.28025
E07000147 447 369.9 122500 1.2084347 11.71587
E07000148 1028 755.5 193500 1.3606883 12.17303
E07000149 744 622.9 200900 1.1944132 12.21056
E07000150 574 437.8 121200 1.3111010 11.70520
E07000151 540 388.2 125800 1.3910355 11.74245
E07000152 590 447.7 146800 1.3178468 11.89683
E07000153 757 567.3 165000 1.3343910 12.01370
E07000154 2312 1455.8 383100 1.5881302 12.85605
E07000155 692 421.4 143300 1.6421452 11.87270
E07000156 691 447.1 131800 1.5455155 11.78904
E07000163 347 218.4 76700 1.5888278 11.24766
E07000164 591 373.9 124400 1.5806365 11.73126
E07000165 1176 700.5 243500 1.6788009 12.40287
E07000166 404 255.8 76600 1.5793589 11.24635
E07000167 366 212.0 72400 1.7264151 11.18996
E07000168 941 472.2 144600 1.9927997 11.88173
E07000169 706 441.0 132700 1.6009070 11.79585
E07000170 873 668.5 195100 1.3059088 12.18127
E07000171 790 565.6 170800 1.3967468 12.04825
E07000172 620 541.7 159400 1.1445450 11.97917
E07000173 636 562.6 173500 1.1304657 12.06393
E07000174 788 585.5 163400 1.3458582 12.00396
E07000175 687 567.8 178400 1.2099331 12.09178
E07000176 611 514.4 176100 1.1877916 12.07881
E07000177 976 839.0 240100 1.1632896 12.38881
E07000178 930 812.1 224900 1.1451792 12.32341
E07000179 785 714.2 220200 1.0991319 12.30229
E07000180 682 682.0 205400 1.0000000 12.23271
E07000181 550 541.2 166800 1.0162602 12.02455
E07000187 712 526.4 174100 1.3525836 12.06739
E07000188 674 585.0 182100 1.1521368 12.11231
E07000189 951 783.8 243200 1.2133197 12.40164
E07000192 599 500.8 149300 1.1960863 11.91371
E07000193 894 652.9 187200 1.3692755 12.13993
E07000194 655 453.3 147500 1.4449592 11.90158
E07000195 948 559.5 176800 1.6943700 12.08277
E07000196 519 440.6 148000 1.1779392 11.90497
E07000197 832 594.5 188400 1.3994954 12.14632
E07000198 570 388.0 134300 1.4690722 11.80783
E07000199 593 420.4 124300 1.4105614 11.73045
E07000200 545 372.8 131000 1.4619099 11.78295
E07000202 1489 855.6 224200 1.7402992 12.32029
E07000203 613 414.4 145800 1.4792471 11.88999
E07000207 928 811.3 244000 1.1438432 12.40492
E07000208 519 444.1 133400 1.1686557 11.80111
E07000209 838 697.6 217500 1.2012615 12.28995
E07000210 494 385.5 132300 1.2814527 11.79283
E07000211 880 837.6 240600 1.0506208 12.39089
E07000212 372 437.4 123400 0.8504801 11.72319
E07000213 503 577.8 154500 0.8705434 11.94795
E07000214 386 435.4 141800 0.8865411 11.86217
E07000215 416 447.4 138000 0.9298167 11.83501
E07000216 665 611.4 209000 1.0876676 12.25009
E07000217 547 612.5 172100 0.8930612 12.05583
E07000218 295 299.7 91900 0.9843177 11.42846
E07000219 597 726.4 203300 0.8218612 12.22244
E07000220 589 579.0 173100 1.0172712 12.06162
E07000221 645 532.9 174100 1.2103584 12.06739
E07000222 772 685.3 200500 1.1265139 12.20857
E07000223 400 329.6 94600 1.2135922 11.45741
E07000224 987 692.1 207400 1.4260945 12.24240
E07000225 781 501.4 161600 1.5576386 11.99288
E07000226 542 731.6 194300 0.7408420 12.17716
E07000227 601 624.6 212500 0.9622158 12.26670
E07000228 964 740.8 237800 1.3012959 12.37919
E07000229 504 531.6 158900 0.9480813 11.97603
E07000234 437 441.8 144300 0.9891354 11.87965
E07000235 315 286.9 105200 1.0979435 11.56362
E07000236 600 486.4 138400 1.2335526 11.83790
E07000237 748 559.0 155200 1.3381038 11.95247
E07000238 610 539.7 175100 1.1302576 12.07311
E07000239 611 486.4 142800 1.2561678 11.86920
E07000240 716 861.7 263800 0.8309156 12.48295
E07000241 444 633.2 183800 0.7012003 12.12160
E07000242 521 762.1 239900 0.6836373 12.38798
E07000243 456 536.4 145900 0.8501119 11.89068
E07000244 1583 1036.0 348600 1.5279923 12.76168
E07000245 1164 984.6 274300 1.1822060 12.52198
E07000246 778 692.0 214400 1.1242775 12.27560
E08000001 3131 1742.6 493600 1.7967405 13.10948
E08000002 1840 1089.6 315500 1.6886931 12.66191
E08000003 5382 3474.9 877000 1.5488215 13.68426
E08000004 3010 1513.0 431200 1.9894250 12.97433
E08000005 2672 1359.5 380700 1.9654285 12.84977
E08000006 2375 1600.8 403400 1.4836332 12.90768
E08000007 2081 1578.2 459300 1.3185908 13.03746
E08000008 2351 1322.4 364500 1.7778282 12.80628
E08000009 1760 1301.0 403000 1.3528055 12.90669
E08000010 1578 1681.7 502200 0.9383362 13.12675
E08000011 1160 889.9 242900 1.3035172 12.40041
E08000012 1448 2647.0 688200 0.5470344 13.44183
E08000013 1458 930.0 269800 1.5677419 12.50544
E08000014 1345 1302.6 397100 1.0325503 12.89194
E08000015 3033 1643.0 499300 1.8460134 13.12096
E08000016 1836 1291.6 369800 1.4214927 12.82072
E08000017 1934 1663.8 484500 1.1623993 13.09087
E08000018 1602 1422.0 418900 1.1265823 12.94539
E08000019 3301 2989.7 864300 1.1041242 13.66968
E08000021 1925 1528.8 426300 1.2591575 12.96290
E08000022 1573 1040.3 301400 1.5120638 12.61619
E08000023 997 743.9 219200 1.3402339 12.29774
E08000024 2221 1366.0 403600 1.6259151 12.90818
E08000025 9075 7674.3 2106100 1.1825183 14.56035
E08000026 2245 2082.7 566900 1.0779277 13.24794
E08000027 2379 1737.0 505800 1.3696028 13.13390
E08000028 2877 2162.2 586500 1.3305892 13.28193
E08000029 1122 1078.2 341500 1.0406233 12.74110
E08000030 1632 1735.0 488200 0.9406340 13.09848
E08000031 2092 1601.9 437200 1.3059492 12.98815
E08000032 5770 3611.5 1034600 1.5976741 13.84953
E08000033 1954 1143.6 339200 1.7086394 12.73435
E08000034 3558 2493.0 731800 1.4271961 13.50326
E08000035 3319 4594.1 1216400 0.7224484 14.01141
E08000036 2382 1872.1 522800 1.2723679 13.16695
E08000037 1426 1009.4 296400 1.4127204 12.59947
E09000001 405 1947.4 461100 0.2079696 13.04137
E09000002 1013 1751.4 447400 0.5783944 13.01121
E09000003 1363 2410.5 661400 0.5654429 13.40211
E09000004 1747 1425.3 414400 1.2257069 12.93459
E09000005 1425 2256.0 557500 0.6316489 13.23122
E09000006 1348 1918.2 536800 0.7027422 13.19338
E09000007 836 1249.8 350500 0.6689070 12.76712
E09000008 2123 2593.0 689300 0.8187428 13.44343
E09000009 1437 2331.4 595500 0.6163678 13.29716
E09000010 2437 2227.2 611900 1.0941990 13.32432
E09000011 2195 2004.0 489600 1.0953094 13.10134
E09000012 405 1947.4 461100 0.2079696 13.04137
E09000013 669 1060.2 256900 0.6310130 12.45644
E09000014 1512 1751.9 442800 0.8630630 13.00087
E09000015 1166 1603.1 422900 0.7273408 12.95489
E09000016 1015 1498.6 408800 0.6772988 12.92098
E09000017 1540 2013.8 524900 0.7647234 13.17096
E09000018 1335 1887.6 462800 0.7072473 13.04505
E09000019 923 1240.5 299800 0.7440548 12.61087
E09000020 336 784.4 207100 0.4283529 12.24096
E09000021 1398 1016.8 278500 1.3749017 12.53717
E09000022 1545 1874.6 457800 0.8241758 13.03419
E09000023 1263 2067.0 495300 0.6110305 13.11292
E09000024 1560 1414.4 339400 1.1029412 12.73493
E09000025 1055 2623.6 620300 0.4021192 13.33796
E09000026 1271 2106.1 552600 0.6034851 13.22239
E09000027 1248 1183.0 325800 1.0549451 12.69404
E09000028 1680 1948.7 465600 0.8621132 13.05108
E09000029 1451 1232.9 343200 1.1769000 12.74607
E09000030 355 2028.5 488900 0.1750062 13.09991
E09000031 1353 2033.8 480800 0.6652572 13.08321
E09000032 2044 2102.3 449500 0.9722685 13.01589
E09000033 385 1201.9 323800 0.3203262 12.68788 
# Create colour pallette
library(RColorBrewer)
colors <- c("#D0F0C0", brewer.pal(4, "Reds"))

# Create buckets for SAR
cases_agg$SAR.cut <- cut(
  cases_agg$SAR,
  breaks = c(min(cases_agg$SAR[!is.na(cases_agg$SAR)]), 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6, max(cases_agg$SAR[!is.na(cases_agg$SAR)])),
)

# Produce spatial map of aggregatged SAR
map_agg <- left_join(map, cases_agg, by = join_by(lad19cd == Level))

# Plot SAR using ggplot
library(ggplot2)
ggplot() +
  geom_sf(data = map_agg, aes(fill = SAR.cut)) +
  theme_bw() +
  scale_fill_viridis_d() +
  theme(legend.position = "bottom")

# Fit hierarchical Poisson log-linear model in INLA
library(INLA)
LA <- seq(1, 317)
formula <- O ~ f(LA, model = "bym2", graph = map.adj, hyper = list(
  prec = list(
    prior = "pc.prec",
    param = c(0.5 / 0.31, 0.01)
  ),
  phi = list(
    prior = "pc",
    param = c(0.5, 2 / 3)
  )
))

sBYM.model <- inla(formula = formula, family = "poisson", data = cases_agg, E = E, control.compute = list(waic = TRUE))

# Obtain the posterior summary statistics
summary(sBYM.model, summary.stat = c("mean", "p0.975", "p0.025", "p0.5in", "p0.5out", "waic"))
## 
## Call:
##    c("inla.core(formula = formula, family = family, contrasts = contrasts, 
##    ", " data = data, quantiles = quantiles, E = E, offset = offset, ", " 
##    scale = scale, weights = weights, Ntrials = Ntrials, strata = strata, 
##    ", " lp.scale = lp.scale, link.covariates = link.covariates, verbose = 
##    verbose, ", " lincomb = lincomb, selection = selection, control.compute 
##    = control.compute, ", " control.predictor = control.predictor, 
##    control.family = control.family, ", " control.inla = control.inla, 
##    control.fixed = control.fixed, ", " control.mode = control.mode, 
##    control.expert = control.expert, ", " control.hazard = control.hazard, 
##    control.lincomb = control.lincomb, ", " control.update = 
##    control.update, control.lp.scale = control.lp.scale, ", " 
##    control.pardiso = control.pardiso, only.hyperparam = only.hyperparam, 
##    ", " inla.call = inla.call, inla.arg = inla.arg, num.threads = 
##    num.threads, ", " blas.num.threads = blas.num.threads, keep = keep, 
##    working.directory = working.directory, ", " silent = silent, inla.mode 
##    = inla.mode, safe = FALSE, debug = debug, ", " .parent.frame = 
##    .parent.frame)") 
## Time used:
##     Pre = 8.37, Running = 0.199, Post = 0.0186, Total = 8.59 
## Fixed effects:
##              mean    sd 0.025quant 0.5quant 0.975quant  mode kld
## (Intercept) 0.131 0.004      0.122    0.131       0.14 0.131   0
## 
## Random effects:
##   Name     Model
##     LA BYM2 model
## 
## Model hyperparameters:
##                   mean    sd 0.025quant 0.5quant 0.975quant  mode
## Precision for LA 9.631 0.914      7.877     9.62     11.472 9.646
## Phi for LA       0.962 0.031      0.882     0.97      0.996 0.988
## 
## Watanabe-Akaike information criterion (WAIC) ...: 3272.77
## Effective number of parameters .................: 158.24
## 
## Marginal log-Likelihood:  -1978.99 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')
RR_sBYM <- c()

for (i in 1:317) {
  RR_sBYM[i] <- inla.emarginal(
    function(x) exp(x),
    sBYM.model$marginals.random$LA[[i]]
  )
}

# Posterior probabilities
RR_sBYM_marg <- sBYM.model$marginals.random$LA[1:317]
PP_sBYM <- lapply(RR_sBYM_marg, function(x) {
  1 - inla.pmarginal(0, x)
})


# Obtain the posterior estimates from the spatial model to be plotted, that is (i) the area level posterior mean of the residual RRs and (ii) the posterior probability (PP) that the residual RRs > 1.
resRR_PP <- data.frame(
  resRR = RR_sBYM,
  PP = unlist(PP_sBYM),
  Level = cases_agg[, 1]
)

# Using ggplot2 package, produce a map of the posterior mean of the residual RRs and the posterior probabilities that the residual RRs are > 1
# For the map of the posterior mean of the residual RRs, use the following breakpoints [min,0.4], (0.4-0.6], (0.6-0.8], (0.8,1], (1,1.2], (1.2-1.4], (1.4-1.6], (1.6-max].
resRR_PP$resRRcat <- cut(resRR_PP$resRR, breaks = c(
  min(resRR_PP$resRR),
  0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6,
  max(resRR_PP$resRR)
), include.lowest = T)

# breakpoints
resRR_PP$PPcat <- cut(resRR_PP$PP, c(0, 0.2, 0.8, 1.00), include.lowest = TRUE)
map_RR_PP <- left_join(map, resRR_PP, by = c("lad19cd" = "Level"))

# Produce the maps of the posterior mean of the residual RRs and the posterior probabilities PP using ggplot2 package.
library(viridis)
p1 <- ggplot() +
  geom_sf(data = map_RR_PP) +
  aes(fill = resRRcat) +
  theme_bw() +
  scale_fill_brewer(palette = "PuOr") +
  guides(fill = guide_legend(title = "RR")) +
  ggtitle("RR Spatial model") +
  theme(
    text = element_text(size = 15),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
  )

p2 <- ggplot() +
  geom_sf(data = map_RR_PP) +
  aes(fill = PPcat) +
  theme_bw() +
  scale_fill_viridis(
    option = "plasma", name = "PP",
    discrete = T,
    direction = -1,
    guide = guide_legend(
      title.position = "top",
      reverse = T
    )
  ) +
  ggtitle("PP Spatial model") +
  theme(
    text = element_text(size = 15),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(), plot.title = element_text(size = 12, face = "bold")
  )

print(p1)

print(p2)

sBYM.model$summary.hyperpar
##                       mean         sd 0.025quant  0.5quant 0.975quant      mode
## Precision for LA 9.6307030 0.91409535  7.8773681 9.6182571 11.4716396 9.6461504
## Phi for LA       0.9616735 0.03059103  0.8818594 0.9697759  0.9958556 0.9876124