From a8001bdcc3c969b45caf30e804a4ac4c2bb69a64 Mon Sep 17 00:00:00 2001 From: Marcelo Ponce Date: Sun, 5 Jul 2020 15:44:56 -0400 Subject: [PATCH] adding plotting of time deriv and force of infection=beta*I(t) for SIR model --- DESCRIPTION | 4 ++-- NEWS | 5 ++++- R/covid19_models.R | 47 +++++++++++++++++++++++++++++++++++----------- R/covid19_plts.R | 10 ++++++---- 4 files changed, 48 insertions(+), 18 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index e1e623fbf..e214c3d41 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: covid19.analytics Type: Package Title: Load and Analyze Live Data from the CoViD-19 Pandemic -Version: 1.1 -Date: 2020-05-02 +Version: 1.1.1 +Date: 2020-07-02 Author: Marcelo Ponce [aut, cre] Maintainer: Marcelo Ponce Description: Load and analyze updated time series worldwide data of reported cases for the Novel CoronaVirus Disease (CoViD-19) from the Johns Hopkins University Center for Systems Science and Engineering (JHU CSSE) data repository . The datasets are available in two main modalities, as a time series sequences and aggregated for the last day with greater spatial resolution. Several analysis, visualization and modelling functions are available in the package that will allow the user to compute and visualize total number of cases, total number of changes and growth rate globally or for an specific geographical location, while at the same time generating models using these trends; generate interactive visualizations and generate Susceptible-Infected-Recovered (SIR) model for the disease spread. diff --git a/NEWS b/NEWS index a27a75bd5..b1d3b428d 100644 --- a/NEWS +++ b/NEWS @@ -1,8 +1,11 @@ May 2020: XXXX _ ver 1.1.1 + - added tests cases + - geneate.SIR.model()/plt.SIR.model(): added plotting of time derivatives and "force of infection" + - covid19.Toronto.data(): updated to handle new format of the data, 3 categories in columns - covid19.Toronto.data(): can return original data as reported by the city in google-docs in a list format - covid19.Toronto.data(): improved reading of Toronto data to protect possible variations in the reported data - added new argument 'interactive.display' for functions that generate interactive figures, eg. live.map(), totals.plt(), plt.SIR.model(), itrends(); for turning on/off the display of the figure - - fixed a few bugs in itrends: failing when not arguments or only one geo.locn. were specified + - fixed a few bugs in itrends: failing when not arguments or only one geo.locn, were specified - fixed data integrity functions not being exported May 2020: Several new features _ ver 1.1 diff --git a/R/covid19_models.R b/R/covid19_models.R index 2544661b9..bbeb13ce5 100644 --- a/R/covid19_models.R +++ b/R/covid19_models.R @@ -278,7 +278,7 @@ generate.SIR.model <- function(data=NULL, geo.loc="Hubei", header("=") # group fit and data in an object to be returned - SIR.model <- list(Infected=Infected, model=fit) + SIR.model <- list(Infected=Infected, model=fit, params=list(beta=Opt_par["beta"],gamma=Opt_par["gamma"],R0=R0)) # display plots if requested if (staticPlt | interactiveFig) @@ -292,7 +292,8 @@ generate.SIR.model <- function(data=NULL, geo.loc="Hubei", ####################################################################### plt.SIR.model <- function(SIR.model, geo.loc="", - interactiveFig=FALSE, fileName=NULL, interactive.display=TRUE) { + interactiveFig=FALSE, fileName=NULL, interactive.display=TRUE, + add.extras=TRUE) { #' function to plot the results from the SIR model fn #' #' @param SIR.model model resulting from the generate.SIR.model() fn @@ -300,6 +301,7 @@ plt.SIR.model <- function(SIR.model, geo.loc="", #' @param interactiveFig optional flag to activate interactive plot #' @param fileName file where to save the HTML version of the interactive figure #' @param interactive.display boolean flag to indicate whether the interactive plot will be displayed (pushed) to your browser +#' @param add.extras boolean flag to add extra indicators, such as, the "force of infection" and time derivatives #' #' @export plt.SIR.model #' @@ -325,14 +327,19 @@ plt.SIR.model <- function(SIR.model, geo.loc="", par(mfrow = c(2, 2)) plot(Day, Infected, type ="b") + if (add.extras) lines(Day[-1],diff(Infected), lty=2) plot(Day, Infected, log = "y") abline(lm(log10(Infected) ~ Day)) + if (add.extras) lines(Day[-1],log10(diff(Infected)), lty=2) + title(paste("Confirmed Cases 2019-nCoV:",toupper(geo.loc)), outer = TRUE, line = -2) col <- c("blue","red","green") matplot(fit$time, fit[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col) + if (add.extras) matlines(fit$time[-1], sapply(fit[ , 2:4],diff), type = "l", lwd = 1, lty = 2, col = col) matplot(fit$time, fit[ , 2:4], type = "l", xlab = "Day", ylab = "Number of subjects", lwd = 2, lty = 1, col = col, log = "y") + if (add.extras) matlines(fit$time[-1], sapply(fit[ , 2:4],diff), type = "l", lwd = 1, lty = 2, col = col, log="y") ## Warning in xy.coords(x, y, xlabel, ylabel, log = log): 1 y value <= 0 ## omitted from logarithmic plot @@ -341,6 +348,7 @@ plt.SIR.model <- function(SIR.model, geo.loc="", title(paste("SIR model 2019-nCoV:", toupper(geo.loc)), outer = TRUE, line = -22) #axis.Date(1,as.Date(names(data)[colOffset:ncol(data)])) + ### INTERACTIVE PLOTS if (interactiveFig) { # load/check plotly @@ -351,28 +359,45 @@ plt.SIR.model <- function(SIR.model, geo.loc="", length(Infected) <- length(fit[,1]) - loc.data <- cbind(Infected,fit[,1:4]) + loc.data <- cbind(Infected,fit[,1:4], Force=fit[,3]*SIR.model$params$beta) #print(loc.data) #model.ifig <- model.ifig %>% add_trace(y = ~Infected, name="Actual data", type='scatter', mode='markers', visible=TRUE) # add traces - model.ifig <- add.N.traces(model.ifig,loc.data, c("data","Susceptible", "Infected", "Recovered"), vis=TRUE) + model.ifig <- add.N.traces(model.ifig,loc.data, c("data","Susceptible", "Infected", "Recovered","Force"), vis=TRUE) # extra traces for activating log-scale - model.ifig <- add.N.traces(model.ifig, log10(loc.data), c("data","Susceptible", "Infected", "Recovered"), vis=FALSE) - + model.ifig <- add.N.traces(model.ifig, log10(loc.data), c("data","Susceptible", "Infected", "Recovered", "Force"), vis=FALSE) + # log-scale menu based on nbr of traces... - updatemenues <- log.sc.setup(4) + updatemenues <- log.sc.setup(5) # add a menu for switching log/linear scale model.ifig <- model.ifig %>% layout(updatemenus=updatemenues) + ### DERIVS ### + derivs <- cbind(time=fit$time[-1],as.data.frame(lapply(fit[,2:4],diff)) ,force=diff(loc.data$Force)) + #print(str(derivs)) + derivs.ifig <- plot_ly(data=derivs, x=time) + derivs.data <- cbind(diff(Infected),derivs) + #print(str(derivs.data)) + derivs.ifig <- add.N.traces(derivs.ifig, derivs.data, c("Deriv.data","Deriv.Susceptible", "Deriv.Infected", "Deriv.Recovered","Deriv.Force"), vis=TRUE) + derivs.ifig <- add.N.traces(derivs.ifig, log10(derivs.data), c("Deriv.data","Deriv.Susceptible", "Deriv.Infected", "Deriv.Recovered","Deriv.Force"), vis=FALSE) + updatemenues <- log.sc.setup(5) + derivs.ifig <- derivs.ifig %>% layout(updatemenus=updatemenues) + ### + ifigs <- subplot(list(model.ifig, derivs.ifig), nrows = 2, shareX = TRUE, titleX = TRUE) + #print(ifigs) + + + # activate interactive figure + if (interactive.display) { + print(model.ifig) + print(ifigs) + } - # activate interactive figure - if (interactive.display) print(model.ifig) - -return(model.ifig) + return(model.ifig) if (!is.null(fileName)) { FileName <- paste0(fileName,".html") diff --git a/R/covid19_plts.R b/R/covid19_plts.R index 71a3968a0..b941d3ed9 100644 --- a/R/covid19_plts.R +++ b/R/covid19_plts.R @@ -68,10 +68,12 @@ add.N.traces <- function(i.fig, traces, tr.names=rep("",length(traces)), vis=TRUE) { tr.x <- traces[,2] - i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,3], name=tr.names[2], type='scatter', mode='lines+markers', visible=vis) - i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,4], name=tr.names[3], type='scatter', mode='lines+markers', visible=vis) - i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,5], name=tr.names[4], type='scatter', mode='lines+markers', visible=vis) - i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,1], name=tr.names[1], type='scatter', mode='markers', visible=vis) + i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,3], name=tr.names[3-1], type='scatter', mode='lines+markers', visible=vis) + i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,4], name=tr.names[4-1], type='scatter', mode='lines+markers', visible=vis) + i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,5], name=tr.names[5-1], type='scatter', mode='lines+markers', visible=vis) + i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,1], name=tr.names[2-1], type='scatter', mode='markers', visible=vis) + # for adding the FORCE... + i.fig <- i.fig %>% add_trace(x=~tr.x, y = ~traces[,6], name=tr.names[6-1], type='scatter', mode='lines', visible=vis) return(i.fig) }