R/util.r

In simglm: Simulate Models Based on the Generalized Linear Model

"%ni%" <- Negate("%in%")

`%||%` <- function(x, y) if (is.null(x)) y else x


is_odd <- function(x) x %% 2 != 0

whole_number <- function(x) x %% 1 == 0

comp_list <- function(x) length(unique.default(x)) == 1L

prop_limits <- function(prop) {
  if (prop > 0.5) {
    u_diff <- 1 - prop
    u_max <- 1
    u_min <- prop - u_diff
  } else {
    u_diff <- 0 + prop
    u_max <- prop + u_diff
    u_min <- 0
  }
  return(c(u_min, u_max))
}

standardize <- function(x, mean, sd) {
  new <- (x - mean) / sd
  return(new)
}

search_factors <- function(x) {
  x <- strsplit(x, ':')
  loc <- lapply(seq_along(x), function(xx) {
    grep("\\.f$|\\.c$|_f$|_c$|\\.k$|_k$", x[[xx]], ignore.case = TRUE)
  })
  len <- lapply(x, length)
  for (i in seq_along(x)) {
    x[[i]][loc[[i]]] <- paste0('factor(', x[[i]][loc[[i]]], ')')
    if (len[[i]] > 1) {
      x[[i]] <- paste(x[[i]], collapse = ':')
    }
  }
  x <- as.formula(paste('~', paste(x, collapse = '+')))
  x
}

sample_sizes <- function(sample_size) {
  if (length(sample_size) == 1) {
    sample_size <- list(level1 = sample_size)
  } else {
    if (length(sample_size) == 3) {
      if (length(sample_size[['level2']]) == sample_size[['level3']]) {
        level2 <- sample_size[['level2']]
      } else {
        level2 <- rep(sample_size[['level2']], sample_size[['level3']])
      }
      sample_size['level2'] <- list(level2)
    }
    total_level2_samplesize <- sum(sample_size[['level2']])

    if (length(sample_size[['level1']]) == total_level2_samplesize) {
      level1 <- sample_size[['level1']]
    } else {
      level1 <- rep(sample_size[['level1']], total_level2_samplesize)
    }
    sample_size['level1'] <- list(level1)
    if (length(sample_size) == 3) {
      sample_size['level3_total'] <- list(sample_size_level3(sample_size))
    }
  }
  sample_size
}

sample_size_level3 <- function(sample_size) {
  end <- cumsum(sample_size[['level2']])
  beg <- c(1, cumsum(sample_size[['level2']]) + 1)
  beg <- beg[-length(beg)]

  lvl3ss <- unlist(lapply(
    lapply(seq_along(beg), function(xx) {
      sample_size[['level1']][beg[xx]:end[xx]]
    }),
    sum
  ))

  lvl3ss
}

create_ids <- function(sample_size_list, id_names) {
  if (length(id_names) == 3) {
    id_vars <- data.frame(
      unlist(lapply(seq_along(sample_size_list[['level1']]), function(xx) {
        1:sample_size_list[['level1']][xx]
      })),
      rep(
        1:sum(sample_size_list[['level2']]),
        times = sample_size_list[['level1']]
      ),
      rep(
        1:sample_size_list[['level3']],
        times = sample_size_list[['level3_total']]
      )
    )
  } else {
    if (length(id_names) == 2) {
      id_vars <- data.frame(
        unlist(lapply(seq_along(sample_size_list[['level1']]), function(xx) {
          1:sample_size_list[['level1']][xx]
        })),
        rep(
          1:sum(sample_size_list[['level2']]),
          times = sample_size_list[['level1']]
        )
      )
    } else {
      id_vars <- data.frame(1:sample_size_list[['level1']])
    }
  }

  names(id_vars) <- id_names
  id_vars
}

compute_samplesize <- function(data, sim_args) {
  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']]]

  samp_size <- lapply(seq_along(id_vars), function(xx) {
    as.numeric(table(data[id_vars[xx]]))
  })

  if (length(id_vars) == 2) {
    level2_ss <- aggregate(
      as.formula(paste0(id_vars[1], "~", id_vars[2])),
      data = data,
      FUN = length_unique
    )[[2]]

    list(
      level1 = samp_size[[1]],
      level2 = level2_ss,
      level3 = sim_args[['sample_size']][['level3']],
      level3_total = samp_size[[2]]
    )
  } else {
    if (length(id_vars) == 1) {
      list(
        level1 = samp_size[[1]],
        level2 = sim_args[['sample_size']][['level2']]
      )
    } else {
      list(level1 = sim_args[['sample_size']])
    }
  }
}

length_unique <- function(x) length(unique(x))

