adding a new fn: estimateRRs -- for computing rolling estimates for f…

…atality and recovery rates

NAMESPACE

@@ -14,17 +14,20 @@ export(covid19.US.data)
 export(covid19.data)
 export(covid19.genomic.data)
 export(data.checks)
+export(estimateRRs)
 export(generate.SIR.model)
 export(geographicalRegions)
 export(growth.rate)
 export(integrity.check)
 export(itrends)
 export(live.map)
+export(mrollingRates)
 export(mtrends)
 export(nullify.data)
 export(plt.SIR.model)
 export(preProcessingData)
 export(report.summary)
+export(rollingRate)
 export(single.trend)
 export(sweep.SIR.models)
 export(totals.plt)

NEWS

@@ -1,8 +1,9 @@
 Aug 2020:  Several new features added _ ver 1.2
 		- New genomic functions supporting the retrieval of nucleotides, proteins, FASTA and phylogenetic tree from NCBI.
 		- Adding a new function, sweep.SIR.models for generating a series of models based on a range for the starting date t0
-		- Added more checks for intregrity/consistency checks: critical columns, quantifiable columns, etc.
-		- Added a nullify.data fn for removing inconsistent/incomplete entries in the data
+		- Added more checks for integrity/consistency checks: critical columns, quantifiable columns, etc.
+		- Added a nullify.data fn, for removing inconsistent/incomplete entries in the data
+		- Added a new estimateRRs fn, for computing rolling estimates for fatality and recovery rates
 
 Jul 2020:  Several new features and bugs fixes _ ver 1.1.1
 		- added new geographical categories to select data per continents, implemented via geographicalRegions() fn

R/covid19_rollingEstimates.R

