hi
I expect people have seen stuff like this, but I've learned how to run the
first stage of Alison Hill's SEIR model, which makes it possible to set
parameters to fit what we know re UK growth rates in the exponential
period.
attached is a chart and the R code for it.
we knew (last week) that the log growth rate is around 0.26. I've set the
transmission rate for infectious (mild) cases to try to achieve that.
I'd guess that whatever mitigation is now in place will not have affected
the deaths yet - as it is over 20 days from infection to death.
look at the 3rd panel, which is a close-up of the deaths. We know there
were 281 deaths as of 22 March. Therefore, from the graph, we are on day
52, i.e. the UK story began 52 days earlier, on 31 Jan - exactly when they
first found it.
So I've added grey dotted vertical lines at day 52 in the other panels.
This shows you where we are, and what's just about to happen if the
mitigation doesn't work.
Unlike my previous comments predicting from exponential growth, these are
predictions from the SEIR model. But I've only modified the transmission
rate (mild) - leaving all other parameters at the defaults from alhill
site.
I'm not interested in interactive display (shiny etc) at this stage (maybe
never), just the model, model fitting, and parameter values in relation to
published data incl clinical. And some day, the impact of mitigation.
please email me off list if you want to talk about the code.
stay well
Greg
******************************************************
Please note that if you press the 'Reply' button your
message will go only to the sender of this message.
If you want to reply to the whole list, use your mailer's
'Reply-to-All' button to send your message automatically
to [log in to unmask]
Disclaimer: The messages sent to this list are the views of the sender and cannot be assumed to be representative of the range of views held by subscribers to the Radical Statistics Group. To find out more about Radical Statistics and its aims and activities and read current and past issues of our newsletter you are invited to visit our web site www.radstats.org.uk.
*******************************************************
#setwd("Desktop/ incid/COVID")
library(deSolve)
N=6.6*10^7
a=0.2
b1=0.98/N
b2=0/N
b3=0/N
g1=0.133
g2=0.188
g3=0.06
p1=0.033
p2=0.062
u=0.04
parm<-c(N,a,b1,b2,b3,g1,g2,g3,p1,p2,u)
fun<-function(times,y,parm)
{
N=parm[1]
a=parm[2]
b1=parm[3]
b2=parm[4]
b3=parm[5]
g1=parm[6]
g2=parm[7]
g3=parm[8]
p1=parm[9]
p2=parm[10]
u=parm[11]
S<-y[1]
E<-y[2]
I1<-y[3]
I2<-y[4]
I3<-y[5]
R<-y[6]
D<-y[7]
N<-y[8]
dS<--b1*I1*S-b2*I2*S-b3*I3*S
dE<--dS-a*E
dI1<-a*E-(p1+g1)*I1
dI2<-p1*I1-(p2+g2)*I2
dI3<-p2*I2-(u+g3)*I3
dR<-g1*I1+g2*I2+g3*I3
dD<-u*I3
dN<-0
#ignoring births, deaths, migration
list(c(dS,dE,dI1,dI2,dI3,dR,dD,dN))
}
times<-seq(0,300,by=1)
yinit<-c(N-1,1,0,0,0,0,0,N)
names(yinit)<-c("S","E","I1","I2","I3","R","D","N")
yout<-ode(yinit,times,fun,parm)
png(filename="mod1.png",width=1200,height=500,pointsize=16)
par(mfrow=c(1,3))
ym=1.2*N
plot(times,yout[,2],type="l",col="turquoise4",ylim=c(0,ym),ylab="compartment values",xlab="day",main="UK unmitigated")
lines(times,yout[,3],type="l",col="orange")
lines(times,yout[,4],type="l",col="blue4")
lines(times,yout[,5],type="l",col="purple2")
lines(times,yout[,6],type="l",col="brown")
lines(times,yout[,7],type="l",col="gold3")
lines(times,yout[,8],type="l",col="red")
lines(c(100,110),c(ym,ym),col="turquoise4")
text(110,ym,"Susceptible",pos=4)
lines(c(200,210),c(ym,ym),col="orange")
text(210,ym,"Infected",pos=4)
lines(c(100,110),c(0.95*ym,0.95*ym),col="blue4")
text(110,0.95*ym,"Infectious (mild)",pos=4)
lines(c(200,210),c(0.95*ym,0.95*ym),col="purple2")
text(210,0.95*ym,"Infectious (severe)",pos=4)
lines(c(100,110),c(0.9*ym,0.9*ym),col="brown")
text(110,0.9*ym,"Infectious (critical)",pos=4)
lines(c(200,210),c(0.9*ym,0.9*ym),col="gold3")
text(210,0.9*ym,"Recovered",pos=4)
lines(c(200,210),c(0.85*ym,0.85*ym),col="red")
text(210,0.85*ym,"Dead",pos=4)
r<-log(yout[50,3]/yout[45,3])/5
text(200,0.5*ym,paste("r =",round(r,3)))
abline(v=52,lty=2,col="grey")
ym=3e6
plot(times,yout[,2],type="l",col="white",ylim=c(0,ym),ylab="compartment values",xlab="day",main="UK unmitigated (detail)")
#lines(times,yout[,3],type="l",col="orange")
#lines(times,yout[,4],type="l",col="blue4")
lines(times,yout[,5],type="l",col="purple2")
lines(times,yout[,6],type="l",col="brown")
#lines(times,yout[,7],type="l",col="gold3")
lines(times,yout[,8],type="l",col="red")
#lines(c(100,110),c(ym,ym),col="turquoise4")
#text(110,ym,"Susceptible",pos=4)
#lines(c(200,210),c(ym,ym),col="orange")
#text(210,ym,"Infected",pos=4)
#lines(c(100,110),c(0.95*ym,0.95*ym),col="blue4")
#text(110,0.95*ym,"Infectious (mild)",pos=4)
lines(c(200,210),c(0.95*ym,0.95*ym),col="purple2")
text(210,0.95*ym,"Infectious (severe)",pos=4)
lines(c(200,210),c(0.9*ym,0.9*ym),col="brown")
text(210,0.9*ym,"Infectious (critical)",pos=4)
#lines(c(200,210),c(0.9*ym,0.9*ym),col="gold3")
#text(210,0.9*ym,"Recovered",pos=4)
lines(c(200,210),c(0.85*ym,0.85*ym),col="red")
text(210,0.85*ym,"Dead",pos=4)
text(200,0.5*ym,paste("r =",round(r,3)))
abline(v=52,lty=2,col="grey")
ym=1000
plot(times,yout[,2],type="l",col="white",ylim=c(0,ym),xlim=c(0,75),ylab="compartment values",xlab="day",main="UK unmitigated (detail)")
#lines(times,yout[,3],type="l",col="orange")
#lines(times,yout[,4],type="l",col="blue4")
#lines(times,yout[,5],type="l",col="purple2")
#lines(times,yout[,6],type="l",col="brown")
#lines(times,yout[,7],type="l",col="gold3")
lines(times,yout[,8],type="l",col="red")
#lines(c(100,110),c(ym,ym),col="turquoise4")
#text(110,ym,"Susceptible",pos=4)
#lines(c(200,210),c(ym,ym),col="orange")
#text(210,ym,"Infected",pos=4)
#lines(c(100,110),c(0.95*ym,0.95*ym),col="blue4")
#text(110,0.95*ym,"Infectious (mild)",pos=4)
#lines(c(200,210),c(0.95*ym,0.95*ym),col="purple2")
#text(210,0.95*ym,"Infectious (severe)",pos=4)
#lines(c(200,210),c(0.9*ym,0.9*ym),col="brown")
#text(210,0.9*ym,"Infectious (critical)",pos=4)
#lines(c(200,210),c(0.9*ym,0.9*ym),col="gold3")
#text(210,0.9*ym,"Recovered",pos=4)
lines(c(200,210),c(0.85*ym,0.85*ym),col="red")
text(210,0.85*ym,"Dead",pos=4)
text(20,0.5*ym,paste("r =",round(r,3)))
abline(h=281)
abline(v=52)
par(mfrow=c(1,1))
dev.off()
******************************************************
Please note that if you press the 'Reply' button your
message will go only to the sender of this message.
If you want to reply to the whole list, use your mailer's
'Reply-to-All' button to send your message automatically
to [log in to unmask]
Disclaimer: The messages sent to this list are the views of the sender and cannot be assumed to be representative of the range of views held by subscribers to the Radical Statistics Group. To find out more about Radical Statistics and its aims and activities and read current and past issues of our newsletter you are invited to visit our web site www.radstats.org.uk.
*******************************************************
|