is_factor_var <- function(sim_args) {
  unlist(lapply(seq_along(sim_args[['fixed']]), function(yy) {
    isTRUE(sim_args[['fixed']][[yy]][['var_type']] == 'factor')
  }))
}

factor_names <- function(sim_args, fixed_vars) {
  which_factor <- is_factor_var(sim_args)
  factor_names <- names(sim_args[['fixed']][which_factor])

  fixed_vars_continuous <- fixed_vars[
    !grepl(paste(factor_names, collapse = "|"), fixed_vars)
  ]
  fixed_vars_cat <- fixed_vars[grepl(
    paste0('^', paste(factor_names, collapse = "|"), "$"),
    fixed_vars
  )]

  if (any(grepl(":|^I", fixed_vars_cat))) {
    int_loc <- grep(":|^I", fixed_vars_cat)
    fixed_vars_cat <- fixed_vars_cat[-int_loc]
  }

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

  loc <- num_levels > 1
  fixed_levels_gt2 <- fixed_vars_cat[loc]
  num_levels_gt2 <- num_levels[loc]

  fixed_levels_gt2_names <- lapply(seq_along(fixed_levels_gt2), function(xx) {
    paste0(fixed_levels_gt2[xx], '_', 1:(num_levels_gt2[[xx]] - 1))
  })
  names(fixed_levels_gt2_names) <- fixed_vars_cat

  if (any(grepl(":|^I", fixed_vars))) {
    int_loc <- grep(":|^I", fixed_vars)
    var_loc <- lapply(seq_along(fixed_levels_gt2), function(xx) {
      grep(paste0('^', fixed_levels_gt2[xx], "$"), fixed_vars[-int_loc])
    })
  } else {
    var_loc <- lapply(seq_along(fixed_levels_gt2), function(xx) {
      grep(paste0('^', fixed_levels_gt2[xx], "$"), fixed_vars)
    })
  }

  updated_names <- lapply(seq_along(fixed_levels_gt2), function(ii) {
    lapply(seq_along(fixed_levels_gt2_names[[ii]]), function(xx) {
      gsub(
        fixed_levels_gt2[ii],
        fixed_levels_gt2_names[[ii]][xx],
        fixed_vars[var_loc[[ii]]]
      )
    })
  })

  reordered_names <- lapply(seq_along(fixed_levels_gt2), function(ii) {
    reorder_names(updated_names[[ii]])
  })

  if (any(grepl(":|^I", fixed_vars))) {
    fixed_vars_cat_rename <- fixed_vars[grepl(
      paste(factor_names, collapse = "|"),
      fixed_vars
    )]

    if (any(grepl("_post$", fixed_vars_cat_rename[int_loc]))) {
      fixed_vars_cat_rename <- fixed_vars_cat_rename[
        !grepl("_post$", fixed_vars_cat_rename)
      ]
    }
    int_loc <- grep(":|^I", fixed_vars_cat_rename)

    new_interaction_names <- interaction_names(
      fixed_vars_cat_rename,
      fixed_levels_gt2_names,
      sim_args
    )
    for (ii in seq_along(int_loc)) {
      # int_loc <- grep(":|^I", fixed_vars)
      int_loc_rename <- grep(
        paste0('^', fixed_vars_cat_rename[int_loc[ii]], '$', collapse = ""),
        fixed_vars
      )
      fixed_vars[int_loc_rename] <- new_interaction_names[ii]
    }
  }

  # gsub(fixed_vars_cat[1], reordered_names[[1]], fixed_vars)

  # fixed_vars <- c(fixed_vars_continuous, fixed_vars_cat_rename)

  for (ii in seq_along(var_loc)) {
    fixed_vars[[var_loc[[ii]]]] <- reordered_names[ii]
  }

  unlist(fixed_vars)
}