@@ -0,0 +1,146 @@
+# module for computing rolling estimates for fatality and recovery rates
+# part of covid19.analytics package
+#
+# M.Ponce
+
+#########
+
+rollingRate <- function(data, fn=mean, period=NULL) {
+#' function to compute a rolling fn of a TS data
+#'
+#' @param  data  TS data
+#' @param  fn  function to compute rolling
+#' @param  period  length of window
+#'
+#' @return a vector with rolling values
+#'
+#' @export
+#'
+
+	ini.col <- 5
+	lst.col <- ncol(data)
+
+	X <- unlist(data[ini.col:lst.col])
+
+	if (is.null(period))
+		period <- length(X)
+
+	Y <- movingFn(X,fn=fn,period=period)
+
+	names(Y) <- names(data[ini.col:lst.col])
+
+	return(Y)
+}
+
+#########
+
+mrollingRates <- function(data=NULL, geo.loc=NULL, fn=mean, period) {
+#' function to compute a rolling fn (rate) of multiple quantities from TS data, eg. fatality and recovery rate
+#'
+#' @param  data  time series dataset to consider
+#' @param  geo.loc  country/region to analyze
+#' @param  fn  function to compute rolling
+#' @param  period  length of window
+#'
+#' @export
+#'
+
+        if (is.null(data) | is.null(geo.loc))
+                stop("Arguments needed!")
+
+        # check that the data is time series data
+        chk.TS.data(data,xtp=TRUE)
+
+        # check on the location
+        geo.loc <- checkGeoLoc(data,geo.loc)
+
+        cases.per.loc <- select.per.loc(data,geo.loc)
+
+        if ("status" %in% names(cases.per.loc)) {
+                sts <- unique(cases.per.loc$status)
+		rolls <- data.frame()
+		for (j in sts) {
+			rolls <- rbind(rolls, rollingRate(cases.per.loc[cases.per.loc$status==j,1:ncol(cases.per.loc)-1], fn=fn, period=period))
+		}
+		rownames(rolls) <- sts
+		names(rolls) <- names(cases.per.loc)[5:(ncol(cases.per.loc)-1)]
+	} else {
+		rolls <- rollingRate(cases.per.loc, fn=fn, period=period)
+	}
+
+	return(rolls)
+}
+
+#########
+
+estimateRRs <- function(data=NULL, geo.loc=NULL, period=NULL, graphics.ON=TRUE, splitG=TRUE) {
+#' estimate rolling rates for a given geographical location for an specific TS data
+#'
+#' @param  data  time series dataset to consider
+#' @param  geo.loc  country/region to analyze
+#' @param  period  length of window
+#' @param  graphics.ON  boolean flag to activate/deactivate graphical output
+#' @param  splitG  boolean flag for having the graphical output separated or not
+#'
+#' @export
+#'
+#' @examples
+#' estimateRRs(covid19.data("TS-all"), geo.loc='Peru', period=7)
+#' estimateRRs(covid19.data("TS-all"),
+#'		geo.loc=c('Peru','Argentina','Uruguay','US','Spain','Japan'), period=7)
+#'
+	if (length(geo.loc) > 1) {
+		rtn <- list()
+		for (locn in geo.loc) {
+			rtn <- list(rtn,c(locn,estimateRRs(data, geo.loc=locn, period=period, graphics.ON=graphics.ON, splitG=splitG)))
+		}
+		#rtn <- rtn[[1]][[-1]]
+	} else {
+
+	rrs.m <- mrollingRates(data, geo.loc, fn=mean, period=period)
+	fat.rate <- as.numeric(rrs.m[2,]/rrs.m[1,])*100
+	rec.rate <- as.numeric(rrs.m[3,]/rrs.m[1,])*100
+	cnst.m <- mrollingRates(data, geo.loc, fn=mean, period=NULL)
+	cnst.fat.rate <- as.numeric(cnst.m[2,]/cnst.m[1,])*100
+	cnst.rec.rate <- as.numeric(cnst.m[3,]/cnst.m[1,])*100
+
+	rrs.sd <- mrollingRates(data, geo.loc, fn=sd, period=period)
+	fat.sd <- as.numeric(rrs.sd[2,]/rrs.sd[1,])
+	rec.sd <- as.numeric(rrs.sd[3,]/rrs.sd[1,])
+
+#print(rrs.m)
+
+	if (graphics.ON) {
+        ### preserve user graphical env.
+        # save the state of par() before running the code
+        oldpar <- par(no.readonly = TRUE)
+        # restore the previous state after the fn is done, even if it fails, so the user environment is not altered
+        on.exit(par(oldpar))
+        #########
+
+		if (splitG) par(mfrow=c(2,1))
+		minY <- min(c(fat.rate,rec.rate), na.rm=TRUE)
+		maxY <- max(c(fat.rate,rec.rate), na.rm=TRUE)
+		xrange <- 1:length(rrs.m)
+		plot(xrange,fat.rate, type='b', pch=20, col='red', axes=FALSE)	#, ylim=c(minY,maxY))
+		axis.Date(1,names(rrs.m)); axis(2)
+		confBand(xrange,fat.rate, 1,length(fat.rate),0,max(fat.rate,na.rm=TRUE), windowsNbr=10, lcolour='red')
+		lines(xrange,cnst.fat.rate, lty='dashed', col='red')
+
+		title(main=geo.loc)
+		if (!splitG) par(new=TRUE)
+		plot(xrange,rec.rate, type='b', pch=20, col='blue', axes=FALSE)	#, ylim=c(minY,maxY))
+		axis.Date(1,names(rrs.m)); axis(4)
+		confBand(xrange,rec.rate, 1,length(fat.rate),0,max(fat.rate,na.rm=TRUE), windowsNbr=10, lcolour='blue')
+		lines(xrange,cnst.rec.rate, lty='dashed', col='blue')
+	}
+
+	rtn <- list(RollingRateEstimates=rrs.m,
+			FatalityRate=fat.rate,
+			RecoveryRate=rec.rate)
+	}
+
+	return(rtn)
+}
+
+#########