1. 2. R script
   3. Covid-19 number of cases predict/prediction/forecast
   4. Gompertz curve
   5. finland
   7. 07.12.2021
   8. 0000.0004

1. 1. choices
   2. 1 finnish wiki data, 2 aggregated cases data
   3. 3 solanpaa finnish data 4 thl cube json data

load_data_from=3

beginday1='25/06/2022'

1. 1. limits of input data

datelimits1=c('2022/06/25', '2022/07/07')

1. 1. dates of prediction

datelimits2=c('2022/06/25', '2022/10/01')

1. beginday1='21/03/2021'
   1. limits of input data
2. datelimits1=c('2021/03/21', '2021/04/24')
   1. dates of prediction
3. datelimits2=c('2021/03/21', '2021/06/15')

1. install.packages("growthmodels", repos ="https://ftp.acc.umu.se/mirror/CRAN/")
2. install.packages("minpack.lm", repos ="https://ftp.acc.umu.se/mirror/CRAN/")

library(stats) library(growthmodels) library(minpack.lm) library(ggplot2) library(svglite)

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

1. library(datamart)

library(XML) library(jsonlite) library(rjstat)

library(tibble) library(caTools) library(mgcv) library(repmis) library(lubridate) library(tidyverse) library(tidyr) library(dplyr)

1. library(covid19.analytics)
2. library(R0)

library(EpiEstim) library(prophet) library(scales)

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

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)

}

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)

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

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 cases<-j1$cases

   #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, "./cases1.csv", row.names=T)

xy0<-data.frame(dates, cases) 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

# compute a MA(7)
ma1<-moving_average(cases1, 14,mean)

1. fit data to the exponential model

1. x=c(1,2,3,4,5,6,7,8,9,10)
2. y=c(1,2,3,4,5,6,7,8,9,10)

pitu=length(cases1)-1

1. print (cases1)

x<-seq(0,pitu)-1

y<-cases1

print (length(x)) print (length(y))

cases=cases1 days=x

dates3=seq(as.Date(datelimits2[1]), as.Date(datelimits2[2]), "days") len2=length(dates3)

len1=length(cases)

days=seq(1,len1)

days.predict=seq(1, len2)

1. print (days)

1. print (length(days))
2. print (length(cases))

1. days<-c(1,2,3,4,5,6)
2. cases<-c(1,2,4,9,11,12)

1. days.predict<-c(1,2,3,4,5,6)

data<-data.frame(days, cases)

1. str(data)

1. stop(-1)

alpha = 9526 beta = 9.1618 k = 0.0028 nls.gompertz <- minpack.lm::nlsLM(data$cases~alpha*exp(-beta*exp(-k*data$days)), data = data, start = list(alpha = alpha, beta = beta, k = k), control = list(maxiter = 500)) coef(nls.gompertz) ## alpha = 9437, beta = 59.24, k = 0.0219

1. 1. Now fit Geompertz model

alpha0 = coef(nls.gompertz)"alpha" beta0 = coef(nls.gompertz)"beta" k0 = coef(nls.gompertz)"k"

alpha1=alpha0 beta1=beta0+0.00 k1=k0/2.5

alpha2=alpha0 beta2=beta0+0.00 k2=k0*2.5

alpha3=alpha0 beta3=beta0+0.00 k3=k0/1.25

alpha4=alpha0 beta4=beta0+0.00 k4=k0*1.25

growth.gompertz <- growthmodels::gompertz(data$days, alpha = coef(nls.gompertz)"alpha", beta = coef(nls.gompertz)"beta", k = coef(nls.gompertz)"k") growth.gompertz

1. 1. Predict

predict.gompertz <-growthmodels::gompertz(days.predict, alpha = coef(nls.gompertz)"alpha", beta = coef(nls.gompertz)"beta", k = coef(nls.gompertz)"k")

predict.gompertz1 <-growthmodels::gompertz(days.predict, alpha = alpha1, beta = beta1, k=k1) predict.gompertz2 <-growthmodels::gompertz(days.predict, alpha = alpha2, beta = beta2, k=k2)

predict.gompertz3 <-growthmodels::gompertz(days.predict, alpha = alpha3, beta = beta3, k=k3) predict.gompertz4 <-growthmodels::gompertz(days.predict, alpha = alpha4, beta = beta4, k=k4)

print("abk") print (alpha) print (beta) print (k)

1. stop(-1)

1. predict.gompertz

1. str(predict.gompertz)

x=days y=cases

x2=days.predict y2=predict.gompertz

ya1=predict.gompertz1 ya2=predict.gompertz2 ya3=predict.gompertz3 ya4=predict.gompertz4

deltalen1=len2-len1

cy3=1:deltalen1 cy3=cy3*0

y3=c(y,cy3)

1. print ("qqqq")
2. print (length(y3))

1. print (x2)
2. print (y2)
3. print (dates3)

1. print (length(x2))
2. print (length(y2))
3. print (length(dates3))

1. print (y3)

xy=data.frame(x,y) xy2=data.frame(x2,y2)

xy3=data.frame(x2,y3)

xy4=data.frame(x2,dates3, y2,y3, ya1, ya2)

1. plot(x2,y2)
2. lines(x,y)

1. plot(dates3,y2)
2. lines(dates2,y)

1. stop(-1)

print(ya1)

print ("ggplot svg ...")

#dev.off()

svg(filename="./suomen_koronavirustapauksia_ennuste_gompertz_1.svg", width=10, height=5, pointsize=12)

spanni=0.1 metodi="loess"

ggplot(xy4, aes(x =dates3 , y = y2))+ ggtitle(" Ennuste - koronatapausten summa") + xlab("Kuukausi") + ylab("Tapauksia")+ # scale_y_continuous(labels=NotFancy)+ scale_y_continuous(breaks= pretty_breaks(), labels = comma_format(big.mark = " ", decimal.mark = ",") ) +

theme(title=element_text(size=20), axis.text=element_text(size=18,face="bold"),axis.title=element_text(size=20,face="bold"))+ geom_line()+ geom_smooth( fill="#a0a0ff",span=spanni, method=metodi, level=0.9999, size=2)+ geom_smooth( fill="#9090ff", span=spanni,method=metodi, level=0.7) + geom_smooth( fill="#8a08af", span=spanni, method=metodi,level=0.5) +

   geom_point(y=y3) +
   geom_line(y=ya1,linetype="dotted", color="blue", size=0.5)+
   geom_line(y=ya2,linetype="dotted", color="blue", size=0.5)+
   
   geom_ribbon( aes(ymin=ya1,ymax=ya2), fill="blue", alpha=0.5) +
   geom_ribbon( aes(ymin=ya3,ymax=ya4), fill="blue", alpha=0.25) 
  
   	

dev.off()

print(".")