Function that returns the last estimated rt from a time series of cummulative number of COVID-19 cases

get_last_covid_rt.r

#based on this method: https://github.com/k-sys/covid-19/blob/master/Realtime%20R0.ipynb
# In fact it is a translation from python to R of some parts of the code explained in the link above.

get_last_covid_rt<- function(df_days_case){ #this dataset must have a structre with a collumn no_cases (accumulated number of cases on a given sequence day)
  
  
  #smooth the data
  
  x <- 1:NROW(df_days_case)
  y<-df_days_case$no_cases 
  lo <- loess(y~x)
  #plot(x,y)
  #lines(predict(lo), col='red', lwd=2)
  
  z= predict(lo)
  
  k= diff(z)
  
  
  k<- k[k>=1]
  
  
  gama_series<- 1/4
  
  R_T_MAX <- 12

  rt_range<- seq(0,R_T_MAX,0.1)
  
  
  #Calculate the priors using poisson distribution, focus on the likelihood for a set of lambdas

  prior<-
    map_dfr(2:NROW(k),function(i){
      
      
      lambda<-k[i-1]*exp(gama_series *(rt_range -1))
      #likelihood_r_t <- (lambda^k[i]*exp(-lambda))/factorial(k[i])
      likelihood_r_t <- dpois(x = round(k[i]), lambda = lambda)
      #likelihood_r_t <-  likelihood_r_t/sum(likelihood_r_t)
      data.frame(day = i, k= rep(k[i],length(lambda)), likelihood_r_t= likelihood_r_t, rt_range = rt_range)
      
      
      
    })
  
  

  spread_prior<-
    prior%>%
    spread(key = rt_range, value = likelihood_r_t)
  
  
  
  cum_prior<-  
    cumprod(spread_prior[-c(1,2)])
  
  
  cum_prior$dia <- spread_prior$day
  
  cum_prior<-
    cum_prior%>%
    gather(key = "rt", value= "valor", -dia)
  
  
  posteriors<-
    cum_prior %>%
    group_by(dia)%>%
    summarise(
      soma_like_rt= sum(valor)
    ) %>% 
    inner_join(cum_prior) %>%
    mutate(likelihood_r_t = valor/soma_like_rt)
  
  

  most_likely_values <-
    posteriors%>%
    group_by(dia) %>%
    summarise(
      maximo=max(likelihood_r_t)
    ) %>%
    inner_join(posteriors) %>%
    filter(likelihood_r_t==maximo ) %>%
    select(dia, rt, maximo)
  
  #Choose the las most likely value
  
  R0_choice<-
    as.numeric((most_likely_values%>%
                  arrange(dia) %>%
                  top_n(1) %>%
                  select(rt))$rt)
  
  R0_choice
  
  
}


#using the function
download.file("https://usafactsstatic.blob.core.windows.net/public/data/covid-19/covid_confirmed_usafacts.csv" ,destfile = "covid_confirmed_usafacts.csv", mode='w')


library(readr)

covid_confirmed_usafacts <- read_csv("covid_confirmed_usafacts.csv")


covid_confirmed_usafacts$county <- paste(covid_confirmed_usafacts$`County Name`,covid_confirmed_usafacts$State, sep = "-")

data_county<-
  covid_confirmed_usafacts[-1] %>%
  filter(	county== "Dallas County-TX") %>%
  gather(key = "date", value= "no_cases", -c(1:3,(NCOL(covid_confirmed_usafacts)-1)) ) %>%
  filter(no_cases>0)



get_last_covid_rt(data_county)