interaction_names <- function(fixed_vars, renamed_vars, sim_args) {
  int_loc <- grep(":|^I", fixed_vars)

  int_names <- lapply(int_loc, function(ii) {
    unlist(strsplit(fixed_vars[ii], split = ":"))
  })

  factor_vars <- unlist(lapply(seq_along(sim_args[['fixed']]), function(xx) {
    sim_args[['fixed']][[xx]][['var_type']] == 'factor'
  }))
  factor_names <- names(sim_args[['fixed']])[factor_vars]

  # int_vars_location_f <- lapply(seq_along(int_names), function(ii)
  #   grep(paste(int_names[[ii]][factor_vars], collapse = "|"), fixed_vars[-int_loc[ii]]))
  #
  renamed_int_vars <- lapply(seq_along(int_names), function(ii) {
    if (length(int_names[[ii]][int_names[[ii]] %in% factor_names]) == 1) {
      cont_vars <- int_names[[ii]][!int_names[[ii]] %in% factor_names]
      if (length(cont_vars) > 1) {
        cont_vars <- paste0(cont_vars, collapse = ":")
      }
      do.call(
        paste,
        c(
          expand.grid(
            cont_vars,
            renamed_vars[[int_names[[ii]][int_names[[ii]] %in% factor_names]]],
            stringsAsFactors = FALSE
          ),
          sep = ":"
        )
      )
    } else {
      if (all(int_names[[ii]] %in% factor_names)) {
        do.call(
          paste,
          c(
            expand.grid(
              renamed_vars[int_names[[ii]][int_names[[ii]] %in% factor_names]],
              stringsAsFactors = FALSE
            ),
            sep = ":"
          )
        )
      } else {
        cont_vars <- int_names[[ii]][!int_names[[ii]] %in% factor_names]
        factor_expand <- expand.grid(
          renamed_vars[int_names[[ii]][int_names[[ii]] %in% factor_names]],
          stringsAsFactors = FALSE
        )
        do.call(paste, c(cbind(cont_vars, factor_expand), sep = ":"))
      }
    }
  })

  renamed_int_vars
}

reorder_names <- function(names) {
  unlisted_names <- do.call('c', names)
  if (any(grepl(":|^I", unlisted_names))) {
    int_loc <- grep(":|^I", unlisted_names)
    c(unlisted_names[-int_loc], unlisted_names[int_loc])
  } else {
    unlisted_names
  }
}


poly_ns_names <- function(sim_args, fixed_vars) {
  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")

  ns_loc <- grep("^ns|^bs", fixed_vars)
  if (any(ns_loc)) {
    ns_new_names <- ns_df_names(fixed_vars[ns_loc])
    if (grepl("df", fixed_vars[ns_loc])) {
      ns_new_names <- lapply(seq_along(ns_loc), function(xx) {
        ns_df_names(fixed_vars[ns_loc[xx]])
      })
      for (i in seq_along(ns_new_names)) {
        fixed_vars[ns_loc[i]] <- ns_new_names[i]
      }
      fixed_vars <- unlist(fixed_vars)
    } else {
      # placeholder for bs names
    }
  }
  poly_loc <- grep("^poly", fixed_vars)
  if (any(poly_loc)) {
    poly_new_names <- lapply(seq_along(poly_loc), function(xx) {
      poly_names(fixed_vars[poly_loc[xx]])
    })
    for (i in seq_along(poly_new_names)) {
      fixed_vars[poly_loc[i]] <- poly_new_names[i]
    }
    fixed_vars <- unlist(fixed_vars)
  }

  fixed_vars
}

ns_df_names <- function(x) {
  name <- gsub("ns\\(|bs\\(|\\,.+\\)$", "", x)

  func_arg <- unlist(regmatches(x, regexec("df\\s+=\\s+[0-9]+", x)))
  num <- unlist(regmatches(func_arg, regexec("[0-9]+", func_arg)))

  paste(name, 1:as.numeric(num), sep = "_")
}

poly_names <- function(x) {
  degree_arg <- unlist(regmatches(x, regexec("degree\\s+=\\s+[0-9]+", x)))
  num <- unlist(regmatches(degree_arg, regexec("[0-9]+", degree_arg)))
  name <- gsub("poly\\(|\\,.+\\)", "", x)

  paste(name, 1:as.numeric(num), sep = "_")
}

basis_term_names <- function(x) {
  if (grepl("^poly\\(", x)) {
    return(poly_names(x))
  }
  if (grepl("^ns\\(|^bs\\(", x)) {
    return(ns_df_names(x))
  }

  x
}

uses_spline_basis <- function(formula) {
  any(grepl("^ns\\(|^bs\\(", attr(terms(formula), "term.labels")))
}

attach_spline_env <- function(formula) {
  if (is.list(formula)) {
    return(lapply(formula, attach_spline_env))
  }

  if (!inherits(formula, "formula") || !uses_spline_basis(formula)) {
    return(formula)
  }

  formula_env <- new.env(parent = environment(formula))
  formula_env$ns <- splines::ns
  formula_env$bs <- splines::bs
  environment(formula) <- formula_env

  formula
}

