1. 3. experimental calculation of R0
   4. R code
   6. calculate R0 from COVID-19 epidemic data
   7. additional script only, meybe needs adjust
   9. v 0001.0000
   10. 11.1.2021

1. install.packages("ggplot2", "plotly", repos ="https://ftp.acc.umu.se/mirror/CRAN/")
2. install.packages("R0",repos ="https://ftp.acc.umu.se/mirror/CRAN/")

1. library(plotly)

library (R0) library (ggplot2)

plottaa=1 ## 0, 1 or 2 -1 debug plotname="R0_Suomessa_1.svg" polku<-"/Users/himot/akor1/"

library(svglite)

library(rvest) library(readtext) library(stringi) library(stringr) library(datamart) library(XML)

beginday1='01/04/2020'

1. 1. limits of input data

datelimits1=c('2020/04/01', '2020/01/07')

1. 1. dates of prediction

datelimits2=c('2020/04/01', '2021/02/15')

load_data_from=2

1. 1. 0, 1 or 2 , 0 none plot, 1 plot, 2 ggplot, -1 debug

plottaa<-1 levylle<-1 ## 0 to display, 1 to disc

plotname<-"R0_Suomessa_4.svg" ## name of plot file

polku<-"/Users/himot/akor1/" ## path of plot file

today=Sys.Date()

1. 1. jos ylläoleva ei toimi, niin tää
   2. if above not func, this
2. today=Sys.Date()-1

1. print(today)

today1=format(today, "%d/%m/%Y") today2=format(today, "%Y/%m/%d")

1. print(today1)
2. print(today2)

1. stop(-1)

datelimits1=c(beginday1, today1)

paivat1=seq(as.Date("2020/4/1"), as.Date(today2), "days")

load_data_from_finnish_wiki<-function() {

url1="https://fi.wikipedia.org/wiki/Suomen_koronaviruspandemian_aikajana" destfile1="./ward0.txt"

download.file(url1, destfile1) texti000<-readtext(destfile1) texti0<-texti000$text

etsittava1="1. huhtikuuta 2020 alkaen" len1=nchar(texti0) k1=regexpr(pattern=etsittava1, texti0) k1b=len1-k1 texti1=strtail(texti0,k1b) sink("out1.txt") print (texti1) sink()

k2=regexpr(pattern=etsittava2, texti1) texti2=strhead(texti1,k2)

sample1<-minimal_html(texti2) tabu1 <- html_table(sample1, fill=TRUE)1 colnames(tabu1) <- c("V1","V2", "V3","V4", "V5","V6", "V7","V8" )

1. print(tabu1)

sairaalassa00<-tabu1$V4 sairaalassa=as.integer(sairaalassa00)

teholla00<-tabu1$V5 teholla=as.integer(teholla00)

uusiatapauksia00<-tabu1$V3 uusiatapauksia0<-gsub(" ", "", uusiatapauksia00) uusia_tapauksia=as.integer(uusiatapauksia0)

uusiakuolleita00<-tabu1$V7 uusiakuolleita1=as.integer(uusiakuolleita00)

uusiakuolleita2<-uusiakuolleita1 uusiakuolleita2[uusiakuolleita2<0]<-0 uusia_kuolleita<-uusiakuolleita2

toipuneita00<-tabu1$V8 toipuneita01<-gsub(" ", "", toipuneita00) toipuneita0<-gsub("[^0-9.-]", "", toipuneita01) toipuneita=as.integer(toipuneita0)

tapauksia00<-tabu1$V2 tapauksia01<-gsub(" ", "", tapauksia00) tapauksia0<-gsub("[^0-9.-]", "", tapauksia01)

tapauksia=as.integer(tapauksia0)

kuolleita00<-tabu1$V6 kuolleita=as.integer(kuolleita00)

aktiivisia_tapauksia=tapauksia-kuolleita-toipuneita

1. print (paivat1)
2. print (teholla)
3. print (sairaalassa)
4. print (tapauksia)
5. print (kuolleita)
6. print (toipuneita)
7. print (uusia_tapauksia)
8. print (uusia_kuolleita)
9. plot(paivat1,aktiivisia_tapauksia)

1. xy<-data.frame(paivat1, sairaalassa)

xy<-data.frame(paivat1, uusia_tapauksia)

1. xy<-data.frame(paivat1, tapauksia)

xyz<-data.frame(paivat1, sairaalassa, teholla) dfout1<-data.frame(paivat1, aktiivisia_tapauksia, uusia_tapauksia, sairaalassa, teholla, uusia_kuolleita )

names(dfout1)<-c("Pvm", "Aktiivisia_tapauksia","Uusia_tapauksia", "Sairaalassa", "Teholla", "Uusia_kuolleita")

write.csv2(dfout1, "./sairaalassa.csv",row.names=FALSE )

return(xy) }

