# this file contains the R code used in

# Data presented by the UK government as lockdown was eased shows the transmission of COVID-19 had already increased.

# by Mike Lonergan <mlonergan@dundee.ac.uk>

# the analysis was carried out in R 4.0.0


library(readxl)

cases<-as.data.frame(read_excel("C:/Mike/covid/regional/2020-06-23_COVID-19_Press_Conference_Data.xlsx",5,skip=4))[1:95,1:2]
cases$Date<-as.numeric(cases$Date)

# 1 Jan is DATE 43810
cases$d2020<-cases$Date-43830

# day of week for cyclic smooth
cases$dow<-(cases$d2020) %%7
# it is easier if one 0 is coded 8 to force the cycle to be weekly (otherwise days 0 and 6 coincide)
cases$dow[cases$dow==0][6]<-8


library(mgcv)

g1<-gam(Cases~s(d2020,k=20)+s(dow,bs='cc',k=7),data=cases,family=nb())
summary(g1)
plot(g1,residuals=TRUE,cex=4,ylim=c(-2,2))


g2<-gam(Cases~s(d2020,k=20)+s(dow,bs='cc',k=7),data=cases[1:(nrow(cases)-14),],family=nb())
summary(g2)

plot(g2,residuals=TRUE,cex=4,ylim=c(-2,2))


pg0<-predict(g1,newdata=cases,se=TRUE)

# show everything without day of week cycle
cases2<-cases
cases2$dow<-1 # Wednesday
pg1<-predict(g1,newdata=cases2,se=TRUE)
pg2<-predict(g2,newdata=cases2,se=TRUE)

#tiff("C:/mike/covid/regional/ukjunegamfig2.tif",width=720)

plot(cases$d2020,cases$Cases,pch=20,xlab="days since 1 Jan 2020",ylab="number of newly diagnosed cases")

lines(cases$d2020,g1$fit,col="grey",lwd=2)
lines(cases$d2020,exp(pg0$fit+2*pg0$se.fit),lty="dashed",lwd=2,col="grey")
lines(cases$d2020,exp(pg0$fit-2*pg0$se.fit),lty="dashed",lwd=2,col="grey")

lines(cases$d2020,exp(pg1$fit),lwd=3,col=2)
lines(cases$d2020,exp(pg1$fit+2*pg1$se.fit),lty="dashed",lwd=3,col=2)
lines(cases$d2020,exp(pg1$fit-2*pg1$se.fit),lty="dashed",lwd=3,col=2)

lines(cases$d2020,exp(pg2$fit),col=1,lwd=3)
lines(cases$d2020,exp(pg2$fit+2*pg2$se.fit),lty="dashed",col=1,lwd=3)
lines(cases$d2020,exp(pg2$fit-2*pg2$se.fit),lty="dashed",col=1,lwd=3)

#dev.off()


# estimate 4 exponential rates
# first up to day 92, second to 129 (uninteresting), third to 160, third from then

g3<-gam(Cases~I(pmin(d2020,92.5))+I(pmin(129.5,pmax(d2020,92.5))) + I(pmin(160.5,pmax(d2020,129.5))) + 
                         I(pmax(d2020,160.5))  + s(dow,bs='cc',k=7), data=cases,family=nb())

pg30<-predict(g3,newdata=cases,se=TRUE)
pg3<-predict(g3,newdata=cases2,se=TRUE)

#tiff("C:/mike/covid/regional/ukjunegamfig3.tif",width=720)

plot(cases$d2020,cases$Cases,pch=20,xlab="days since 1 Jan 2020",ylab="number of newly diagnosed cases")

lines(cases$d2020,exp(pg30$fit),col="grey",lwd=1)
lines(cases$d2020,exp(pg30$fit+2*pg30$se.fit),lty="dashed",lwd=1,col="grey")
lines(cases$d2020,exp(pg30$fit-2*pg30$se.fit),lty="dashed",lwd=1,col="grey")

lines(cases$d2020,exp(pg3$fit),lwd=2)
lines(cases$d2020,exp(pg3$fit+2*pg3$se.fit),lty="dashed",lwd=2)
lines(cases$d2020,exp(pg3$fit-2*pg3$se.fit),lty="dashed",lwd=2)

abline(v=c(92.5,129.5,160.5),lty="dotted")

summary(g3)

# calculate R0
library(R0)

# Nishiura et al 2020: lognormal mean 4.7 sd 2.9
mGT <- generation.time("lognormal", c(4.7,2.9))

# because these are only a sample of the people, R0 doesn't like working directly off them
# so go via the slope estimates


text(85,6110,paste("R=",round(R0:::R.from.r(summary(g3)$p.table[2,1],GT=mGT),2),sep=""))
text(85,5900,paste("(95%CI:", round(R0:::R.from.r(summary(g3)$p.table[2,1]-2*summary(g3)$p.table[2,2],GT=mGT),2),
                      ",",round(R0:::R.from.r(summary(g3)$p.table[2,1]+2*summary(g3)$p.tabl[2,2],GT=mGT),2),")",
                          sep=""))


text(145,6110,paste("R=",round(R0:::R.from.r(summary(g3)$p.table[4,1],GT=mGT),2),sep=""))
text(145,5900,paste("(95%CI:", round(R0:::R.from.r(summary(g3)$p.table[4,1]-2*summary(g3)$p.table[4,2],GT=mGT),2),
                      ",",round(R0:::R.from.r(summary(g3)$p.table[4,1]+2*summary(g3)$p.table[4,2],GT=mGT),2),
                            ")",sep=""))


text(170,6110,paste("R=",round(R0:::R.from.r(summary(g3)$p.table[5,1],GT=mGT),2),sep=""))
text(170,5900,paste("(95%CI:", round(R0:::R.from.r(summary(g3)$p.table[5,1]-2*summary(g3)$p.table[5,2],GT=mGT),2),
                      ",",round(R0:::R.from.r(summary(g3)$p.table[5,1]+2*summary(g3)$p.table[5,2],GT=mGT),2),
                             ")",sep=""))


#dev.off()

# definition of periods
cases[cases$d2020 %in% c(81,92,129,160,175),]


# each period
# slopes
for(i in c(2:5))
print(c(i,summary(g3)$p.table[i,1]+c(0,-2,2)*summary(g3)$p.table[i,2]))

# R
for(i in c(2:5))
   print(paste(i,R0:::R.from.r(summary(g3)$p.table[i,1],GT=mGT),
                    R0:::R.from.r(summary(g3)$p.table[i,1]-2*summary(g3)$p.table[i,2],GT=mGT),
                       R0:::R.from.r(summary(g3)$p.table[i,1]+2*summary(g3)$p.table[i,2],GT=mGT)))