1. 3. calculate growth percent and doubling time from covid data
   4. from COVID-19 epidemic data
   5. this code needs adjust if you wanna use it
   6. you must find manually peak date of epidemy
   7. additional script only, meybe needs adjust
   10. 12.01.2021
   11. v 0000.0001
   14. install.packages("ggplot2", "outbreaks", "incidence", "R0", "magrittr", "plotrix")

1. install.packages("ggplot2", "outbreaks", "incidence", "R0", "magrittr", "plotrix", repos ="https://ftp.acc.umu.se/mirror/CRAN/")

Sys.setenv(LANGUAGE="fi") Sys.setlocale("LC_TIME", "fi_FI.UTF-8") Sys.setlocale("LC_ALL","Finnish") options(encoding = "UTF-8")

library(outbreaks) library(ggplot2) library(incidence) library (R0) library(magrittr) library(plotrix) library(svglite)

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

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

plotname<-"Koronavirusepidemian 1 aalto kasvuluvut 1.svg" ## name of plot file

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

1. 1. 1 wiki data, 2 aggregated cases data

load_data_from=2

beginday1='1/03/2020'

1. beginday1='1/03/2020'

endday1='1/07/2020'

1. 1. 1 from finnish wiki, 2 cases from net, 3 from net 2

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

tulosfilee1="/Users/himot/akor1/R0_Suomessa_2.svg"

1. 1. ggolot smooth curves pars

spanni=0.05

1. spanni=0.5

metodi="loess"

1. 1. NOTE date limits change this

1. datelimits1=c('1/3/2020', '9/11/2020')
   1. display date limits

1. today=Sys.Date()-1

1. 1. if above not func try this

today=Sys.Date()

1. print(today)

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

print(today1)

1. print(today2)

1. stop(-1)

1. datelimits1=c(beginday1, today1)

datelimits1=c(beginday1, endday1)

1. 1. dates of dataset

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

paivat1=seq(as.Date("2020/3/1"), as.Date("2020/7/1"), "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)

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)

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'

dfine000 <- read.csv(file=srkurl)

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

1. head(dfinland)

1. print(dfinland)

1. quit(-2)

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)

}

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)

}

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)

# quit(-1)

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

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

#stop(-1)

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

len1=length(cases1)

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

num1<-cases1
names1<-as.Date(dates1)
  

num<- num1
names<-names1

 datess1<-dates2

1. 1. r0 from last week!

start_lok<-len1-7 end_lok<-len1

1. print (names[start_lok])

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

1. lena=length(num)

1. print (lena)

1. num<-num[num<0]=1

df1 <- setNames(num, names)

str(df1)

lena=length(names1)

1. 1. indates<-names1

indates<-1:lena incases<-as.integer(num1) indates<-as.integer(indates)

1. indates<-datess1

indates incases

dat<-0

n <- 1 for (paiva in indates) {

kase<-incases[n] if(kase<0) kase=0 #print (paiva) dat1<-rep(paiva, times = kase) #print (dat1) dat<-append(dat, dat1)

n = n+1 }

1. dat<-as.integer(dat)

class(dat)

1. > [1] "Date"

head(dat)

print("Dunkku")

i <- incidence(dat) i

class(i)

str(i) plot(i)

i.7 <- incidence(dat, interval = 1)

1. fittis<-fit_optim_split(i.7)

1. fitti1<-get_fit(fittis$fit)

Pvm<-datess1[0:140]

d1<-datess1[0:38]

earlylate<-i.7[0:140] early<-i.7[0:38] late<-i.7[38:140]

early.fit <- fit(early) late.fit <- fit(late)

1. plot(fitti1)

plot(earlylate) %>% add_incidence_fit(early.fit) %>% add_incidence_fit(late.fit)

rearly<-round(get_info(early.fit, "r"),3)*100 rlate<-round(get_info(late.fit, "r"),3)*100

doublingearly<-round(get_info(early.fit, "doubling.conf"),1) halvinglate<-round(get_info(late.fit, "halving.conf"),1)

print(rearly) print(doublingearly) print(rlate) print(halvinglate)

str(early.fit)

Tapauksia<-earlylate$counts

1. predik<-get_info(fit(earlylate), "pred", groups = 2)

1. print(predik)

1. plot(predik)

xy<-data.frame(Pvm, Tapauksia)

gr1=rearly/100 td1=log(2.0)/gr1 tau1=5.0 k1=log(2.0)/td1 r01=round(exp(k1*tau1),2)

gr2=rlate/100 td2=log(2.0)/gr2 tau2=5.0 k2=log(2.0)/td2 r02=round(exp(k2*tau2),2)

print(r02)

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

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

plot(Pvm,Tapauksia, type="lines", col="black", lwd=5, ltyt=2, main="Koronaviruksen 1. aalto: kasvu ja hiipuminen", xlab="Kuukausi vuonna 2020", ylab="Tapauksia",

cex.lab=1.5, cex.axis=1.5, cex.main=1.5, cex.sub=1.5)

text(18335, 200, paste0(" Kasvoi \n",rearly, " % pv R0 ",r01), cex = 1.8, col='red') text(18430, 180, paste0(" Hiipui \n",rlate, " % pv"), cex = 1.8, col='blue') text(18335, 150, paste0(" Kaksinkertaistui \n",doublingearly[1],"-",doublingearly[2], " pv"), cex = 1.8, col='darkred') text(18430, 130, paste0(" Puolittui \n",halvinglate[1],"-",halvinglate[2], " pv"), cex = 1.8, col='darkblue') text(18400, 240, paste0(" Alkoi 26. helmikuuta, nousi huippuunsa 38 pvssa"), cex = 1.4, col='brown')

if(levylle==1)

{

	dev.off()

}

}

if(plottaa==2) {

if(levylle==1) { plotname1<-paste0(polku, plotname) # svg(filename=plotname1, width=6, height=3, pointsize=12) }

gee1<- ggplot(xy, aes(x =Pvm , y = Tapauksia)) +geom_point()+

ggtitle("Suomen 1. koronavirusaalto alkuvuosi 2020") +
xlab("Kuukausi vuonna 2020") + ylab("Koronavirustapauksia")+
theme(title=element_text(size=15), axis.text=element_text(size=14,face="bold"),axis.title=element_text(size=14,face="bold"))+ 
geom_line(lwd=1.2,lty=1, colour="#00003f")

gee2=gee1 + geom_text(x=18400, y=250, size=5, label=paste0(" Alkoi 26. helmikuuta, nousi huippuunsa 38 pv"), color="#7f3f2f") gee2=gee2 + geom_text(x=18400, y=220, size=5, label=paste0(" Huippu 4. huhtikuuta , 267 tapausta"), color="#7f3f2f")

gee2=gee2 + geom_text(x=18335, y=200, size=6.0, label= paste0(" Kasvoi ",rearly, " % pv\nR0 ", r01), color="#af0000") gee2=gee2 + geom_text(x=18330, y=130, size=6.0, label= paste0(" Kaksinkertaistui \n",doublingearly[1],"-",doublingearly[2], " pv"), color="#7f0000")

gee2=gee2 + geom_text(x=18430, y=180, size=6, label= paste0(" Hiipui ",rlate, " % pv, R0 ",r02), color="#0000af") gee2=gee2 + geom_text(x=18430, y=130, size=6, label=paste0(" Puolittui ",halvinglate[1],"-",halvinglate[2], " pv") , color="#0000af")

plot(gee2)

if(levylle==1)

{

ggsave(file=plotname1, plot=gee2, width=10, height=4)		
#	dev.off()

} }