R/mix_models.R

In PatrickEslick/MixedModelLoads: Mixed Model Loads

#' Calculate loads from a model list
#'
#' @param model_list a list of modelEquation objects, see \code{modelEquation}
#' @param time_series a time series of data including the datetime and flow
#' @param datetime_heading the name of the datetime heading
#' @param flow_heading the name of the flow heading
#' @param retransformation a list of functions used to retransform the model calculations before
#' using them to calculate the load. NULL if no retransformation needed. Should be the same
#' length as model list
#' @param bias_correction_factor a vector of bias correction factors to multiply model calculations
#' by after retransforming the model, if applicable. Should be the same ling as model_list
#' @param conversion_factor the conversion factor to use
#' @param combine a logical indicating whether to combine the models, favoring
#' the models in the order of \code{model_list}
#'
#' @return a data frame of loads calculated from each model and a column mixing
#' the loads using the order of the model_list for preference.
#'
#' @export
#'

mix_models <- function(model_list, time_series,
                       datetime_heading = NULL, flow_heading = NULL,
                       retransformation = NULL, bias_correction_factor = NULL,
                       conversion_factor = 1,
                       combine = FALSE) {

  if(!is.null(datetime_heading)) {
    mixed <- data.frame(datetime = time_series[,datetime_heading])
  } else {
    mixed <- data.frame()
  }
  if(!is.null(flow_heading)) {
    flow <- time_series[,flow_heading]
  } else {
    flow <- NULL
  }
  if(length(retransformation) == 1) {
    if(is.null(retransformation)) {
      retransformation <- rep(NULL, length(model_list))
    }
  }
  for(i in seq_along(model_list)) {
    mixed[,paste0("Model", i)] <-
      apply_model(model_list[[i]], time_series,
                  conversion_factor = conversion_factor,
                  flow = flow, retransformation = retransformation[[i]],
                  bias_correction_factor = bias_correction_factor[i])
  }
  if(combine) {
    mixed$Combined <- rep(NA, nrow(mixed))
    for(i in seq_along(model_list)) {
      mixed$Combined[is.na(mixed$Combined)] <-
        mixed[is.na(mixed$Combined),paste0("Model", i)]
    }
  }
  return(mixed)

}

#' Align models
#'
#' @param mixed_models a data frame of output from \code{mixed_models}, with two models.
#' More than two models not currently supported.
#' @param datetime_heading the heading of the datetime column
#' @param datetime_format the format of the datetimes
#'
#' @return the model data frame with an additional column aligning the two models by shifting
#' the results of the second model to match up with the results of the first model
#'
#' @export
#'

align_models <- function(mixed_models, datetime_heading, datetime_format) {

  mixed_models[,datetime_heading] <- as.POSIXct(mixed_models[,datetime_heading], format = datetime_format)
  mixed_models$Difference <- mixed_models$Model1 /  mixed_models$Model2
  mixed_models$Used <- "None"
  mixed_models$Used[!is.na(mixed_models$Model1)] <- "Model1"
  mixed_models$Used[is.na(mixed_models$Model1) & !is.na(mixed_models$Model2)] <- "Model2"

  has_model <- mixed_models[mixed_models$Used != "None",]
  used2 <- which(has_model$Used == "Model2")

  #Loop through the row numbers to find which are consecutive
  gap <- 1
  gapno <- vector()
  for(i in 2:length(used2)) {
    if(used2[i] == used2[i - 1] + 1) {
      gapno[i-1] <- gap
    } else {
      gapno[i-1] <- gap
      gap <- gap + 1
    }
  }
  if(used2[length(used2)] == (used2[length(used2) - 1] + 1)) {
    gapno[length(used2)] <- gapno[length(used2) - 1]
  } else {
    gapno[length(used2)] <- gapno[length(used2) - 1] + 1
  }
  gap_df <- data.frame(used2, gapno)
  g <- start <- end <- vector()
  for(i in min(gap_df$gapno):max(gap_df$gapno)) {
    g[length(g) + 1] <- i
    start[length(start) + 1] <- min(gap_df$used2[gap_df$gapno == i])
    end[length(end) + 1] <- max(gap_df$used2[gap_df$gapno == i])
  }
  gap_info <- data.frame(Gap = g, Start = start, End = end)

  has_model$Adjustment <- 1

  for(i in 1:nrow(gap_info)) {
    start_adjustment <- has_model$Difference[gap_info$Start[i] - 1]
    start_datetime <- has_model[gap_info$Start[i] - 1, datetime_heading]
    end_adjustment <- has_model$Difference[gap_info$End[i] + 1]
    end_datetime <- has_model[gap_info$End[i] + 1, datetime_heading]
    dt <- has_model[gap_info$Start[i]:gap_info$End[i], datetime_heading]
    if(!any(is.na(c(start_adjustment, end_adjustment))) & !any(is.na(c(start_datetime, end_datetime)))) {
      adj <- approx(x = c(start_datetime, end_datetime),
                    y = c(start_adjustment, end_adjustment),
                    xout = dt)$y
    } else {
      adj <- NA
    }
    has_model$Adjustment[gap_info$Start[i]:gap_info$End[i]] <- adj
  }

  has_model$Aligned <- has_model$Combined * has_model$Adjustment
  has_model <- select(has_model, -Difference, -Used, -Adjustment)

  output <- mixed_models[,datetime_heading, drop = FALSE]
  output <- left_join(output, has_model, by = datetime_heading)
  return(output)


}