R/plotextra.R

In simulatr/simrel: Simulation of Multivariate Linear Model Data

#' @title Tidy Functions to make plotting easy
#' @name tidy_beta
#' @param obj A Simrel Object
#' @return A tibble with three columns: Predictor, Response and BetaCoef
#' @importFrom purrr pluck modify_at modify
#' @importFrom tibble as_tibble
#' @importFrom tidyr gather_
#' @examples
#' sobj <- multisimrel()
#' beta_df <- tidy_beta(sobj)
#' beta_df
#' @rdname tidy_beta
#' @export
tidy_beta <- function(obj) {
  beta <- obj %>% pluck("beta")
  dimnames(beta) <- list(
    c(1:nrow(beta)), 
    c(paste0("_", 1:ncol(beta)))
  )
  beta_df <- beta %>% 
    as_tibble(rownames = "Predictor")
  beta_df %>% 
    gather_('Response', 'BetaCoef', names(beta_df)[-1]) %>% 
    modify_at("Response", ~gsub("_", "", ..1)) %>% 
    modify_at(1:2, as.integer)
}

#' @title Extract various sigma matrices
#' @name cov_mat
#' @param obj A simrel object
#' @param which A character string to specify which covariance matrix to extract, possible values are "xy", "zy" and "zw"
#' @param use_population A boolean whether to use compute population values or to estimate from sample
#' @return A matrix of covariances with column equals to the number of response and row equals to the number of predictors
#' @importFrom purrr pluck
#' @examples
#' set.seed(1983)
#' sobj <- multisimrel()
#' cov_mat(sobj, which = "xy", use_population = TRUE)
#' cov_mat(sobj, which = "xy", use_population = FALSE)
#' @rdname cov_mat
#' @export
cov_mat <- function(obj, which = c("xy", "zy", "zw"), use_population = TRUE) {
  switch (which,
          xy = cov_xy(obj, use_population),
          zy = cov_zy(obj, use_population),
          zw = cov_zw(obj)
  )
}

#' @name tidy_sigma
#' @title Tidy covariance matrix
#' @param covs A sigma matrix obtained from cov_mat function
#' @return A tibble with three columns: Predictor, Response and Covariance
#' @importFrom purrr pluck modify_at
#' @importFrom tibble as_tibble
#' @importFrom tidyr gather_
#' @examples
#' sobj <- multisimrel()
#  sobj %>%
#     cov_mat("zy") %>%
#     tidy_sigma()
#' @rdname tidy_sigma
#' @export
tidy_sigma <- function(covs) {
  n_pred <- nrow(covs)
  n_resp <- ncol(covs)
  dimnames(covs) <- list(1:n_pred, paste0("_", 1:n_resp))
  covs_df <- covs %>% 
    as_tibble(rownames = "Predictor")
  covs_df %>% 
    gather_('Response', 'Covariance', names(covs_df)[-1]) %>% 
    modify_at("Response", ~gsub("_", "", ..1)) %>% 
    modify_at(1:2, as.integer)
}

#' @title Absolute value of sigma scaled by the overall maximum absolute value
#' @name abs_sigma
#' @param sigma_df A tidy covariance data frame generated by tidy_sigma function
#' @return Another data.frame (tibble) of same dimension with absolute covarinace scaled by overall maximum absolute values
#' @importFrom purrr modify_at
#' @examples
#' sobj <- multisimrel()
#' sobj %>% 
#'     cov_mat("zy") %>% 
#'     tidy_sigma() %>% 
#'     abs_sigma()
#' @rdname tidy_beta
#' @export
abs_sigma <- function(sigma_df) {
  sigma_df %>% 
    modify_at("Covariance", ~abs(..1)/max(abs(..1)))
}

#' @name tidy_lambda
#' @title Extract Eigenvalues of predictors
#' @param obj A simrel Object
#' @param use_population A boolean to specify where to use population value or calculate from sample
#' @return A dataframe of eigenvalues for each predictors
#' @importFrom purrr pluck modify_at
#' @importFrom tibble tibble
#' @examples
#' sobj <- multisimrel()
#' sobj %>% 
#'     tidy_lambda()
#' @rdname tidy_lambda
#' @export
#' 
tidy_lambda <- function(obj, use_population = TRUE) {
  if (use_population) {
    lmd <- obj %>% pluck("lambda")
  } else {
    svdres <- svd(obj %>% pluck("X"))
    lmd <- (svdres$d ^ 2)/(obj$n - 1)
    lmd <- lmd/lmd[1]
  }
  tibble(
    Predictor = seq_along(lmd),
    lambda = lmd
  )
}

#' Helper Functions
#' @keywords intgernal
#' @title Covariance between Z and W
#' @name cov_zw
#' @param obj A simrel object
#' @return A covariance matrix of Z and W
cov_zw <- function(obj) {
  type <- obj %>% pluck("type")
  if (type != "multivariate") {
    stop("\n\nUnivariate  and bivariate simulation does not have response component.",
         "\nUse cov_zy function instead.")
  } else {
    idx <- obj$m
    covs <-obj$SigmaWZ[-c(1:idx), 1:idx, drop = FALSE] 
    return(covs)
  }
}

#' @keywords intgernal
#' @title Covariance between Z and Y
#' @name cov_zy
#' @param obj A simrel object
#' @param use_population A boolean to specify wheather to use population or sample
#' @return A covariance matrix of Z and Y
cov_zy <- function(obj, use_population = TRUE) {
  if (use_population) {
    type <- obj %>% pluck("type")
    if (type == "multivariate") {
      covs <- t(obj$SigmaYZ)
    } else {
      m <- ifelse(type == "univariate", 1, 2)
      covs <- obj$Sigma[-c(1:m), 1:m, drop = FALSE]
    }
  } else {
    X <- scale(obj$X, center = TRUE, scale = FALSE)
    Y <- scale(obj$Y, center = TRUE, scale = FALSE)
    svdres <- svd(X)
    Z <- X %*% svdres$v
    covs <- t(cov(Y, Z))
  }
  return(unname(covs))
}

#' @keywords intgernal
#' @title Covariance between X and Y
#' @name cov_xy
#' @param obj A simrel object
#' @param use_population A boolean to specify wheather to use population or sample
#' @return A covariance matrix of X and Y
cov_xy <- function(obj, use_population = TRUE) {
  covs <- if (use_population) {
    covs <- if(obj$type == "multivariate") {
      t(obj$SigmaYX)
    } else {
      m <- ifelse(obj$type == 'univariate', 1, 2)
      p <- obj$p
      rotation <- obj$Rotation
      covs <- obj$Sigma[-c(1:m), c(1:m), drop = FALSE]
      covs <- t(rotation) %*% covs
    }
  } else {
    X <- scale(obj$X, center = TRUE, scale = FALSE)
    Y <- scale(obj$Y, center = TRUE, scale = FALSE)
    covs <- t(cov(Y, X))
  }
  return(unname(covs))
}