R/fixef_sim.r

In simglm: Simulate Models Based on the Generalized Linear Model

#' Simulate categorical or factor variables
#'
#' Function that simulates factor or categorical variables.
#' Is essentially a wrapper around the sample function from base R.
#'
#' @param n A list of sample sizes.
#' @param levels Scalar indicating the number of levels for categorical or
#'   factor variable. Can also specify levels as a character vector.
#' @param var_level The level the variable should be simulated at. This can either
#'      be 1, 2, or 3 specifying a level 1, level 2, or level 3 variable
#'      respectively.
#' @param replace TRUE/FALSE indicating whether levels should be sampled with
#'   replacement. Default is TRUE.
#' @param force_equal TRUE/FALSE indicating if the sample size should be forced
#'     to be equal. Should not be used with the `replace = FALSE` argument.
#' @param ... Additional parameters passed to the sample function.
#' @export
sim_factor2 <- function(
  n,
  levels,
  var_level = 1,
  replace = TRUE,
  force_equal = FALSE,
  ...
) {
  if (force_equal) {
    if (var_level == 1) {
      cat_var <- 1
      num_equal <- sum(n[['level1']]) / length(levels)

      if (whole_number(num_equal)) {
        while (any((table(cat_var) - num_equal) != 0)) {
          cat_var <- base::sample(
            x = levels,
            size = sum(n[['level1']]),
            replace = replace,
            ...
          )
        }
      } else {
        while (any(abs(round(table(cat_var) - num_equal)) > 1)) {
          cat_var <- base::sample(
            x = levels,
            size = sum(n[['level1']]),
            replace = replace,
            ...
          )
        }
      }
    } else {
      if (var_level == 2) {
        cat_var <- 1
        num_equal <- sum(n[['level2']]) / length(levels)

        if (whole_number(num_equal)) {
          while (any((table(cat_var) - num_equal) != 0)) {
            cat_var <- base::sample(
              x = levels,
              size = sum(n[['level2']]),
              replace = replace,
              ...
            )
          }
        } else {
          while (any(abs(round(table(cat_var) - num_equal)) > 1)) {
            cat_var <- base::sample(
              x = levels,
              size = sum(n[['level2']]),
              replace = replace,
              ...
            )
          }
        }

        cat_var <- rep(cat_var, times = n[['level1']])
      } else {
        cat_var <- 1
        num_equal <- n[['level3']] / length(levels)

        if (whole_number(num_equal)) {
          while (any((table(cat_var) - num_equal) != 0)) {
            cat_var <- base::sample(
              x = levels,
              size = n[['level3']],
              replace = replace,
              ...
            )
          }
        } else {
          while (any(abs(round(table(cat_var) - num_equal)) > 1)) {
            cat_var <- base::sample(
              x = levels,
              size = n[['level3']],
              replace = replace,
              ...
            )
          }
        }
        cat_var <- rep(cat_var, times = n[['level3_total']])
      }
    }
  } else {
    if (var_level == 1) {
      cat_var <- base::sample(
        x = levels,
        size = sum(n[['level1']]),
        replace = replace,
        ...
      )
    } else {
      if (var_level == 2) {
        cat_var <- rep(
          base::sample(
            x = levels,
            size = sum(n[['level2']]),
            replace = replace,
            ...
          ),
          times = n[['level1']]
        )
      } else {
        cat_var <- rep(
          base::sample(
            x = levels,
            size = n[['level3']],
            replace = replace,
            ...
          ),
          times = n[['level3_total']]
        )
      }
    }
  }

  cat_var <- factor(cat_var, levels = levels)

  cat_var
}

