#devtools::install_github("dgolicher/usdata") Install data package by uncommenting this line
library(tidyverse)
library(tmap)
library(mapview)
library(sjPlot)
library(usdata)
data(ukdata)
library(giscourse)
library(spdep)
#dir(pattern="*shp")
#dd<-st_read("209e5dd2-cc00-4f2d-a041-a35b0e0e62482020330-1-9o3ovk.cyfal.shp")
#library(readr)
#lookup <- read_csv("Output_Area_to_Lower_Layer_Super_Output_Area_to_Middle_Layer_Super_Output_Area_to_Local_Authority_District__December_2020__Lookup_in_England_and_Wales.csv")[,c(3,5)]
#names(lookup)<-tolower(names(lookup))
#dd %>% left_join(lookup) ->ddd
#str(dd)
#ddd %>% group_by(msoa11cd) %>% summarise(geom=st_union(geometry))-> mm
#st_write(mm,"msoasshape.shp")
con<-sconnect11()

#st_write(mm,con, "msoasshape")
mm<-st_read(con,"msoasshape")
msoa_wave<-merge(msoa_wave,msoas_area)
msoa_wave %>% mutate(density=pop/(st_areasha/10000)) ->msoa_wave
msoa_wave$Region<-as.factor(msoa_wave$RegionNation)

Method and explanation

To be added. Response is the number of deaths divided by the number of people over 60 in each MSOA. House prices and population density are centred and scaled to each region. Mean neighbours is the response averaged over the neighbouring MSOAS, centred and scaled to arnge between 1 and -1 within each of the regions. This helps to enable direct comparisons of the relative effect size of each of the variables though looking at the fitted model coefficients plotted with confidence intervals.

Wave

First

msoa_wave %>% filter(wave=="wave1") ->w

Covid

w %>% mutate(response =1000*(covid/over60)) -> w

Model

scale.center <- function(x){
  (x-mean(x,na.rm=TRUE)) /(max(x,na.rm = TRUE)-min(x,na.rm = TRUE))
}
w %>% group_by(Region) %>% mutate(price=-scale.center(log(price))) %>%  mutate(density=scale.center(density)) -> ww
ww<-na.omit(ww)
mod<-lm(data=ww, response ~Region+ price +density)

plot_model(mod, sort.est = TRUE)

ww$residuals<-residuals(mod)

Moran’s