load_data_from_aggregated<-function() {

1. 2. fetch the data

srkurl='https://datahub.io/core/covid-19/r/countries-aggregated.csv'

dfine <- read.csv(file=srkurl)

1. head(dfine)
2. class(dfine)

1. tail(dfine, 5)

dfinland <- dfine[ which(dfine$Country=='Finland'), ]

1. head(dfinland)

kols <- c("Date", "Confirmed","Recovered","Deaths")

tapaukset <- dfinland[kols]

1. head(tapaukset)

len1=nrow(tapaukset)

1. len1

len2=len1-1

len3=len2

confirmed<-tapaukset$Confirmed deaths<-tapaukset$Deaths

dailycases <- vector() dailycases <- c(dailycases, 0:(len2)) dailydeaths <- vector() dailydeaths <- c(dailydeaths, 0:(len2))

m=0 dailycases[1]<-tapaukset$Confirmed[1]

1. dailydeaths[1]<-tapaukset$Deaths[1]

dailydeaths[1]<-0

1. confirmed
2. deaths

m=1 for(n in 2:(len3+1)) {

a<-confirmed[n] b<-confirmed[m] #print (a) #print (b) cee<- (a-b) #print(cee) dailycases[n]=cee m=m+1 }

mm=1 for(nn in 2:(len3+1)) {

aa<-deaths[nn] bb<-deaths[mm] #print ("_") #print (aa) #print (bb) ceb=aa-bb #if (ceb<0) ceb=0 #print(ceb) dailydeaths[nn]=ceb mm=mm+1 }

1. deaths

1. dailycases

1. dailydeaths

dfout1<-dfinland

1. print(nrow(dfinland))
2. print(length(dailydeaths))

dfout1 <- cbind(dfout1, data.frame(dailycases)) dfout1 <- cbind(dfout1, data.frame(dailydeaths))

1. head(dfout1)

dfout2<-within(dfout1, rm(Country))

names(dfout2) <- c('Date','Confirmed','Recovered','Deaths', 'DailyConfirmed','DailyDeaths')

1. head(dfout2)

write.csv2(dfout2, "/Users/himot/akor1/finland_data1.csv");

daate1<-dfout2$Date dailydeaths1<-dfout2$DailyDeaths dailycases1<-dailycases

1. daate1

1. daate2<-gsub("2020-", "", daate1)

daate2<-daate1

leenu<-length(daate2)

1. alkupvm<-50

alkupvm<-1

daate3<-daate2[alkupvm:leenu] dailydeaths3<-dailydeaths1[alkupvm:leenu] dailycases3<-dailycases1[alkupvm:leenu]

1. daate3
2. dailydeaths3

# barplot(dailydeaths3, main="Koronaviruskuolemat päivittäin vuonna 2020",
# names.arg=daate3) 
 
     
 dataf1 <- data.frame("Date" = daate3, "Paivitt_kuolemat"=dailydeaths3)

1. str(dataf1)

  dataf2 <- data.frame("Date" = daate3, "Paivitt_tapaukset"=dailycases3)

1. str(dataf2)

write.csv(dataf1, "/Users/himot/akor1/dailydeaths1.csv", row.names=T)
write.csv(dataf2, "/Users/himot/akor1/dailycases1.csv", row.names=T)

indf1 <- read.csv(file = '/Users/himot/akor1/dailycases1.csv')
#head(indf1)

cases1<-indf1$Paivitt_tapaukset
dates1<-indf1$Date

len1=length(cases1)

dates2<-as.Date(dates1)
paivat<-1:len1	

xy<-data.frame(daate3, dailycases3)

}

if(load_data_from==1) { xy<-load_data_from_finnish_wiki() }

if(load_data_from==2) { xy<-load_data_from_aggregated() }

names(xy)<-c("Dates","Cases")

print(xy)

#stop(-1)


select_datelimit_begin=as.Date(beginday1,format="%d/%m/%Y")
select_datelimit_end=as.Date(today1)

xy2<-xy[xy$Dates >= select_datelimit_begin,]

#print(xy2)

1. stop(-1)

cases1<-xy2$Cases
dates1<-xy2$Dates

xy3<-data.frame( as.Date(dates1),as.integer(cases1) )
names(xy3)<-c("Dates", "Cases")

len1=length(cases1)

dates2<-as.Date(dates1)
paivat<-1:len1

num1<-cases1 dates1<-dates1 names1=dates1

len1=length(num1)

1. 1. sure negative values to zero

num1[num1<0]<-0

1. 1. r0 from last week!

lensub1=105 lensub1=7*17+2 lensub3=7*24+2

start_lok<-len1-lensub1

1. start_lok<-len1-7

end_lok<-len1

print (names1[start_lok])

num<-num1[start_lok:end_lok] names<-names1[start_lok:end_lok]

lena=length(num)

print (lena)

df1 <- setNames(num, names)

str(df1)

1. generation time distribution

mGT = generation.time ("gamma", c(3, 1.5))

1. 1. r0, exponential method, last week:

est1<-est.R0.EG (df1, mGT, begin=1, end=77)

mGT = generation.time("gamma", c(2.45, 1.38))

1. 1. r0, second method

est2<-est.R0.ML (df1, mGT)

est1 est2

class(est1)

est1[1] est2[1]

r0exp1<-est1[1] r0exp1

estimaatti1<-r0exp1[1]

class(r0exp1) class(estimaatti1)

str(estimaatti1)

vaalu1<-estimaatti1[1]$R

print ("vaalu1") print (vaalu1[1])

1. jono1<-toString(round(as.Numeric(r0exp1),4))

valju1<-round(vaalu1,2)

jono1<-toString(valju1) jono11<-paste("R0_exp (kulunut viikko) ",jono1)

1. barplot(dailycases3, main="Koronavirustapaukset päivittäin vuonna 2020",
2. sub=jono11,
3. names.arg=daate3)

mGT<-generation.time("gamma", c(3, 1.5)) TD <- est.R0.TD(df1, mGT, begin=1, end=lensub1, nsim=200)

TD.weekly <- smooth.Rt(TD, 7) TD.weekly

TD.5D <- smooth.Rt(TD, 5)

paivat1<-TD.5D$epid$t paivat2<-as.Date(paivat1) r0t1<-TD.5D$R conf1<-TD.5D$conf.int

class(TD.5D$conf.int)

if(plottaa==-1) { plot(paivat2, r0t1, pch=4, main="Arvioitu R0 Suomessa", xlab="Päivä", ylab="R0") lines(paivat2,r0t1, col="black", lwd=4) lines(paivat2,conf1$upper, col="red", lwd=1) lines(paivat2,conf1$lower, col="blue") }

if (plottaa==0) { plot(TD.5D, main="R0", xlab="Päivä", ylab="R0") }

if(plottaa==1) { print("Plot 1 ,,,") plotname1<-paste0(polku, plotname) print(plotname1) svg(filename=plotname1, width=9, height=5, pointsize=12)

plot(paivat2, r0t1, pch=20, main="Arvioitu R0 Suomessa", xlab="Kuukausi", ylab="R0", ylim=c(0.3,2.0), cex.lab=1.3, cex.axis=1.3, cex.main=1.3, cex.sub=1.3) abline(h=1.0, col="green", lty=2, lwd=2) lines(paivat2,r0t1, col="black", lwd=4) lines(paivat2,conf1$upper, col="red", lwd=1) lines(paivat2,conf1$lower, col="blue")

dev.off() }

if(plottaa==2) { plotname1<-paste0(polku, plotname) svg(filename=plotname1, width=6, height=3, pointsize=12)

plot(TD.5D, main="R0", xlab="Päivä", ylab="R0")

dev.off() }

1. plot(TD)

1. plot(TD.weekly,type = "o", col = "red", xlab = "Viikko", ylab = "R0", main="Koronaviruksen R0 Suomessa")

1. TD.weekly$R[1]