1. 2. COVID-19 data refch && r0 calculation
      1. with EpiEstim package
      2. data fetch from github COVID-19 data site
   4. additional script only, maybe needs adjust
   7. v 0001.0002
   8. 30.5.2021
   10. cases data smoothed with GAM

1. 1. install.packages("EpiEstim",repos ="https://ftp.acc.umu.se/mirror/CRAN/")

library(EpiEstim) library(mgcv) library(ggplot2) library(svglite)

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

beginday1='01/04/2020'

1. 1. limits of input data

datelimits1=c('2020/04/01', '2021/05/28')

1. 1. dates of prediction

datelimits2=c('2020/04/01', '2021/05/28')

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

load_data_from=3

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

polku<-"../" ## path of plot file

today=Sys.Date()

1. 1. jos ylläoleva ei toimi, niin tää
   2. if above not func, this

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

dfine <- read.csv(file = 'https://datahub.io/core/covid-19/r/countries-aggregated.csv')

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)

}

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)

}

lataa_thl_tapaukset_kuolleet<-function() {

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

cube1 <- fromJSONstat(url1, naming = "label", use_factors = F, silent = T) res01 <- cube11 #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)

}

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

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

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

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

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

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)

alku<-as.integer(len1*0+1) loppu<-as.integer(len1) start_lok<-as.integer(1) end_lok<-as.integer(len1)

paivat<-as.Date(dates1[start_lok:(end_lok)]) indexes1=start_lok:end_lok

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

num<-num1[start_lok:end_lok]

1. names<-names1[start_lok:end_lok]
2. names<-indexes1

dates<-paivat

1. dates<-names1

I<-num1

1. ceises <- setNames(dates, num)

ceises<-data.frame(dates, I) incid<-I

1. smooth data

x=indexes1 y=incid

#3. deg poly

1. lo<- loess(y~x)

1. 1. 10 asteen poly ok

fit3<-gam(y~poly(x,10)) xx <- seq(1,length(incid), length=length(incid))

1. fittedvalues<- fit3$fitted.values

1. str(fit3)

fity1<-predict(fit3, data.frame(x=xx)) fitu1<-fity1 names(fitu1) <- NULL

fittedvalues1<-fitu1 fittedvalues1

str(fity1)

I<-NULL

fittedvalues1<-as.integer(fittedvalues1) fittedvalues1[fittedvalues1<1] <- 1

I<-fittedvalues1

ceises<-data.frame(dates, I) incid<-I

1. fity2<-predict(lo)

1. names(ceises)<-c("dates" "I")

1. 1. dates, I

1. 1. parametric si

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

1. 1. estimete R

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

plot(res)

1. plot(res, "R")

1. 1. uncertain si

1. onfig <- make_config(list(mean_si = 2.6, std_mean_si = 1,
2. min_mean_si = 1, max_mean_si = 4.2,
3. std_si = 1.5, std_std_si = 0.5,
4. min_std_si = 0.5, max_std_si = 2.5))

1. es2 <- estimate_R(incid,method = "uncertain_si",config = config)

1. lot(res2)

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

dayss<-paivat[daydexes]

if(plottaa==1) { if(levylle==1) { plotname1<-paste0(polku, plotname) svg(filename=plotname1, width=16, height=5, pointsize=12) }

1. plot(dayss, meanr, main="Koonaviruksen R0 Suomessa")

time=as.POSIXct(dayss) y=meanr

plot(time, y, xaxt="n", type="lines", col="black", lwd=3, main="Arvioitu koronaviruksen R0 Suomessa", xlab="Kuukausi 2020-2021 ...", ylab="Tartutusluku R0",

ylim=c(0.6,1.8), cex.lab=1.6, cex.axis=1.6, cex.main=1.5, cex.sub=1.5)

axis.POSIXct(side=1,at=cut(time, breaks="month"),format="%b", cex.lab=1.6, cex.axis=1.6, cex.main=1.5, cex.sub=1.5) abline(h=1.0, col="#007f00", lty=2, lwd=3) lines(time,meanr, col="black", lwd=4) lines(time,medianr, col="#004f00", lwd=2) lines(time,quantile95, col="red", lty=2) lines(time,quantile05, col="blue", lty=2) lines(time,quantile75, col="red", lty=1) lines(time,quantile25, col="blue", lty=1)

if(levylle==1)

{

	dev.off()

}

}

if(plottaa==-1) {

plot(x,y,pch=19)

print(length(xx)) print(length(y)) print (length(paivat))

plot(paivat,y,pch=4,ylim=c(0,1000))

lines(paivat,fity1 , col="blue", lwd=3, lty=2)

#lines(paivat,fity2 , col="red", lwd=3)

}