#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)
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
ww_geo[,c("Region","Laname","price","covid", "All","density","response","residuals","geom")] %>% filter(Region=="London") %>%
mapview(zcol="response", col.regions = colors)
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)
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
ww_geo[,c("Region","Laname","price","covid", "All","density","response","residuals","geom")] %>% filter(Region=="North West") %>%
mapview(zcol="response", col.regions = colors)
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)
