Site icon R-bloggers

Backcast a Time Series for Covid-19 Truths

[This article was first published on DataGeeek, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A couple of months ago, Turkey’s Health Minister announced that the positive cases showing no signs of illness were not included in the statistics. This statement made an earthquake effect in Turkey, and unfortunately, the articles about covid-19 I have wrote before came to nothing.

The reason for this statement was the pressure of the Istanbul Metropolitan Mayor. He has said that according to data released by the cemetery administration, a municipal agency, the daily number of infected deaths were nearly that two times the daily number of the death tolls explained by the ministry.

So, I decided to check the mayor’s claims. To do that, I have to do some predictions; but, not for the future, for the past. Fortunately, there is a method for this that is called Backcasting. Let’s take a vector of time series and estimate .

As you can see above, the backcasting coefficients are the same as the forecasting coefficients(). For instance, in this case, the model for new cases is ARIMA(0, 1, 2) with drift:

#Function to reverse the time series
reverse_ts <- function(y)
{
  y %>% 
    rev() %>% 
    ts(start=tsp(y)[1L], frequency=frequency(y))
}

#Function to reverse the forecast
reverse_forecast <- function(object)
{
  h <- object[["mean"]] %>% length()
  
  f <- object[["mean"]] %>% frequency()
  
  object[["x"]] <- object[["x"]] %>% reverse_ts()
  
  object[["mean"]] <- object[["mean"]] %>% rev() %>% 
    ts(end=tsp(object[["x"]])[1L]-1/f, frequency=f)
  
  
  object[["lower"]] <- object[["lower"]][h:1L,]
  object[["upper"]] <- object[["upper"]][h:1L,]
  return(object)
}

We would first reverse the time series and then make predictions and again reverse the forecast results. The data that we are going to model is the number of daily new cases and daily new deaths, between the day the health minister’s explanation was held and the day the vaccine process in Turkey has begun. We will try to predict the ten days before the date 26-11-2020.

#Creating datasets
df <- read_excel("datasource/covid-19_dataset.xlsx")
df$date <- as.Date(df$date)
#The data after the date 25-11-2020:Train set
df_after<- df[df$date > "2020-11-25",]
#The data between 15-11-2020 and 26-11-2020:Test set 
df_before <- df[ df$date > "2020-11-15" & df$date < "2020-11-26",]
#Creating dataframes for daily cases and deaths
df_cases <- bc_cases %>% data.frame()
df_deaths <- bc_deaths %>% data.frame()

#Converting the numeric row names to date object
options(digits = 9)
date <- df_cases %>% 
  rownames() %>% 
  as.numeric() %>% 
  date_decimal() %>% 
  as.Date() 

#Adding date object created above to the data frames
df_cases <- date %>% cbind(df_cases) %>% as.data.frame()
colnames(df_cases)[1] <- "date"

df_deaths <- date %>% cbind(df_deaths) %>% as.data.frame()
colnames(df_deaths)[1] <- "date"

#Convert date to numeric to use in ts function
n <- as.numeric(as.Date("2020-11-26")-as.Date("2020-01-01")) + 1

#Creating time series variables
ts_cases <- df_after$new_cases %>%
  ts(start = c(2020, n),frequency = 365 )

ts_deaths <- df_after$new_deaths %>%
  ts(start = c(2020, n),frequency = 365 )

#Backcast variables
ts_cases %>%
  reverse_ts() %>%
  auto.arima() %>%
  forecast(h=10) %>%
  reverse_forecast() -> bc_cases

ts_deaths %>%
  reverse_ts() %>%
  auto.arima() %>%
  forecast(h=10) %>%
  reverse_forecast() -> bc_deaths


It might be very useful to make a function to plot the comparison for backcast values and observed data.

#Plot function for comparison
plot_fun <- function(data,column){
  ggplot(data = data,aes(x=date,y=Point.Forecast))+
    geom_line(aes(color="blue"))+
    geom_line(data = df_before,aes(x=date,y=.data[[column]],color="red"))+
    geom_line(data = df_after,aes(x=date,y=.data[[column]],color="black"))+
    geom_ribbon(aes(ymin=Lo.95, ymax=Hi.95), linetype=2,alpha=0.1,fill="blue")+
    geom_ribbon(aes(ymin=Lo.80, ymax=Hi.80), linetype=2, alpha=0.1,fill="blue")+
    scale_color_identity(name = "Lines",
                         breaks = c("black", "red", "blue"),
                         labels = c("After", "Real", "Backcast"),
                         guide = "legend")+
    ylab(str_replace(column,"_"," "))+
    theme_light()
}


plot_fun(df_cases, "new_cases")
plot_fun(df_deaths, "new_deaths")

Conclusion

When we examine the graph, the difference in death toll seems relatively close. However, the levels of daily cases are significantly different from each other. Although this estimate only covers ten days, it suggests that there is inconsistency in the numbers given.


References

To leave a comment for the author, please follow the link and comment on their blog: DataGeeek.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.