1. 3. Calculate r0 of Covid-19 in Finland
   4. "R" script
   6. 23.7.2022
   7. v 0000.0013

new_in_skript=0

if (new_in_skript==1) { #install.packages("ggplot2", "plotly", repos ="https://ftp.acc.umu.se/mirror/CRAN/") install.packages("svglite") install.packages("ggplot2")

install.packages("rvest") install.packages("readtext") install.packages("stringi") install.packages("datamart") install.packages("XML")

install.packages("tidyr") install.packages("stringr") install.packages("stringi") install.packages("tibble") install.packages("readr") install.packages("data.table")

install.packages("caTools") install.packages("mgcv") install.packages("repmis") install.packages("lubridate") install.packages("tidyverse") install.packages("R0") install.packages("EpiEstim") install.packages("jsonlite") install.packages("rjstat") }

library(ggplot2) library(svglite)

library(rvest) library(readtext) library(stringi) library(stringr)

1. library(datamart)

library(XML) library(jsonlite)

library(tibble) library(caTools) library(mgcv) library(repmis) library(lubridate) library(tidyverse) library(tidyr) library(readr) library(data.table) library(rjstat)

1. library (R0)

library(EpiEstim)

1. 1. choices
   2. 1 finnish wiki data !!! NOT UPDATED
   3. 2 aggregated cases data OK
   4. 3 solanpaa finnish data OK
   5. 4 thl cube json data NOK
   6. 5 thl cube json hospital data OK 15.4.2022

load_data_from=3

beginday1='28/02/2020'

1. beginday1='1/7/2021'

1. beginday1='1/1/2021'

1. today=Sys.Date()-1

today=Sys.Date()-4

spanni=0.1

1. spanni=0.3

yala=0.6 yyla=2.0

1. yala=0.9
2. yyla=1.35

metodi="loess"

widthi=10 heighti=5

plottaa=1 ## must be 1 tulosta_svg=1 # plot to out svg 0, 1 of 2

tulosfilee1="./R0_Suomessa_2.svg"

beginday0=as.Date(beginday1) beginday2=format(beginday0, "%Y/%m/%d")

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

print(today1)

datelimits1=c(beginday1, today1)

paivat1=seq(as.Date(beginday2), as.Date(today2), "days")

calculate_r0 <- function(time1, time2, val1, val2) { td=time2-time1 gr0<-log(val2/val1) gr=gr0/td td = log(2)/gr tau<-5.0 k<-log(2.0)/td r0<-exp(k*tau) return(r0) }

moving_average <- function(x, w, FUN, ...) {

 if (w < 1) {
   stop("Window length: mustbe greater than 0")
 }
 output <- x
 for (i in 1:length(x)) {
   lower_bound <- i - w + 1
   if (lower_bound < 1) {
     output[i] <- NA_real_
     ## !!! assume NA 0
     output[i] <- 0
   } else {
     output[i] <- FUN(x[lower_bound:i, ...])
   }
 }
 return (output)

}

calculate_multiple_r0 <- function(daata1) {

 lenu1<-length(daata1)
  
 daata2<-1:lenu1
 

for (n in 2:lenu1){ valju1=daata1[n-1] valju2=daata1[n] timex1=0 timex2=1

r0<-calculate_r0(0, 1, valju1, valju2) daata2[n]<-r0

#print (r0) } return(daata2)

}

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)

pv0<-tabu1$V1 len1=length(pv0) daates1 <- vector(mode="character", length=len1)

1. print(pv0)

n=1

