R/sim_data.R

In dalejbarr/funfact: Functions for factorial experiments

#' Generate population for data simulation
#'
#' @param design_args list specifying experimental design (see
#'   \code{\link{stim_lists}}); for best results, n_item should be
#'   2*minimum number of items needed for the design
#' @param n_subj number of subjects (for defining random effects
#'   structure; for best results, should be 2*number of stimulus lists
#' @param fixed_ranges list of 2-element vectors (min, max) defining
#'   ranges of fixed effect parameters
#' @param var_range 2-element vector defining range (min, max) of
#'   random effect variances
#' @param err_range 2-element vector defining range of error variance
#' @return list with parameters for data generation
#' @seealso \code{\link{sim_norm}}
#' @examples
#'
#' dargs <- list(ivs = c(A = 2, B = 2), n_item = 8)
#' gen_pop(dargs, 8)
#' @importFrom clusterGeneration genPositiveDefMat
#' @export
gen_pop <- function(design_args,
                    n_subj = NULL,
                    fixed_ranges = NULL,                    
                    var_range = c(0, 3),
                    err_range = c(0, 6)) {
    tdat <- trial_lists(design_args, n_subj)
    if ("n_rep" %in% colnames(tdat)) {
        design_args[["ivs"]] <- c(as.list(design_args[["ivs"]]), list(n_rep = unique(tdat[["n_rep"]])))
    } else {}
    forms <- design_formula(design_args, n_subj, lme4_format = FALSE)
    tnames <- lapply(forms, function(x) term_names(design_args, x))
    err_var <- runif(1, err_range[1], err_range[2])
    ## generate random effects
    rfx <- lapply(tnames[-1], function(x) {
        mx <- clusterGeneration::genPositiveDefMat(length(x),
                                                   covMethod = "onion",
                                                   rangeVar = var_range)$Sigma
        dimnames(mx) <- list(x, x)
        return(mx)
    })
    ## generate fixed effects
    if (!is.null(fixed_ranges)) {
        if (length(fixed_ranges) != length(tnames[["fixed"]])) {
            stop("fixed_ranges must have ", length(tnames[["fixed"]]),
                 " elements, corresponding to variables ",
                 paste(tnames[["fixed"]], collapse = ", "))
        } else {}
    } else {
        fixed_ranges <- lapply(seq_along(tnames[["fixed"]]), function(x) {c(0, 3)})
        names(fixed_ranges) <- tnames[["fixed"]]        
    }
    if (!is.null(names(fixed_ranges)) &&
            !identical(names(fixed_ranges), tnames[["fixed"]])) {
        warning("names of 'fixed_ranges' elements do not match variable names: ",
                paste(tnames[["fixed"]], collapse = ", "))
    } else {}
    return(list(fixed = sapply(fixed_ranges, function(x) runif(1, x[1], x[2])),
                subj_rfx = rfx[["subj_id"]],
                item_rfx = rfx[["item_id"]],
                err_var = err_var))
}

#' Sample data from population with normal error variance
#'
#' @param design_args experiment design (see
#'   \code{\link{stim_lists}}); must contain sub-element \code{n_item}
#'   defining the number of items.
#' @param n_subj number of subjects
#' @param params list with parameters defining population
#'   (\code{fixed}, \code{subj_rfx}, \code{item_rfx}), and
#'   \code{err_var}; normally generated by a call to
#'   \code{\link{gen_pop}}.
#' @param contr_type name of function defining IV contrast type; see
#'   \code{?contrasts}
#' @param verbose whether to return explicit information about
#'   individual random effects and residual error
#' 
#' @return A data frame with simulated data
#' @seealso \code{\link{gen_pop}}
#' @examples
#' design_args <- list(ivs = c(A = 2, B = 3), n_item = 18)
#' pop_params <- gen_pop(design_args, 12)
#' 
#' dat <- sim_norm(design_args, 12, pop_params)
#' @importFrom MASS mvrnorm
#' @export
sim_norm <- function(design_args,
                     n_subj,
                     params,
                     contr_type = "contr.dev",
                     verbose = FALSE) {
  required_elements <- c("fixed", "subj_rfx", "item_rfx", "err_var")
  missing_elements <- setdiff(required_elements, names(params))
  if (length(missing_elements) > 0) {
    stop("mcr.data was missing element(s): ",
         paste(missing_elements, collapse = ", "))
  } else {}
  if (is.null(design_args[["n_item"]])) {
    stop("'n_item' not specified in 'design_args'")
  } else {}
  rfx <- mapply(function(x, n) {
    if (!is.null(x)) {
      MASS::mvrnorm(n, mu = rep(0, ncol(x)), x)
    } else {
      n
    }
  }, params[c("subj_rfx", "item_rfx")],
  c(n_subj, design_args[["n_item"]]), SIMPLIFY = FALSE)
  dat <- compose_data(design_args,
                      fixed = params[["fixed"]],
                      subj_rmx = rfx[["subj_rfx"]],
                      item_rmx = rfx[["item_rfx"]],
                      contr_type = contr_type, verbose = verbose)
  dat[["err"]] <- rnorm(nrow(dat), sd = sqrt(params[["err_var"]]))
  dat[["Y"]] <- dat[["Y"]] + dat[["err"]]
  return(dat)
}