Adding slope.estimate.quant function

Author: rafalopespx

R/nowcasting_age.R

@@ -9,7 +9,7 @@
 #' @export
 nowcasting_age <- function(dataset,
                            zero_inflated=FALSE){
-  ## Cehcl for zero-inflated
+  ## [Not in use] Check for zero-inflated
   if (zero_inflated){
     family <- "zeroinflatednbinomial1"
     control.family <- list(

R/nowcasting_inla.R

@@ -31,7 +31,8 @@
 #' @param age_col Column for ages
 #' @param date_onset Column of dates of onset of the events, normally date of onset of first symptoms of cases
 #' @param date_report Column of dates of report of the event, normally date of digitation of the notification of cases
-#' @param trajectories Returns the trajectories estimated on the inner 'INLA' model
+#' @param trajectories Returns the trajectories estimated from the inner 'INLA' model
+#' [Default] FALSE.
 #' @param zero_inflated [Experimental] In non-structured models, fit a model that deals with zero-inflated data.
 #' [Default] FALSE. If the [age_col] is not missing this flag is ignored.
 #' @param ... list parameters to other functions
@@ -82,7 +83,7 @@ nowcasting_inla <- function(dataset,
 
   }
   if(missing(date_onset) | missing(date_report)){
-    stop("date_onset or date_report missing! Please give a column name for each of this parameters")
+    stop("'date_onset' or 'date_report' missing! Please give a column name for each of this parameters")
   }
   if(K < 0 ){
     stop("K less than 0, we cannot produce backcasting! \n
@@ -113,34 +114,34 @@ nowcasting_inla <- function(dataset,
     ## Missing trim.data warning
     if(missing(trim.data)){
       trim.data <- 0
-      warning("Using default to trim dates, trim.data = 0")
+      warning("Using default to trim dates, 'trim.data = 0'")
     }else{
       trim.data<-trim.data
       message("Using trim.data inputed")
     }
     ## Missing Dmax warning
     if(missing(Dmax)){
       Dmax <- 15
-      warning("Using default to maximum delay, Dmax = 15")
+      warning("Using default to maximum delay, 'Dmax = 15'")
     }else{
       Dmax<-Dmax
       message("Using Dmax inputed")
     }
     ## Missing wdw warning
     if(missing(wdw)){
       wdw <- 30
-      warning("Using default to window of action, wdw = 30")
+      warning("Using default to window of action, 'wdw = 30'")
     }else{
       wdw<-wdw
       message("Using wdw inputed")
     }
     ## Missing data.by.week warning
     if(missing(data.by.week)){
       data.by.week <- FALSE
-      warning("Using default to returning option for the data, data.by.week = FALSE")
+      warning("Using default to returning option for the data, 'data.by.week = FALSE'")
     }else{
       data.by.week<-data.by.week
-      message("Returning data.by.week")
+      message("Returning 'data.by.week'")
     }
     # ## Missing return.age warning
     # if(missing(return.age)){
@@ -149,9 +150,9 @@ nowcasting_inla <- function(dataset,
     # }
     ## Missing age_col warning
     if(missing(age_col)){
-      warning("Age_col missing, nowcasting with unstructured model")
+      warning("'age_col' missing, nowcasting with unstructured model")
     }else{
-      message("Age col inputed, nowcasting with structured model")
+      message("'age_col' inputed, nowcasting with structured model")
     }
 
     if(missing(trajectories) | trajectories == FALSE){
@@ -162,7 +163,7 @@ nowcasting_inla <- function(dataset,
 
       if(!missing(age_col) & !missing(zero_inflated)){
         zero_inflated<-FALSE
-        warning("age_col parsed, zero_inflated ignored!")
+        warning("'age_col' parsed, 'zero_inflated' ignored!")
       }
     }
   }

R/nowcasting_no_age.R

@@ -6,11 +6,16 @@
 #' @param dataset data pre formatted in to age classes and delays by week for each cases,
 #' delay triangle format
 #' @param zero_inflated zero-inflated model
+#' [Default] FALSE.
 #'
 #' @return Trajectories from the inner 'INLA' model
 #' @export
 nowcasting_no_age <- function(dataset,
                               zero_inflated=FALSE){
+  ## Safe test
+  if(missing(dataset)){
+    stop("'dataset' is missing in 'nowcasting_no_age()'.")
+  }
 
   ## Check for the zero-inflated
   if (zero_inflated){

R/slope.estimate.quant.R

@@ -1,53 +1,74 @@
-#' @title Slope.estimate
+#' @title Slope.estimate.quant
+#'
+#' @description
+#' Fits a linear model to trajectories returned from 'nowcasting_inla()' within a given window [Default] is 3 weeks. If 'end.week' is missing uses the maximum date in 'trajectories'.
 #'
-#' @param end.week [in weeks] The end of the week wanted to the slope estimate
 #' @param trajectories Data.frame with the predicted or nowcasted estimate
+#' @param end.week [in weeks] The end of the week wanted to the slope estimate
+#' [Default] max. date in 'trajectories'.
 #' @param window [in weeks] Window of how much time will be used to calculate the slope estimate
+#' [Default] 3 weeks.
 #'
-#' @return The slope of the estimate
+#' @return The numerical value of the slope of the estimate
 #' @export
 #'
 #' @examples
+#' # Loading Belo Horizonte SARI dataset
+#' data(sragBH)
+#' now <- nowcasting_inla(dataset = sragBH,
+#'                 date_onset = DT_SIN_PRI,
+#'                 date_report = DT_DIGITA,
+#'                 trajectories = TRUE,
+#'                 silent = T)
+#' slope.estimate.quant(trajectories = now$trajectories)
 slope.estimate.quant <- function(end.week, trajectories, window=3){
 
+  if(missing(trajectories)){
+    stop("'trajectories' is missing!")
+  }
+
+  ## If 'end.week' is missing uses the max. date from trajectories
   if(missing(end.week)){
     end.week<-max(trajectories$dt_event)
+    warning("Using max. date in 'trajectories' as 'end.week'")
   }
 
+  ## Handling trajectories object
   trajectories<-trajectories |>
     dplyr::rename(Date = dt_event,
-                  Casos = Y)
+                  Cases = Y)
 
-  ## Testing if the amount of weeks in the data encompasse the size of the window
+  ## Testing if the amount of weeks in the data encompass the size of the window
   base.week <- end.week - (window*7 - 7)
   if (!base.week %in% trajectories$Date){
-    return(NA)
+    base.week <- min(trajectories$Date)
+    warning("'base.week' wasn't present in 'trajectories', using min. of date.")
   }
 
   ## Normalizing the cases column name
-  norm.casos <- trajectories |>
+  norm.cases <- trajectories |>
     dplyr::filter(Date == base.week) |>
-    dplyr::rename(valorbase = Casos) |>
-    dplyr::group_by(sample, valorbase) |>
-    dplyr::select(-Time, -Date)
-
-  norm.casos <- norm.casos |>
+    dplyr::rename(database_value = Cases) |>
+    dplyr::group_by(sample, database_value) |>
+    dplyr::select(-Time, -Date) |>
     dplyr::right_join(trajectories |>
                         dplyr::filter(Date >= base.week & Date <= end.week),
                       by='sample') |>
-    dplyr::mutate(Casos = case_when(
-      valorbase > 0 ~ Casos/valorbase,
-      TRUE ~ Casos))
+    dplyr::mutate(Cases = dplyr::case_when(
+      database_value > 0 ~ Cases/database_value,
+      TRUE ~ Cases))
 
   ## Testing the amount of samples in the trajectories, if it unique, fit the line and get the slope
   ## if it is more than one, fit for each sample the line and does statistical estimates for the slope
-  if (length(norm.casos$sample |>  unique()) == 1){
+  if (length(norm.cases$sample |>  unique()) == 1){
 
     ## The line model fitting
-    tmp <- stats::lm(Casos ~ Date, data = norm.casos)
+    tmp <- stats::lm(Cases ~ Date, data = norm.cases)
+
     ## Confidence Intervals
     l <- stats::confint(tmp, parm = 'Date', level = .9)
     q <- stats::confint(tmp, parm = 'Date', level = .5)
+
     ## Slope estimate
     slope <- dplyr::case_when(
       l[1] > 0 ~ 1,
@@ -57,9 +78,12 @@ slope.estimate.quant <- function(end.week, trajectories, window=3){
       l[2] < 0 ~ -1
     ) |>
       as.numeric()
+
   } else {
+
     ## More than one sample
-    tmp <- lme4::lmList(Casos ~ Date | sample, data = norm.casos)
+    tmp <- lme4::lmList(Cases ~ Date | sample, data = norm.cases)
+
     slope <- stats::coefficients(tmp) |>
       dplyr::select(Date) |>
       ## Quantiles calculations
@@ -78,6 +102,7 @@ slope.estimate.quant <- function(end.week, trajectories, window=3){
       as.numeric()
 
   }
+
   ## Return the slope
   return(slope)
 }