regression_weights <- function(reg_weights, x_names, fixed_formula) {
  if (!is.list(reg_weights) || is.null(names(reg_weights))) {
    return(reg_weights)
  }

  if (any(names(reg_weights) == "")) {
    return(reg_weights)
  }

  term_labels <- attr(terms(fixed_formula), "term.labels")
  weights <- stats::setNames(rep(NA_real_, length(x_names)), x_names)

  for (term in names(reg_weights)) {
    term_weights <- unlist(reg_weights[[term]], use.names = FALSE)
    if (term %in% c("(Intercept)", "Intercept", "X.Intercept.")) {
      term_cols <- "X.Intercept."
    } else if (term %in% x_names) {
      term_cols <- term
    } else if (term %in% term_labels) {
      term_cols <- basis_term_names(term)
    } else {
      term_cols <- paste(term, seq_along(term_weights), sep = "_")
    }

    if (length(term_cols) != length(term_weights)) {
      stop(
        "Regression weights for `",
        term,
        "` must have one value per generated model column.",
        call. = FALSE
      )
    }

    missing_cols <- term_cols[term_cols %ni% x_names]
    if (length(missing_cols) > 0) {
      stop(
        "Regression weights for `",
        term,
        "` refer to generated columns that were not found: ",
        paste(missing_cols, collapse = ", "),
        call. = FALSE
      )
    }

    weights[term_cols] <- term_weights
  }

  if (any(is.na(weights))) {
    stop(
      "`reg_weights` must specify weights for all fixed-effect columns: ",
      paste(names(weights)[is.na(weights)], collapse = ", "),
      call. = FALSE
    )
  }

  unname(weights)
}

regression_weight_terms <- function(reg_weights, fixed_formula) {
  if (!is.list(reg_weights) || is.null(names(reg_weights))) {
    return(reg_weights)
  }

  if (any(names(reg_weights) == "")) {
    return(reg_weights)
  }

  term_labels <- attr(terms(fixed_formula), "term.labels")
  weights <- numeric()

  for (term in names(reg_weights)) {
    term_weights <- unlist(reg_weights[[term]], use.names = FALSE)
    if (term %in% c("(Intercept)", "Intercept", "X.Intercept.")) {
      term_names <- "(Intercept)"
    } else if (term %in% term_labels && grepl("^ns\\(|^bs\\(", term)) {
      term_names <- paste0(term, seq_along(term_weights))
    } else if (term %in% term_labels && grepl("^poly\\(", term)) {
      term_names <- paste0(term, seq_along(term_weights))
    } else if (length(term_weights) == 1) {
      term_names <- term
    } else {
      term_names <- paste(term, seq_along(term_weights), sep = "_")
    }

    if (length(term_names) != length(term_weights)) {
      stop(
        "Regression weights for `",
        term,
        "` must have one value per fitted model coefficient.",
        call. = FALSE
      )
    }

    weights <- c(weights, stats::setNames(term_weights, term_names))
  }

  weights
}


unique_columns <- function(x, y) {
  sapply(names(y), function(nm) {
    idx <- grep(paste0("^", nm, "$"), names(x))
    if (length(idx) == 0L) {
      # Column only in y: keep it
      TRUE
    } else {
      any(y[[nm]] != x[[idx[1]]])
    }
  })
}


dataframe2matrix <- function(data, corr_variable, var_names) {
  if (is.null(data)) {
    NULL
  } else {
    n <- (sqrt(1 + 8 * nrow(data)) + 1) / 2

    corr_mat <- matrix(NA, nrow = n, ncol = n)
    lower <- lower.tri(corr_mat, diag = FALSE)
    upper <- upper.tri(corr_mat, diag = FALSE)

    corr_mat[lower] <- data[[corr_variable]]
    corr_mat[upper] <- data[[corr_variable]]
    diag(corr_mat) <- 1

    mat_names <- unique(unlist(data[var_names]))

    colnames(corr_mat) <- mat_names
    rownames(corr_mat) <- mat_names

    corr_mat
  }
}

density_quantile <- function(data, quantiles) {
  data_density <- data.frame(density(data, na.rm = TRUE, n = 10000)[c(
    'x',
    'y'
  )])

  matches <- lapply(quantiles, closest_match, data = data_density$x)

  data_density[unlist(matches), ]
}

closest_match <- function(data, value) {
  which.min(abs(data - value))
}

list_select <- function(list, names, exclude = TRUE, simplify = FALSE) {
  if (exclude) {
    index <- which(!(names(list) %in% names))
  } else {
    index <- which(names(list) %in% names)
  }
  if (simplify) {
    list[[index]]
  } else {
    list[index]
  }
}

# Horrible hack to keep CRAN happy and suppress NOTES about
# parts of the code that use non-standard evaluation.
# See:
# http://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
# https://github.com/smbache/magrittr/issues/29
utils::globalVariables(c(
  'test_stat',
  'reject',
  'estimate',
  'term',
  'std.error',
  '.',
  'sd_estimate',
  'avg_se',
  't1e',
  'power_args',
  'adjusted_teststat',
  't1e_args',
  'param_estimate_sd',
  'avg_standard_error',
  'statistic',
  'n',
  'miss_prob',
  'miss_prop',
  'test_statistic',
  'group',
  'ID',
  'p.value'
))