First wave vs second wave

#devtools::install_github("dgolicher/usdata") Install data package by uncommenting this line
library(tidyverse)
library(tmap)
library(mapview)
library(sjPlot)
library(usdata)
data(ukdata)

First wave (up to 2020-8-15) deaths

All cause deaths per thousand people over 60

Total over all msoas

House prices

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)
msoa_wave %>% filter(wave=="wave1") ->w
w %>% mutate(response =1000*(All/over60)) -> w
g0<-ggplot(w,aes(y=response, x=log10(price))) + ylab("All cause deaths in the first wave \n Per thousand people over 60") + xlab("Log10 house prices")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("All cause deaths in the first wave \n Per thousand people over 60") + xlab("Population density")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm")

By region

House prices

g0<-ggplot(w,aes(y=response, x=log10(price)), legend=FALSE) + ylab("All cause deaths in the first wave \n Per thousand people over 60") + xlab("Log10 house prices") +facet_wrap(~Region)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("All cause deaths in the first wave \n Per thousand people over 60") + xlab("Population density") +facet_wrap(~Region)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Statistical models

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

Region

mod<-lm(data=ww, response~Region)
plot_model(mod, sort.est = TRUE)

Region + price

mod<-lm(data=ww, response ~Region+ price)
plot_model(mod, sort.est = TRUE)

Region + price + density

mod<-lm(data=ww, response ~Region+ price +density)
plot_model(mod, sort.est = TRUE)

Region and price interactions

mod<-lm(data=ww, response ~Region* price)
plot_model(mod)

Explanation of model fitting

The first model shows differences between regions, ignoring any trend associated with the within region house price for the MSOA. This is equivalent to plotting out the mean deaths per region over the msoas with confidence intervals.

The second model uses logged house prices. However these are centred and scaled to lie between 0 and 1 for each of the regions. An additive model is then built using this scaled relative house prices as input. The interpretation of the coefficient for price is that the effect (slope of the regression line) is uniform across the country after scaling the variability within each region. An unscaled model would be unsuitable as it would pick up on differences in house prices between regions. This approach avoids that artefact and separates regional differences from differences in their mean house prices.

The third model includes an interaction term. This allows any differences for within region trends associated with house prices to be detected. So, if the slope of the line is steeper of flatter in some regions this will show up as an interaction. This is the most complete moel and will tend to have the best statistical properties when diagnostics are investigated. However it the hardest to interpret directly.

Mapped

All cause mortality

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

Centred and scaled house prices

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

qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("price", palette = colors[7:1]) + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

Centred and scaled population density

qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("density", palette = colors) + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

Regions

qtm(backhex) +
tm_shape(w_geo) +
  tm_fill("Region") + tm_borders(lwd = 0) + tm_legend(legend.outside =TRUE)

Covid deaths deaths per thousand of the population over 60

Total over all msoas

House prices

msoa_wave %>% filter(wave=="wave1") ->w
w %>% mutate(response =1000*(covid/over60)) -> w

g0<-ggplot(w,aes(y=response, x=log10(price))) + ylab("Covid deaths in the first wave \n Per thousand people over 60") + xlab("Log10 house prices")
g0 + geom_point(size=0.1, aes(col=Region)) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("All cause deaths in the first wave \n Per thousand people over 60") + xlab("Population density")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm")

By region

House prices

g0<-ggplot(w,aes(y=response, x=log10(price))) + ylab("Covid deaths in the first wave \n Per thousand people over 60") + xlab("Log10 house prices") +facet_wrap(~RegionNation)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("All cause deaths in the first wave \n Per thousand people over 60") + xlab("Population density") +facet_wrap(~Region)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Statistical models

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

Region

mod<-lm(data=ww, response~Region)
plot_model(mod, sort.est = TRUE)

Region + price

mod<-lm(data=ww, response ~Region+ price)
plot_model(mod, sort.est = TRUE)

Region + price + density

mod<-lm(data=ww, response ~Region+ price +density)
plot_model(mod, sort.est = TRUE)

Region and price interactions

mod<-lm(data=ww, response ~Region* price)
plot_model(mod)

Mapped

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

Second wave (after 2020-8-15) deaths

All cause deaths per thousand people over 60

