d<-read.csv("https://www.mortality.org/Public/STMF/Outputs/stmf.csv",skip=1)
save(d, file="data.rob")
load("data.rob")
filter(d, Year >1999) -> d
d$date<-paste(d$Year,d$Week,1,sep="-")
d$date<-as.Date(d$date, "%Y-%U-%u")

filter(d, Sex=="b") -> both
both[,c(1,20,11:16)] -> death_rate
pivot_longer(death_rate,cols=3:8) -> death_rate
death_rate$name<-gsub("R","",death_rate$name)


both[,c(1,20,5:10)] -> deaths
pivot_longer(deaths,cols=3:8) -> deaths
deaths$name<-gsub("D","",deaths$name)

Fit a gamm to total

Smooth on week for all years. Allow some different form by modelling week conditioned on country code as the random effect.

death_rate %>% filter(name == "Total") -> total
total$year<-year(total$date)
total$week<-week(total$date)

gamm.mod<-gamm4(value~s(week,k=20),data=total,random = ~ (week|CountryCode))

Extract the fixed, random and residuals

total$fixed<-predict(gamm.mod$gam)
total$rand<-predict(gamm.mod$mer)
total$residfixed<-residuals(gamm.mod$gam)
total$residrand<-residuals(gamm.mod$mer)

Plot out with a green line for the random effect and red for the fixed

ct <- levels(total$CountryCode)[1:18]
total %>% filter (CountryCode %in% ct) -> ttl
g0<-ggplot(ttl,aes(x=week,y=value))
g1<-g0+geom_point(alpha=0.5,cex=0.5)
g1<-g1+geom_line(aes(x=week,y=fixed),colour=2,lwd=0.5) 
g1<-g1+geom_line(aes(x=week,y=rand),colour=3,lwd=0.5) 
#g1<-g1+labs(y = ylabel,x=xlabel,title="Fixed and random splines")
g1+facet_wrap(~CountryCode,ncol = 6) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

65 to 74

death_rate %>% filter(name == "65_74") -> total
total$year<-year(total$date)
total$week<-week(total$date)

gamm.mod<-gamm4(value~s(week,k=20),data=total,random = ~ (week|CountryCode))

Extract the fixed, random and residuals

total$fixed<-predict(gamm.mod$gam)
total$rand<-predict(gamm.mod$mer)
total$residfixed<-residuals(gamm.mod$gam)
total$residrand<-residuals(gamm.mod$mer)

Plot out with a green line for the random effect and red for the fixed

ct <- levels(total$CountryCode)[1:18]
total %>% filter (CountryCode %in% ct) -> ttl
g0<-ggplot(ttl,aes(x=week,y=value))
g1<-g0+geom_point(alpha=0.5,cex=0.5)
g1<-g1+geom_line(aes(x=week,y=fixed),colour=2,lwd=0.5) 
g1<-g1+geom_line(aes(x=week,y=rand),colour=3,lwd=0.5) 
#g1<-g1+labs(y = ylabel,x=xlabel,title="Fixed and random splines")
g1+facet_wrap(~CountryCode,ncol = 6) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Where to go next

The idea is that excess deaths can be regarded from two perspective. The residuals on the country level random effect represents excess over the expected value for that country. However the excess can be measured from the pooled fixed effect. Some countries always have excesses on that measure, others tend tolie below the line

library(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(dygraphs)
cntrs<-c("GBRTENW","ESP","USA","SWE")
cntrs<-c("ESP")
total %>% filter(total$CountryCode %in% cntrs)  -> cntrs

dt<-xts(cntrs$residrand ,order.by = cntrs$date)
dygraph(dt)
library(mgcv)
max(cntrs$date)
## [1] "2020-09-14"
#lubridate::month(cntrs$date)
acf(residuals(lm(cntrs$residfixed~cntrs$date)))

#cntrs$date
# pivot_wider(cntrs,,names_from="CountryCode",values_from="residfixed") -> cntrs
# 
# dd<-data.frame(date=GB$date,fixed=GB$residfixed,rand=GB$residrand)
# 
# library(xts)
# library(dygraphs)
# dt<-xts(dd[2:3],order.by = dd$date)
# dygraph(dt)

https://robjhyndman.com/hyndsight/excess-deaths/

# dd %>% pivot_longer(2:3) %>% ggplot(aes(x=date,y=value,col=name)) +geom_line()