#' Simulate discrete variables
#'
#' Function that simulates discrete variables.
#' Is essentially a wrapper around the sample function from base R.
#'
#' @param n A list of sample sizes.
#' @param levels Scalar indicating the number of levels for discrete variable.
#'      Can also specify levels as a character vector.
#' @param var_level The level the variable should be simulated at. This can either
#'      be 1, 2, or 3 specifying a level 1, level 2, or level 3 variable
#'      respectively.
#' @param replace TRUE/FALSE indicating whether levels should be sampled with
#'   replacement. Default is TRUE.
#' @param ... Additional parameters passed to the sample function.
#' @export
sim_ordinal2 <- function(n, levels, var_level = 1, replace = TRUE, ...) {
  if (var_level == 1) {
    ord_var <- base::sample(
      x = levels,
      size = sum(n[['level1']]),
      replace = replace,
      ...
    )
  } else {
    if (var_level == 2) {
      ord_var <- rep(
        base::sample(
          x = levels,
          size = sum(n[['level2']]),
          replace = replace,
          ...
        ),
        times = n[['level1']]
      )
    } else {
      ord_var <- rep(
        base::sample(x = levels, size = n[['level3']], replace = replace, ...),
        times = n[['level3_total']]
      )
    }
  }
  ord_var
}

#' Simulate continuous variables
#'
#' Function that simulates continuous variables. Any distribution function in
#' R is supported.
#'
#' @param n A list of sample sizes.
#' @param dist A distribution function. This argument takes a quoted
#'      R distribution function (e.g. 'rnorm'). Default is 'rnorm'.
#' @param var_level The level the variable should be simulated at. This can either
#'      be 1, 2, or 3 specifying a level 1, level 2, or level 3 variable
#'      respectively.
#' @param variance The variance for random effect simulation.
#' @param ther_sim A TRUE/FALSE flag indicating whether the error simulation
#'  function should be simulated, that is should the mean and standard deviation
#'  used for standardization be simulated.
#' @param ther_val A vector of 2 that should include the theoretical mean and
#'  standard deviation of the generating function.
#' @param ceiling A numeric value that specifies the ceiling (maximum) of an
#' attribute being generated. Defaults to NULL meaning no ceiling effect.
#' If a value is specified, any data larger than integer is rounded to
#' that ceiling value.
#' @param floor A numeric value that specifies the floor (minimum) of an
#' attribute being generated. Defaults to NULL meaning no floor effect.
#' If a value is specified, any data larger than integer is rounded to
#' that floor value.
#' @param ... Additional parameters to pass to the dist_fun argument.
#' @export
sim_continuous2 <- function(
  n,
  dist = 'rnorm',
  var_level = 1,
  variance = NULL,
  ther_sim = FALSE,
  ther_val = NULL,
  ceiling = NULL,
  floor = NULL,
  ...
) {
  if (var_level == 1) {
    cont_var <- unlist(lapply(n[['level1']], FUN = dist, ...))
  } else {
    if (var_level == 2) {
      cont_var <- rep(
        unlist(lapply(sum(n[['level2']]), FUN = dist, ...)),
        times = n[['level1']]
      )
    } else {
      cont_var <- rep(
        unlist(lapply(n[['level3']], FUN = dist, ...)),
        times = n[['level3_total']]
      )
    }
  }

  if (!is.null(variance)) {
    if (ther_sim) {
      ther_val <- do.call(dist, c(list(n = 10000000), ...))
      ther <- c(mean(ther_val), sd(ther_val))

      cont_var <- standardize(cont_var, ther[1], ther[2])
    }
    if (!is.null(ther_val)) {
      cont_var <- standardize(cont_var, ther_val[1], ther_val[2])
    }
    cont_var <- cont_var %*% chol(c(variance))
  }

  if (!is.null(ceiling)) {
    cont_var <- ifelse(cont_var > ceiling, ceiling, cont_var)
  }
  if (!is.null(floor)) {
    cont_var <- ifelse(cont_var < floor, floor, cont_var)
  }
  cont_var
}

#' Simulate knot locations
#'
#' Function that generates knot locations. An example of usefulness of this function
#' would be with generation of interrupted time series data. Another application may
#' be with simulation of piecewise linear data structures.
#'
#' @param data Data to pass to the sim_knot2 function to determine knot
#'   location.
#' @param sim_args A named list with special model formula syntax. See details and examples
#'   for more information. The named list may contain the following:
#'   \itemize{
#'     \item fixed: This is the fixed portion of the model (i.e. covariates)
#'     \item random: This is the random portion of the model (i.e. random effects)
#'     \item error: This is the error (i.e. residual term).
#'   }
#' @param data Mostly internal argument.
#' @export
simulate_knot <- function(data, sim_args) {
  do.call(
    "cbind.data.frame",
    lapply(seq_along(sim_args[['knot']]), function(ii) {
      purrr::exec(sim_knot2, !!!sim_args[['knot']][[ii]], data = data)
    })
  )
}