for(n in 1:len1) { it1<-pv0[n] #print(it1) qq1<-str_split(it1, "\\[")1 qq2<-qq1[1] qq3<-gsub(" ", "", qq2, fixed = TRUE) daates1[n]=qq3 }

daates2=as.Date(daates1, format="%d.%m.%Y")

print(daates2)

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(daates2, sairaalassa)

xy<-data.frame(daates2, uusia_tapauksia) names(xy)<-c("Dates", "Cases")

xyz<-data.frame(daates2, sairaalassa, teholla) dfout1<-data.frame(daates2, 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

print("Aggreg")

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

srkurl='https://raw.githubusercontent.com/datasets/covid-19/main/data/countries-aggregated.csv'

dfine <- read.csv(file=srkurl)

1. print(dfine000)

1. str(dfine000)

1. head(dfine000,20)

1. stop(-1)

1. dfine<-as.data.frame(dfine000)

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

1. print(dfine)

1. tail(dfine)

1. stop(-1)

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

head(dfinland)

1. print(dfinland)

1. stop(-1)

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')

datez1<-dfout2['Date'] dailyz1<-dfout2['DailyConfirmed'] deathz1<-dfout2['DailyDeaths']

dfout3<-cbind(datez1, dailyz1)

return(dfout3)

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

pv0<-tabu1$V1 len1=length(pv0) daates1 <- vector(mode="character", length=len1)

1. print(pv0)

n=1

for(n in 1:len1) { it1<-pv0[n] #print(it1) qq1<-str_split(it1, "\\[")1 qq2<-qq1[1] qq3<-gsub(" ", "", qq2, fixed = TRUE) daates1[n]=qq3 }

daates2=as.Date(daates1, format="%d.%m.%Y")

print(daates2)

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

1. str(dataf1)

  dataf2 <- data.frame("Date" = daates2, "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(daates2, dailycases3)


return(xy)

}

download_solanpaa_finnish_data<-function() { solanpaa_fi="https://covid19.solanpaa.fi/data/fin_cases.json" cache_file="solanpaa_fi.json"

download.file(solanpaa_fi, cache_file)

j1 <- fromJSON(cache_file)

 ## maybe errori

dates<-as.Date(j1$date)

dailycases<-j1$new_cases dailydeaths<-j1$new_deaths

   dataf1 <- data.frame("Date" = dates, "Paivitt_kuolemat"=dailydeaths) 
   dataf2 <- data.frame("Date" = dates, "Paivitt_tapaukset"=dailycases)

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

xy0<-data.frame(dates, dailycases) names(xy0)<-c("Dates", "Cases") xy<-na.omit(xy0)

return(xy)

}

calculate_r0_with_r0<-function(xy2) {

## calculate r0 w/r0 package dates<-as.Date(xy2$Dates) cases<-as.integer(xy2$Cases)

cases[is.na(cases)] <- 1 cases[(cases<0)] <- cases*-1 cases[cases==0] <- 1

nummeros<-1:length(dates) num<-cases #names<-nummeros names<-dates lenu=length(dates)

bekini=as.Date(dates[1]) enti=as.Date(dates[lenu])

#print(bekini) #print(enti)

#stop(-1)

#enti=lenu

#bekini=enti*0+1

#enti=as.integer(enti) #bekini=as.integer(bekini)

df1 <- setNames(num, names)

mGT<-generation.time("gamma", c(3, 1.5)) #TD <- est.R0.TD(df1, mGT, begin=1, end=length(dates), nsim=200) #TD <- est.R0.TD(df1, mGT, begin=bekini, end=enti, nsim=200) TD <- est.R0.TD(df1, mGT, begin=bekini, end=enti, nsim=200) TD.5D <- smooth.Rt(TD, 5) paivat1<-TD.5D$epid$t paivat2<-as.Date(paivat1) r0t1<-TD.5D$R conf1<-TD.5D$conf.int

xypaluu<-data.frame(paivat1,r0t1) names(xypaluu)<-c("paivat","r0") return(xypaluu) }

calculate_r0_with_epiestim<-function(xy2) {

## calculate r0 w/r0 package dates<-as.Date(xy2$Dates) cases<-as.integer(xy2$Cases) nummeros<-1:length(dates) num<-cases #names<-nummeros names<-dates lenu=length(dates)

cases[is.na(cases)] <- 1 cases[(cases<0)] <- cases*-1 cases[cases==0] <- 1

incid<-cases

bekini=as.Date(dates[1]) enti=as.Date(dates[lenu])

config<-make_config( list(mean_si = 2.6,std_si = 1.5) )

res<-estimate_R(incid,method="parametric_si", config = config)

plot(res) resr<-res$R

str(resr)

meanr<-resr$Mean medianr<-resr$Median quantile95<-resr$Quantile.0.95 quantile05<-resr$Quantile.0.05 quantile75<-resr$Quantile.0.75 quantile25<-resr$Quantile.0.25 meanr

daydexes<-resr$t_start

daydexes

plot(daydexes, meanr)

dayss<-as.Date(dates[daydexes])

print (dayss) #stop(-1) plot(dayss, meanr)

xypaluu<-data.frame(dayss,meanr) names(xypaluu)<-c("paivat","r0") return(xypaluu) }

calculate_r0_with_simple_exponent_moving_average<-function(xy2, madays1, madays2) {

## calculate r0 w/r0 package dates<-as.Date(xy2$Dates) cases<-as.integer(xy2$Cases) nummeros<-1:length(dates) num<-cases #names<-nummeros names<-dates lenu=length(dates)

cases[is.na(cases)] <- 1 cases[(cases<0)] <- cases*-1 cases[cases==0] <- 1

# compute a MA(7) ma1<-moving_average(cases,madays1,mean) r0t1<-calculate_multiple_r0(ma1) r0avg1<-moving_average(r0t1, madays2, mean) xypaluu<-data.frame(dates,r0t1)

plot(r0t1) print (r0t1) #stop(-1)

names(xypaluu)<-c("paivat","r0")

return(xypaluu)

}

lataa_thl_tapaukset_kuolleet<-function() { ## oriko

   ##url1<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=measure-492118&column=dateweek20200101-508804L"

# viikoittain

url1<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=hcdmunicipality2020-445222.&column=dateweek20200101-509030&filter=measure-444833" cube1 <- fromJSONstat(url1, naming = "label", use_factors = F, silent = T) head(cube1)

res01 <- cube11

head(res01,40) print("Pazka") stop(-1)

#res00 url2<-"https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19case/fact_epirapo_covid19case.json?row=measure-444833&column=dateweek20200101-508804L" cube2 <- fromJSONstat(url2, naming = "label", use_factors = F, silent = T) res02 <- cube21 #res02

#stop (-1) paiva=as.Date(res01$dateweek20200101) kuolleet=as.integer(res01$value) tapaukset=as.integer(res02$value)

kuolin_prosentit=kuolleet/tapaukset kuolin_prosentit=kuolin_prosentit*10000 kuolin_prosentit=as.integer(kuolin_prosentit) kuolin_prosentit=as.double(kuolin_prosentit) kuolin_prosentit=kuolin_prosentit/100.0

#print (paiva) #print (kuolleet) #stop(-1) #print (tapaukset) #print (kuolin_prosentit )

df1<-data.frame(paiva,tapaukset, kuolleet, kuolin_prosentit)

names(df1)<-c("Paiva", "Tapauksia", "Kuolleita", "Kuolinprosentti") #write.csv2(df1, "./kuolleet_ikaryhmittain.csv", sep = ";" )

write.csv(df1, "./thl_tapaukset_kuolleet.csv")

xy0<-data.frame(paiva, tapaukset) names(xy0)<-c("Dates", "Cases") xy<-na.omit(xy0)

#return(df1)

}

nth_element <- function(vector, starting_position, n) {

 vector[seq(starting_position, length(vector), n)] 
 }   

get_thl_hospital_data<-function() {

url_base2="https://sampo.thl.fi/pivot/prod/fi/epirapo/covid19care/fact_epirapo_covid19care.json" request2 <- "?row=dateweek20200101-509093L&column=measure-547523.547516.547531.456732.&fo=1"

url2 <- paste0(url_base2, request2) cube2 <- fromJSONstat(url2, naming = "label", use_factors = F, silent = T)

1. 1. print (cube2)

1. head(cube2, 40)

1. stop(-1)

res02 <- cube21

head(res02, 40)

daates00<-res021

1. print(daates00)

days0<-as.Date(daates00) days1<-nth_element(days0, 1, 4)

1. print(days1)

gecko1<-as.integer(res023)

1. print(head(gecko1,40))

gecko2<-matrix(gecko1,nrow=4)

1. print (head(gecko2,20))

sairaalassa1=gecko2[1,]+gecko2[2,]+gecko2[4,]

teholla1=gecko2[4,]

sairaalassa1[is.na(sairaalassa1)]<-0 teholla1[is.na(teholla1)] <- 0

print (head(sairaalassa1,50))

1. 1. print (teholla1)

df1<-data.frame(days1, sairaalassa1, teholla1) names(df1)<-c("Paiva", "Sairaalassa", "Teholla") df0<-data.frame(days1, sairaalassa1) names(df0)<-c("Paiva", "Sairaalassa")

write.csv(df0, "./thl_sairaalassa.csv")

xy0<-data.frame(days1, sairaalassa1) names(xy0)<-c("Dates", "Cases") xy<-na.omit(xy0)

return(xy)

1. return(df0)

}

1. 1. 1. main program

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

if(load_data_from==2) { xy<-load_data_from_aggregated() names(xy)<-c("Dates","Cases")

print("Aggreg data") ##stop(-1)

}

if(load_data_from==3) { xy<-download_solanpaa_finnish_data() }

if(load_data_from==4) { xy<-lataa_thl_tapaukset_kuolleet() }

if(load_data_from==5) { xy<-get_thl_hospital_data() names(xy)<-c("Dates","Cases") print("From THL hospital data")

#print(xy) #head(xy)

#stop(-1)

}

1. print (xy)

## quit(-1)

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

#format(select_datelimit_begin, "%Y-%m-%d")

print(select_datelimit_begin)

print(select_datelimit_end)
#2020-12-16

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

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

print("xy2")
print(xy2)

#stop(-1)

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

len1=length(cases1)

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

## test code
arrat0<-calculate_r0_with_simple_exponent_moving_average(xy2, 14,7)

#arrat1<-calculate_r0_with_r0(xy2)
#arrat2<-calculate_r0_with_epiestim(xy2)

#print("calcu ok")

 #plot(arrat$paivat, arrat$r0)

# arrat<-arrat2
arrat<-arrat0



str(arrat)
head(arrat)

 #plot(arrat$paivat, arrat$r0)

1. stop(-1)

if(tulosta_svg==1) {

#svg(filename=tulosfilee1, width=8, height=3, pointsize=12) svg(filename=tulosfilee1, width=widthi, height=heighti, pointsize=12)

}

if(plottaa==1) {

 metodi="loess"

print ("ggplot")

#ggplot(arrat, aes(x =paivat , y = r0)) +ylim(0.6, 1.8)+xlim(as.Date(datelimits1, format="%d/%m/%Y") )+ ggplot(arrat, aes(x =paivat , y = r0)) +ylim(yala, yyla)+xlim(as.Date(datelimits1, format="%d/%m/%Y") )+ ggtitle("Arvioitu koronaviruksen perusuusiutumisluku R0") + xlab("Kuukausi") + ylab("R0")+ theme(title=element_text(size=15), axis.text=element_text(size=12,face="bold"),axis.title=element_text(size=14,face="bold"))+ geom_point() + geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.99, size=3)+ geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) + geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5) + geom_hline(yintercept=1.0, linetype="dashed", color = "red", size=1)

}

if(tulosta_svg==1) { dev.off() }