mm %>% right_join(ww) ->ww_geo
s<-as(ww_geo,"Spatial")
nb <- poly2nb(s, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
moran.test(s$response,lw, alternative="greater", zero.policy=TRUE)
## 
##  Moran I test under randomisation
## 
## data:  s$response  
## weights: lw  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 37.518, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.805510e-01     -1.473188e-04      5.597449e-05

Neighbour model

ww$mean_neighbours <- lag.listw(lw, s$response)
ww %>% group_by(Region) %>% mutate(mean_neighbours=scale.center(mean_neighbours))  -> ww
mod<-lm(data=ww, response ~Region+ price  + mean_neighbours + density )
plot_model(mod, sort.est = TRUE)

Neighbour map

msoahex %>% right_join(ww, by="msoa11cd") ->w_geo
library(RColorBrewer)
colors <- brewer.pal(7, "RdYlGn")[7:1]
qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("mean_neighbours", palette = colors, midpoint = NA, style="quantile") + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

London

 mapviewOptions(basemaps =c("CartoDB.Positron",
                            "OpenStreetMap",
                            "Esri.WorldImagery",
                             "Esri.WorldStreetMap",
                             "Esri.WorldShadedRelief",
                             "Stamen.Terrain"))

ww_geo[,c("Region","Laname","price","covid", "All","density","response","pop", "over60","geom")] %>% filter(Region=="London") %>%
mapview(zcol="response", col.regions = colors) ->mp
mp@map %>% addFullscreenControl()

All cause

w %>% mutate(response =1000*(All/over60)) -> w

Model

scale.center <- function(x){
  (x-mean(x,na.rm=TRUE)) /(max(x,na.rm = TRUE)-min(x,na.rm = TRUE))
}
w %>% group_by(Region) %>% mutate(price=-scale.center(log(price))) %>%  mutate(density=scale.center(density)) -> ww
ww<-na.omit(ww)

#mapview(ww_geo,zcol="residuals")
mod<-lm(data=ww, response ~Region+ price +density )
plot_model(mod, sort.est = TRUE)

ww$residuals<-residuals(mod)

Moran’s

mm %>% right_join(ww) ->ww_geo
s<-as(ww_geo,"Spatial")
nb <- poly2nb(s, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
moran.test(s$response,lw, alternative="greater", zero.policy=TRUE)
## 
##  Moran I test under randomisation
## 
## data:  s$response  
## weights: lw  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 28.248, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.113124e-01     -1.473188e-04      5.603794e-05

Neighbour model

ww$mean_neighbours <- lag.listw(lw, s$response)
ww %>% group_by(Region) %>% mutate(mean_neighbours=scale.center(mean_neighbours))  -> ww
mod<-lm(data=ww, response ~Region+ price  + mean_neighbours + density )
plot_model(mod, sort.est = TRUE)

Neighbour map

msoahex %>% right_join(ww, by="msoa11cd") ->w_geo
library(RColorBrewer)
colors <- brewer.pal(7, "RdYlGn")[7:1]
qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("mean_neighbours", palette = colors, midpoint = NA, style="quantile") + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

London

mapviewOptions(basemaps =c("CartoDB.Positron",
                            "OpenStreetMap",
                            "Esri.WorldImagery",
                             "Esri.WorldStreetMap",
                             "Esri.WorldShadedRelief",
                             "Stamen.Terrain"))

ww_geo[,c("Region","Laname","price","covid", "All","density","response","pop", "over60","geom")] %>% filter(Region=="London") %>%
mapview(zcol="response", col.regions = colors) ->mp
mp@map %>% addFullscreenControl()

Second

msoa_wave %>% filter(wave=="wave2") ->w

Covid

w %>% mutate(response =1000*(covid/over60)) -> w

Model

scale.center <- function(x){
  (x-mean(x,na.rm=TRUE)) /(max(x,na.rm = TRUE)-min(x,na.rm = TRUE))
}
w %>% group_by(Region) %>% mutate(price=-scale.center(log(price))) %>%  mutate(density=scale.center(density)) -> ww
ww<-na.omit(ww)

#mapview(ww_geo,zcol="residuals")
mod<-lm(data=ww, response ~Region+ price +density )
plot_model(mod, sort.est = TRUE)

ww$residuals<-residuals(mod)

Moran’s

mm %>% right_join(ww) ->ww_geo
s<-as(ww_geo,"Spatial")
nb <- poly2nb(s, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
moran.test(s$response,lw, alternative="greater", zero.policy=TRUE)
## 
##  Moran I test under randomisation
## 
## data:  s$response  
## weights: lw  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 54.157, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      4.050812e-01     -1.473188e-04      5.598707e-05

Neighbour model

ww$mean_neighbours <- lag.listw(lw, s$response)
ww %>% group_by(Region) %>% mutate(mean_neighbours=scale.center(mean_neighbours))  -> ww
mod<-lm(data=ww, response ~Region+ price  + mean_neighbours + density)
plot_model(mod, sort.est = TRUE)

Neighbour map

msoahex %>% right_join(ww, by="msoa11cd") ->w_geo
library(RColorBrewer)
colors <- brewer.pal(7, "RdYlGn")[7:1]
qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("mean_neighbours", palette = colors, midpoint = NA, style="quantile") + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

North west

mapviewOptions(basemaps =c("CartoDB.Positron",
                            "OpenStreetMap",
                            "Esri.WorldImagery",
                             "Esri.WorldStreetMap",
                             "Esri.WorldShadedRelief",
                             "Stamen.Terrain"))
ww_geo[,c("Region","Laname","price","covid", "All","density","response","pop", "over60","geom")] %>% filter(Region=="North West") %>%
mapview(zcol="response", col.regions = colors) ->mp
mp@map %>% addFullscreenControl()

All cause

w %>% mutate(response =1000*(All/over60)) -> w

Model

scale.center <- function(x){
  (x-mean(x,na.rm=TRUE)) /(max(x,na.rm = TRUE)-min(x,na.rm = TRUE))
}
w %>% group_by(Region) %>% mutate(price=-scale.center(log(price))) %>%  mutate(density=scale.center(density)) -> ww
ww<-na.omit(ww)

#mapview(ww_geo,zcol="residuals")
mod<-lm(data=ww, response ~Region+ price +density )
plot_model(mod, sort.est = TRUE)

ww$residuals<-residuals(mod)

Moran’s

mm %>% right_join(ww) ->ww_geo
s<-as(ww_geo,"Spatial")
nb <- poly2nb(s, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
moran.test(s$response,lw, alternative="greater", zero.policy=TRUE)
## 
##  Moran I test under randomisation
## 
## data:  s$response  
## weights: lw  n reduced by no-neighbour observations
##   
## 
## Moran I statistic standard deviate = 41.337, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.093342e-01     -1.473188e-04      5.605265e-05

Neighbour model

ww$mean_neighbours <- lag.listw(lw, s$response)
ww %>% group_by(Region) %>% mutate(mean_neighbours=scale.center(mean_neighbours))  -> ww
mod<-lm(data=ww, response ~Region+ price  + mean_neighbours + density )
plot_model(mod, sort.est = TRUE)

Neighbour map

msoahex %>% right_join(ww, by="msoa11cd") ->w_geo
library(RColorBrewer)
colors <- brewer.pal(7, "RdYlGn")[7:1]
qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("mean_neighbours", palette = colors, midpoint = NA, style="quantile") + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

North west

mapviewOptions(basemaps =c("CartoDB.Positron",
                            "OpenStreetMap",
                            "Esri.WorldImagery",
                             "Esri.WorldStreetMap",
                             "Esri.WorldShadedRelief",
                             "Stamen.Terrain"))
ww_geo[,c("Region","Laname","price","covid", "All","density","response","pop", "over60","geom")] %>% filter(Region=="North West") %>%
mapview(zcol="response", col.regions = colors) ->mp
mp@map %>% addFullscreenControl()