sim_knot2 <- function(data, variable, knot_locations, right = FALSE, ...) {
  if (!right) {
    knot_locat <- c(
      min(data[[variable]]),
      knot_locations,
      (max(data[[variable]]) + 1)
    )
  } else {
    knot_locat <- c(
      (min(data[[variable]]) - 1),
      knot_locations,
      (max(data[[variable]]))
    )
  }

  cut(data[[variable]], knot_locat, labels = FALSE, right = right, ...) - 1
}

#' Simulate Time
#'
#' This function simulates data for the time variable of longitudinal data.
#'
#' @param n Sample size of the levels.
#' @param time_levels The values the time variable should take. If NULL (default),
#'   the time values are discrete integers starting at 0 and going to n - 1.
#' @param ... Currently not used.
#'
#' @export
sim_time <- function(n, time_levels = NULL, ...) {
  if (is.null(time_levels)) {
    do.call(
      'c',
      lapply(seq_along(n[['level1']]), function(xx) {
        0:(n[['level1']][xx] - 1)
      })
    )
  } else {
    do.call(
      'c',
      lapply(seq_along(n[['level1']]), function(xx) {
        time_levels[1:n[['level1']][xx]]
      })
    )
  }
}

sim_variable <- function(
  var_type = c("continuous", "factor", "ordinal", 'time'),
  ...
) {
  var_type <- match.arg(var_type)

  switch(
    var_type,
    continuous = sim_continuous2(...),
    factor = sim_factor2(...),
    ordinal = sim_ordinal2(...),
    #knot = sim_knot2(...),
    time = sim_time(...)
  )
}