Total over all msoas

House prices

msoa_wave$Region<-as.factor(msoa_wave$RegionNation)
msoa_wave %>% filter(wave=="wave2") ->w
w %>% mutate(response =1000*(All/over60)) -> w
g0<-ggplot(w,aes(y=response, x=log10(price))) + ylab("All cause deaths in the second wave \n Per thousand people over 60") + xlab("Log10 house prices")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("All cause deaths in the second wave \n Per thousand people over 60") + xlab("Population density")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm")

By region

House prices

g0<-ggplot(w,aes(y=response, x=log10(price)), legend=FALSE) + ylab("All cause deaths in the second wave \n Per thousand people over 60") + xlab("Log10 house prices") +facet_wrap(~Region)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("All cause deaths in the second wave \n Per thousand people over 60") + xlab("Population density") +facet_wrap(~Region)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Statistical models

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

Region

mod<-lm(data=ww, response~Region)
plot_model(mod, sort.est = TRUE)

Region + price

mod<-lm(data=ww, response ~Region+ price)
plot_model(mod, sort.est = TRUE)

Region + price + density

mod<-lm(data=ww, response ~Region+ price +density)
plot_model(mod, sort.est = TRUE)

Region and price interactions

mod<-lm(data=ww, response ~Region* price)
plot_model(mod)

Mapped

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

Covid deaths deaths per thousand of the population over 60

Total over all msoas

House prices

msoa_wave %>% filter(wave=="wave2") ->w
w %>% mutate(response =1000*(covid/over60)) -> w

g0<-ggplot(w,aes(y=response, x=log10(price))) + ylab("Covid deaths in the second wave \n Per thousand people over 60") + xlab("Log10 house prices")
g0 + geom_point(size=0.1, aes(col=Region)) + geom_smooth(se=FALSE)

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("Covid deaths in the second wave \n Per thousand people over 60") + xlab("Population density")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm")

By region

House prices

g0<-ggplot(w,aes(y=response, x=log10(price))) + ylab("Covid deaths in the second wave \n Per thousand people over 60") + xlab("Log10 house prices") +facet_wrap(~RegionNation)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Population density

g0<-ggplot(w,aes(y=response, x=density)) + ylab("Covid deaths in the second wave \n Per thousand people over 60") + xlab("Population density") +facet_wrap(~Region)
g0 + geom_point(size=0.1) + geom_smooth(method="lm")

Statistical models

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

Region

mod<-lm(data=ww, response~Region)
plot_model(mod, sort.est = TRUE)

Region + price

mod<-lm(data=ww, response ~Region+ price)
plot_model(mod, sort.est = TRUE)

Region + price + density

mod<-lm(data=ww, response ~Region+ price + density)
plot_model(mod, sort.est = TRUE)

Region and price interactions

mod<-lm(data=ww, response ~Region* price)
plot_model(mod)

Mapped

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

First wave against second wave

All cause

msoa_wave %>% mutate(response =1000*(All/over60)) %>% filter(response<100)-> ms
ms %>% pivot_wider(id_cols=c("Region","msoa11cd", "price"),names_from = wave,values_from = response) -> w

g0<-ggplot(w,aes(y=wave1, x=wave2)) + ylab("All cause second wave deaths \n Per thousand people over 60") + xlab("All cause first wave deaths \n Per thousand people over 60")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm") +facet_wrap(~Region)

Covid

msoa_wave %>% mutate(response =1000*(covid/over60)) %>% filter(response<12)-> ms
ms %>% pivot_wider(id_cols=c("Region","msoa11cd", "price"),names_from = wave,values_from = response) -> w

g0<-ggplot(w,aes(x=wave1, y=wave2)) + ylab("Covid second wave deaths \n Per thousand people over 60") + xlab("Covid first wave deaths \n Per thousand people over 60")
g0 + geom_point(size=0.1, aes(col= Region)) + geom_smooth(method="lm") +facet_wrap(~Region)

Model for covid deaths in wave 2

House prices and wave 1 deaths both centred and scaled to region.

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(wave1=scale.center(wave1)) -> ww

mod<-lm(data=ww, wave2 ~ wave1+Region+ price)
plot_model(mod)