covid_ita.R

  # install.packages("xts","zoo","growthrates","rio","R0")
  library(xts)
  library(zoo)
  library(growthrates)
  library(rio)
  library(R0)
  
  ###############################################################################################################################################################################
  # retrieve data
  
  ## covid_ita_df=read.csv("https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv")
  
  covid_ita_df=rio::import("https://github.com/pcm-dpc/COVID-19/raw/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv")
  covid_ita_reg_df=rio::import("https://github.com/pcm-dpc/COVID-19/raw/master/dati-regioni/dpc-covid19-ita-regioni.csv")
  covid_ita_pro_df=rio::import("https://github.com/pcm-dpc/COVID-19/raw/master/dati-province/dpc-covid19-ita-province.csv")
  dates=as.Date(covid_ita_df$data)
  datej=as.numeric(format(dates,"%j"))
  
  Ratio_casi_tamponi=covid_ita_df$totale_casi/covid_ita_df$tamponi
  plot(xts(Ratio_casi_tamponi,order.by=as.Date(dates)),main="Ratio casi totali vs tamponi")                    
  
  
  covid_ita_reg_ls_regioni=split(covid_ita_reg_df,covid_ita_reg_df$codice_regione)
  covid_ita_reg_ls_province=split(covid_ita_pro_df,covid_ita_pro_df$codice_provincia)
  
  
  
  
  
  
  
  
  covid_ita_reg_df_Toscana=covid_ita_reg_df[covid_ita_reg_df$denominazione_regione=="Toscana",]
  covid_ita_reg_df_Lomb=covid_ita_reg_df[covid_ita_reg_df$denominazione_regione=="Lombardia",]
  day_biased= 34
  bias_34=115
  
  covid_ita_reg_df_Toscana$totale_casi[(day_biased-8):(day_biased-1)]=covid_ita_reg_df_Toscana$totale_casi[(day_biased-8):(day_biased-1)]-16.4285
                          
                           
    #covid_ita_reg_df_Toscana$totale_casi[1]=1
  
  
  #############################################################################################################################
  ## The reproductive number R0 of COVID-19 based on estimate of a statistical time delay dynamical system
  ## Nian Shao, Jin Cheng, Wenbin Chen
  ## doi: https://doi.org/10.1101/2020.02.17.20023747 
  # n=5,alpha=2/3
  
  mGT <- generation.time("gamma", c(7.5,11.25))
  EG <- est.R0.EG(covid_ita_df$totale_casi, mGT)
  EGtosc <- est.R0.EG(covid_ita_reg_df_Toscana$totale_casi, mGT)
  EGlomb <- est.R0.EG(covid_ita_reg_df_Lomb$totale_casi, mGT)
  SB <- est.R0.SB(covid_ita_df$totale_casi, mGT)
  SBtosc <- est.R0.SB(covid_ita_reg_df_Toscana$totale_casi, mGT)
  SBlomb <- est.R0.SB(covid_ita_reg_df_Lomb$totale_casi, mGT)
  
  # res <- get_R(covid_ita_reg_df_Toscana$totale_casi,si_mean=7.5,si_sd=11.5)
  # plot(res,"lambdas")
  
  kfit <- fit_easylinear(datej, covid_ita_df$totale_casi)
  kfittosc <- fit_easylinear(datej[-1], covid_ita_reg_df_Toscana$totale_casi[-1])
  kfitlomb <- fit_easylinear(datej, covid_ita_reg_df_Lomb$totale_casi)
  
  ###############################################################################################################################################################################
  # Italia
  
  covid_ita_cases <- data.frame(growth_TC=covid_ita_df$totale_casi, days=datej)
  covid.ss.ita <- nls(growth_TC ~ SSlogis(datej, phi1, phi2, phi3), data =covid_ita_cases)
  alpha <- coef(covid.ss.ita)  #extracting coefficients
  
  plot(growth_TC ~ days, data =covid_ita_cases, main = "Covid 19 total cases\n  in Italy since 2020-02-24 \n Italian Civil Protection", 
       xlab = "Day of Year", ylab = "T cases", xlim = c(55, 140), ylim = c(0, alpha[1]+2000))  # 
  curve(SSlogis(x, alpha[1],alpha[2],alpha[3]), add = T, col = "red") 
  #curve(alpha[1]/(1 + exp(-(x - alpha[2])/alpha[3])), add = T, col = "red")  # Fitted model
  abline(h =alpha[1], col="green")
  alpha.ita <- coef(covid.ss.ita)  #extracting coefficients
  
  K_ita=alpha.ita[1]
  r_ita=1/alpha.ita[3]
  doubletime_ita=log(2)/r_ita
  datepeakn=match(1,ifelse(SSlogis(1:130,alpha[1],alpha[2],alpha[3])-as.numeric(K_ita)>-20,1,0))
  datepeak=seq.Date(as.Date("2020-01-01"),as.Date("2020-12-31"),1)[datepeakn]
  points(tail(datej)[6],tail(covid_ita_df$totale_casi)[6],col="red")
  abline(h =alpha[1], col="green")
  text(datepeakn,100,datepeak)
  abline(v =datepeakn, col="green")
  prev=SSlogis(tail(datej)[6]+1, alpha[1],alpha[2],alpha[3])-tail(covid_ita_df$totale_casi)[6]
  text(tail(datej)[6]-10,alpha[1]-3000,paste("Per ",as.Date(tail(covid_ita_df$data)[6])+1 ," +", floor(prev)," nuovi casi"),col="red")
  
  
  fitl <- fit_easylinear(datej, covid_ita_df$totale_casi)
  p     <- c(y0 = as.numeric(coef(fitl)[2]), mumax =as.numeric(coef(fitl)[3]),K = as.numeric(coef(fitl)[4]))
  fitg <- fit_growthmodel(FUN = grow_logistic,p=p, datej, covid_ita_df$totale_casi)
  
  ###############################################################################################################################################################################
  # Toscana
  
  covid_tosc_cases <- data.frame(growth_TC=covid_ita_reg_df_Toscana$totale_casi, days=datej)
  #covid_tosc_cases=covid_tosc_cases[1:25,]
  covid.ss <- nls(growth_TC ~ SSlogis(days, phi1, phi2, phi3), data =covid_tosc_cases)
  alpha <- coef(covid.ss)  #extracting coefficients
  
  plot(growth_TC ~ days, data =covid_tosc_cases, main = "Covid 19 total cases\n  in Tuscany since 2020-02-24 \n Italian Civil Protection", 
       xlab = "Day of Year", ylab = "T cases", xlim = c(55, 130), ylim = c(0, alpha[1]+2000))  # 
  curve(SSlogis(x, alpha[1],alpha[2],alpha[3]), add = T, col = "red") 
  #curve(alpha[1]/(1 + exp(-(x - alpha[2])/alpha[3])), add = T, col = "red")  # Fitted model
  abline(h =alpha[1], col="green")
  
  alpha.tosc <- coef(covid.ss)  #extracting coefficients
  K_tosc=alpha.tosc[1]
  r_tosc=1/alpha.tosc[3]
  doubletime_tosc=log(2)/r_tosc
  datepeakn=match(1,ifelse(SSlogis(1:130,alpha[1],alpha[2],alpha[3])-as.numeric(K_tosc)>-20,1,0))
  datepeak=seq.Date(as.Date("2020-01-01"),as.Date("2020-12-31"),1)[datepeakn]
  points(tail(datej)[6],tail(covid_ita_reg_df_Toscana$totale_casi)[6],col="red")
  abline(h =alpha[1], col="green")
  text(datepeakn,100,datepeak)
  abline(v =datepeakn, col="green")
  abline(v =tail(datej,1), col="orange")
  text(tail(datej,1),100,tail(dates,1))
  
  #text(tail(datej)[6]-10,alpha[1]-300,paste("Per ",as.Date(tail(covid_ita_reg_df_Toscana$data)[6])+1 ," +", floor(prev)," nuovi casi"),col="red")
  
  plot(diff(covid_tosc_cases$growth_TC)~covid_tosc_cases$days[2:nrow(covid_tosc_cases)],
       xlab = "Day of Year", ylab = "T cases", xlim = c(55, 130),
       main = "Covid 19 daily variation in total cases\n  in Tuscany since 2020-02-24 \n Italian Civil Protection")
  lines(c(56:130),diff(SSlogis(c(55:130), alpha[1],alpha[2],alpha[3])), col = "orange")
  text(datepeakn,100,datepeak)
  abline(v =datepeakn, col="green")
  abline(v =tail(datej,1), col="orange")
  text(tail(datej,1),100,tail(dates,1))
  
  
  Ratio_casi_tamponi=covid_ita_reg_df_Toscana$totale_casi/covid_ita_reg_df_Toscana$tamponi
  plot(xts(Ratio_casi_tamponi,order.by=as.Date(dates)),main="Ratio casi totali vs tamponi Toscana")                    
  Tamponi=c(0,diff(covid_ita_reg_df_Toscana$tamponi))
  plot(xts(Tamponi,order.by=as.Date(dates)),main="Tamponi giornalieri in Toscana")                    
  
  
  fitl <- fit_easylinear(datej, covid_ita_reg_df_Toscana$totale_casi)
  p     <- c(y0 = as.numeric(coef(fitl)[2]), mumax =as.numeric(coef(fitl)[3]),K = as.numeric(coef(fitl)[4]))
    fitg <- fit_growthmodel(FUN = grow_logistic,p=p, datej, covid_ita_reg_df_Toscana$totale_casi)
  
  ###############################################################################################################################################################################
  # Lombardia
  
  covid_lomb_cases <- data.frame(growth_TC=covid_ita_reg_df_Lomb$totale_casi, days=datej)
  #covid_lomb_cases=covid_lomb_cases[1:25,]
  covid.ss <- nls(growth_TC ~ SSlogis(days, phi1, phi2, phi3), data =covid_lomb_cases)
  alpha <- coef(covid.ss)  #extracting coefficients
  
  plot(growth_TC ~ days, data =covid_lomb_cases, main = "Covid 19 total cases\n  in Lombardy since 2020-02-24 \n Italian Civil Protection", 
       xlab = "Day of Year", ylab = "T cases", xlim = c(55, 130), ylim = c(0, alpha[1]+2000))  # 
  curve(SSlogis(x, alpha[1],alpha[2],alpha[3]), add = T, col = "red") 
  #curve(alpha[1]/(1 + exp(-(x - alpha[2])/alpha[3])), add = T, col = "red")  # Fitted model
  abline(h =alpha[1], col="green")
  
  alpha.lomb <- coef(covid.ss)  #extracting coefficients
  K_lomb=alpha.lomb[1]
  r_lomb=1/alpha.lomb[3]
  doubletime_lomb=log(2)/r_lomb
  datepeakn=match(1,ifelse(SSlogis(1:130,alpha[1],alpha[2],alpha[3])-as.numeric(K_lomb)>-20,1,0))
  datepeak=seq.Date(as.Date("2020-01-01"),as.Date("2020-12-31"),1)[datepeakn]
  abline(h =alpha[1], col="green")
  text(datepeakn,100,datepeak)
  abline(v =datepeakn, col="green")
  points(tail(datej)[6],tail(covid_ita_reg_df_Lomb$totale_casi)[6],col="red")
  prev=SSlogis(tail(datej)[6]+1, alpha[1],alpha[2],alpha[3])-tail(covid_ita_reg_df_Lomb$totale_casi)[6]
  text(tail(datej)[6]-10,alpha[1]-1000,paste("Per ",as.Date(tail(covid_ita_reg_df_Lomb$data)[6])+1 ," +", floor(prev)," nuovi casi"),col="red")
  
  fitl <- fit_easylinear(datej, covid_ita_reg_df_Lomb$totale_casi)
  p    <- c(y0 = as.numeric(coef(fitl)[2]), mumax =as.numeric(coef(fitl)[3]),K = as.numeric(coef(fitl)[4]))
  fitg <- fit_growthmodel(FUN = grow_logistic,p=p, datej, covid_ita_reg_df_Lomb$totale_casi)
  
  ###############################################################################################################
  covid_ita_reg_df_sard=covid_ita_reg_df[covid_ita_reg_df$denominazione_regione=="Sardegna",]
  
  covid_sard_cases <- data.frame(growth_TC=covid_ita_reg_df_sard$totale_casi, days=datej)
  covid.ss <- nls(growth_TC ~ SSlogis(days, phi1, phi2, phi3), data =covid_sard_cases)
  alpha <- coef(covid.ss)  #extracting coefficients
  
  plot(growth_TC ~ days, data =covid_sard_cases, main = "Covid 19 total cases\n  in Sardinia since 2020-02-24 \n Italian Civil Protection", 
       xlab = "Day of Year", ylab = "T cases", xlim = c(55, 130), ylim = c(0, alpha[1]+alpha[1]*0.3))  # 
  curve(SSlogis(x, alpha[1],alpha[2],alpha[3]), add = T, col = "red") 
  #curve(alpha[1]/(1 + exp(-(x - alpha[2])/alpha[3])), add = T, col = "red")  # Fitted model
  abline(h =alpha[1], col="green")
  
  alpha.sard<- coef(covid.ss)  #extracting coefficients
  K_lomb=alpha.sard[1]
  r_sard=1/alpha.sard[3]
  doubletime_sard=log(2)/r_sard
  datepeakn=match(1,ifelse(SSlogis(1:130,alpha[1],alpha[2],alpha[3])-as.numeric(K_lomb)>-20,1,0))
  datepeak=seq.Date(as.Date("2020-01-01"),as.Date("2020-12-31"),1)[datepeakn]
  abline(h =alpha[1], col="green")
  text(datepeakn,100,datepeak)
  abline(v =datepeakn, col="green")
  points(tail(datej)[6],tail(covid_ita_reg_df_sard$totale_casi)[6],col="red")
  prev=SSlogis(tail(datej)[6]+1, alpha[1],alpha[2],alpha[3])-tail(covid_ita_reg_df_sard$totale_casi)[6]
  text(tail(datej)[6]-10,alpha[1]-alpha[1]*0.2,paste("Per ",as.Date(tail(covid_ita_reg_df_sard$data)[6])+1 ," +", floor(prev)," nuovi casi"),col="red")
  
  fitl <- fit_easylinear(datej, covid_ita_reg_df_sard$totale_casi)
  
  covid.ss <- nls(growth_TC ~ SSlogis(days, phi1, phi2, phi3), data =covid_sard_cases)
  alpha <- coef(covid.ss)  #extracting coefficients
  
  
  
  
  
  
  ###############################################################################################################
  
  # plot of tampons
  covid <- xts(x = covid_ita_df[,5:12], order.by = as.Date(covid_ita_df$data))
  plot(covid[,"Tamponi effettuati"])
  
  
  
  
  
  
  
  
  ###############################################################################################################
  # reference
  
  # https://www.statforbiology.com/nonlinearregression/usefulequations