#' Tidy fixed effect formula simulation
#'
#' This function simulates the fixed portion of the model using a formula syntax.
#'
#' @param data Data simulated from other functions to pass to this function. Can pass
#'  NULL if first in simulation string.
#' @param sim_args A named list with special model formula syntax. See details and examples
#'   for more information. The named list may contain the following:
#'   \itemize{
#'     \item fixed: This is the fixed portion of the model (i.e. covariates)
#'     \item random: This is the random portion of the model (i.e. random effects)
#'     \item error: This is the error (i.e. residual term).
#'   }
#' @param ... Other arguments to pass to error simulation functions.
#' @importFrom purrr modify_if
#' @importFrom stats setNames
#'
#' @export
simulate_fixed <- function(data, sim_args, ...) {
  if (is.null(parse_formula(sim_args)[['fixed']])) {
    list_formula <- parse_formula(sim_args)
    fixed_list <- lapply(seq_along(list_formula), function(xx) {
      as.character(list_formula[[xx]][['fixed']])
    })
    if (comp_list(fixed_list)) {
      fixed_formula <- list_formula[[1]][['fixed']]
    } else {
      NULL
    }
  } else {
    fixed_formula <- parse_formula(sim_args)[['fixed']]
  }

  fixed_vars <- attr(terms(fixed_formula), "term.labels")

  if (any(grepl('^factor\\(', fixed_vars))) {
    fixed_vars <- gsub("factor\\(|\\)$", "", fixed_vars)
  }
  if (any(grepl('^ns\\(|^bs\\(', fixed_vars))) {
    fixed_vars <- gsub("ns\\(|bs\\(|\\,.+\\)$", "", fixed_vars)
  }
  if (any(grepl("^poly\\(", fixed_vars))) {
    fixed_vars <- gsub("poly\\(|\\,.+\\)", "", fixed_vars)
  }
  if (any(grepl("_post$", fixed_vars))) {
    fixed_vars <- fixed_vars[!grepl("_post$", fixed_vars)]
  }

  if (!is.null(sim_args[['propensity']])) {
    propensity_data <- simulate_propensity(sim_args = sim_args)
    prop_formula <- parse_formula(sim_args[['propensity']])
  }

  if (is.null(data)) {
    n <- sample_sizes(sim_args[['sample_size']])
    id_vars <- parse_randomeffect(parse_formula(sim_args)[['randomeffect']])[[
      'cluster_id_vars'
    ]]
    multiple_member <- parse_multiplemember(
      sim_args,
      parse_randomeffect(parse_formula(sim_args)[['randomeffect']])
    )

    id_vars <- id_vars[id_vars %ni% multiple_member[['multiple_member_idvars']]]
    ids <- create_ids(n, c('level1_id', id_vars))
    Xmat <- do.call(
      "cbind.data.frame",
      lapply(seq_along(sim_args[['fixed']]), function(ii) {
        purrr::exec(sim_variable, !!!sim_args[['fixed']][[ii]], n = n)
      })
    )
  } else {
    n <- compute_samplesize(data, sim_args)
    Xmat <- do.call(
      "cbind.data.frame",
      lapply(seq_along(sim_args[['fixed']]), function(ii) {
        purrr::exec(sim_variable, !!!sim_args[['fixed']][[ii]], n = n)
      })
    )
  }
  if (!is.null(sim_args[['knot']])) {
    num_knot_vars <- length(sim_args[['knot']])

    knot_names <- names(sim_args[['knot']])
    knot_loc <- unlist(lapply(seq_along(seq_len(num_knot_vars)), function(xx) {
      grep(knot_names[xx], fixed_vars)
    }))

    fixed_vars_knot <- fixed_vars[-knot_loc]

    if (any(grepl(":|^I", fixed_vars_knot))) {
      int_loc <- grep(":|^I", fixed_vars_knot)
      if (!is.null(sim_args[['propensity']])) {
        prop_vars <- c(
          attr(
            terms(parse_formula(sim_args[['propensity']])[['fixed']]),
            "term.labels"
          ),
          parse_formula(sim_args[['propensity']])[['outcome']]
        )
        fixed_vars <- fixed_vars[!(fixed_vars %in% prop_vars)]
      }
      colnames(Xmat) <- fixed_vars_knot[-int_loc]
    } else {
      if (!is.null(sim_args[['propensity']])) {
        prop_vars <- c(
          attr(
            terms(parse_formula(sim_args[['propensity']])[['fixed']]),
            "term.labels"
          ),
          parse_formula(sim_args[['propensity']])[['outcome']]
        )
        fixed_vars <- fixed_vars[!(fixed_vars %in% prop_vars)]
      }
      colnames(Xmat) <- fixed_vars_knot
    }

    Xmat_knot <- do.call(
      "cbind.data.frame",
      lapply(seq_along(sim_args[['knot']]), function(ii) {
        purrr::exec(sim_knot2, !!!sim_args[['knot']][[ii]], data = Xmat)
      })
    )

    Xmat <- cbind(Xmat, Xmat_knot)
  }

  if (any(grepl(":|^I", fixed_vars))) {
    int_loc <- grep(":|^I", fixed_vars)
    if (!is.null(sim_args[['propensity']])) {
      prop_vars <- c(
        attr(
          terms(parse_formula(sim_args[['propensity']])[['fixed']]),
          "term.labels"
        ),
        parse_formula(sim_args[['propensity']])[['outcome']]
      )
      fixed_vars <- fixed_vars[!(fixed_vars %in% prop_vars)]
    }
    colnames(Xmat) <- fixed_vars[-int_loc]
  } else {
    if (!is.null(sim_args[['propensity']])) {
      prop_vars <- c(
        attr(
          terms(parse_formula(sim_args[['propensity']])[['fixed']]),
          "term.labels"
        ),
        parse_formula(sim_args[['propensity']])[['outcome']]
      )
      fixed_vars <- fixed_vars[!(fixed_vars %in% prop_vars)]
    }
    colnames(Xmat) <- fixed_vars
  }

  if (!is.null(sim_args[['propensity']])) {
    if (nrow(Xmat) == 0) {
      Xmat <- data.frame(propensity_data)
    } else {
      Xmat <- cbind.data.frame(Xmat, propensity_data)
    }
  }

  # Place holder for post-process effects
  if (any(grepl("_post$", attr(terms(fixed_formula), "term.labels")))) {
    detect_ifelse <- unlist(lapply(seq_along(sim_args[['post']]), function(ii) {
      grepl('ifelse', sim_args[['post']][[ii]][['fun']])
    }))

    post_names <- names(sim_args[['post']])

    sim_args_post_ifelse <- sim_args[['post']][detect_ifelse]
    sim_args_post_other <- sim_args[['post']][!detect_ifelse]

    if (length(sim_args_post_ifelse) > 0) {
      Xmat_post_ifelse <- do.call(
        "cbind.data.frame",
        lapply(seq_along(sim_args_post_ifelse), function(ii) {
          unlist(
            lapply(
              X = eval(parse(
                text = paste0(
                  'Xmat[["',
                  sim_args_post_ifelse[[ii]][['variable']],
                  '"]]',
                  sim_args_post_ifelse[[ii]][['condition']]
                )
              )),
              FUN = sim_args_post_ifelse[[ii]][['fun']],
              yes = sim_args_post_ifelse[[ii]][['yes']],
              no = sim_args_post_ifelse[[ii]][['no']]
            )
          )
        })
      )
      names(Xmat_post_ifelse) <- names(sim_args_post_ifelse)

      Xmat <- cbind(Xmat, Xmat_post_ifelse)
    }

    if (length(sim_args_post_other) > 0) {
      Xmat_tmp <- data.frame(Xmat, ids)

      Xmat_post_other <- #do.call("cbind.data.frame",
        lapply(seq_along(sim_args_post_other), function(ii) {
          setNames(
            aggregate(
              Xmat_tmp[[sim_args_post_other[[ii]][['variable']]]] ~
                Xmat_tmp[[sim_args_post_other[[ii]][['by']]]],
              FUN = sim_args_post_other[[ii]][['fun']]
            ),
            list(
              sim_args_post_other[[ii]][['by']],
              names(sim_args_post_other)[ii]
            )
            #)
          )
        })
      #names(Xmat_post_other) <- names(sim_args_post_other)

      Xmat_tmp2 <- do.call(
        "cbind.data.frame",
        lapply(seq_along(Xmat_post_other), function(ii) {
          merge(
            Xmat_tmp,
            Xmat_post_other[[ii]],
            by = sim_args_post_other[[ii]][['by']],
            all.x = TRUE
          )
        })
      )[names(sim_args_post_other)]

      Xmat <- cbind(Xmat, Xmat_tmp2)
    }
  }

  if (
    any(
      unlist(lapply(seq_along(sim_args[['fixed']]), function(xx) {
        sim_args[['fixed']][[xx]]$var_type
      })) ==
        'factor'
    ) |
      any(grepl("^ns\\(|^bs\\(|^poly\\(", attr(terms(fixed_formula), "term.labels")))
  ) {
    if (any(grepl('^factor\\(', fixed_vars))) {
      fixed_vars <- gsub("factor\\(|\\)$", "", fixed_vars)
    }

    num_levels <- lapply(seq_along(sim_args[['fixed']]), function(xx) {
      sim_args[['fixed']][[xx]][['levels']]
    })
    num_levels <- purrr::modify_if(num_levels, is.character, length)

    if (
      any(unlist(lapply(seq_along(sim_args[['fixed']]), function(xx) {
        num_levels[[xx]] > 1 &
          sim_args[['fixed']][[xx]][['var_type']] == 'factor'
      })))
    ) {
      fixed_vars_factornames <- attr(terms(fixed_formula), "term.labels")
      if (any(grepl('^factor\\(', fixed_vars_factornames))) {
        fixed_vars_factornames <- gsub(
          "factor\\(|\\)$",
          "",
          fixed_vars_factornames
        )
      }

      fixed_vars <- factor_names(sim_args, fixed_vars_factornames)
    }
    if (any(grepl("^ns\\(|^bs\\(|^poly\\(", attr(terms(fixed_formula), "term.labels")))) {
      if (any(grepl("^ns\\(|^bs\\(|^poly\\(", fixed_vars))) {
        fixed_vars <- poly_ns_names(sim_args, fixed_vars)
      } else {
        fixed_vars <- poly_ns_names(
          sim_args,
          attr(terms(fixed_formula), "term.labels")
        )
      }
    }

    Omat <- Xmat
    fixed_formula <- attach_spline_env(fixed_formula)
    Xmat <- data.frame(model.matrix(fixed_formula, Xmat, ...))
    colnames(Xmat)[2:ncol(Xmat)] <- fixed_vars

    if (
      any(
        unlist(lapply(seq_along(sim_args[['fixed']]), function(xx) {
          sim_args[['fixed']][[xx]]$var_type
        })) ==
          'factor'
      )
    ) {
      Omat_factor <- Omat[unique_columns(Xmat, Omat)]
      Omat_factor_names <- attr(terms(fixed_formula), "term.labels")
      if (any(grepl(":|^I", Omat_factor_names))) {
        Omat_factor_names <- Omat_factor_names[-int_loc]
      }
    } else {
      Omat_factor <- NULL
    }
    if (any(grepl("^ns\\(|^bs\\(|^poly\\(", attr(terms(fixed_formula), "term.labels")))) {
      fixed_vars_new <- attr(terms(fixed_formula), "term.labels")
      fixed_vars_poly_ns <- gsub(
        "poly\\(|ns\\(|bs\\(|\\,.+\\)$",
        "",
        fixed_vars_new[grepl("^poly\\(|^ns\\(|^bs\\(", fixed_vars_new)]
      )
      if (any(grepl(":", fixed_vars_poly_ns))) {
        fixed_vars_poly_ns <- fixed_vars_poly_ns[
          !grepl(":", fixed_vars_poly_ns)
        ]
      }
      Omat_poly_ns <- Omat[, fixed_vars_poly_ns, drop = FALSE]
    } else {
      Omat_poly_ns <- NULL
    }

    Omat <- dplyr::bind_cols(Omat_factor, Omat_poly_ns)

    Xmat <- dplyr::bind_cols(Xmat, Omat)
  } else {
    fixed_formula <- attach_spline_env(fixed_formula)
    Xmat <- data.frame(model.matrix(fixed_formula, Xmat, ...))

    if (grepl("^0", as.character(fixed_formula)[2])) {
      colnames(Xmat)[seq_len(ncol(Xmat))] <- attr(
        terms(fixed_formula),
        "term.labels"
      )
    } else {
      colnames(Xmat)[2:ncol(Xmat)] <- attr(terms(fixed_formula), "term.labels")
    }
  }

  # if (!is.null(sim_args[['propensity']])) {
  #   # avoid duplicating column names
  #   keep_cols <- setdiff(names(propensity_data), names(Xmat))
  #   if (length(keep_cols) > 0L) {
  #     Xmat <- cbind(Xmat, propensity_data[keep_cols])
  #   }
  # }

  if (is.null(data)) {
    data.frame(Xmat, ids)
  } else {
    data.frame(data, Xmat)
  }
}

#' Tidy heterogeneity of variance simulation
#'
#' This function simulates heterogeneity of level one error variance.
#'
#' @param data Data simulated from other functions to pass to this function.
#' This function needs to be specified after `simulate_fixed` and `simulate_error`.
#' @param sim_args A named list with special model formula syntax. See details and examples
#'   for more information. The named list may contain the following:
#'   \itemize{
#'     \item fixed: This is the fixed portion of the model (i.e. covariates)
#'     \item random: This is the random portion of the model (i.e. random effects)
#'     \item error: This is the error (i.e. residual term).
#'   }
#' @param ... Other arguments to pass to error simulation functions.
#' @export
simulate_heterogeneity <- function(data, sim_args, ...) {
  heterogeneity_error <- heterogeneity(
    variance = sim_args[['heterogeneity']][['variance']],
    fixef = data,
    variable = sim_args[['heterogeneity']][['variable']],
    err = data[['error']]
  )

  data.frame(
    data[, !(names(data) %in% 'error')],
    error = heterogeneity_error,
    orig_error = data[['error']]
  )
}