vignettes/workflow.R

In AFialkowski/SimMultiCorrData: Simulation of Correlated Data with Multiple Variable Types

## ---- echo = FALSE-------------------------------------------------------
#knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
knitr::opts_chunk$set(fig.width = 6, fig.height = 4.5) 

## ---- warning = FALSE, message = FALSE-----------------------------------
library("SimMultiCorrData")
library("printr")

# Turn off scientific notation
options(scipen = 999)

# Set seed and sample size
seed <- 11
n <- 10000

# Continuous Distributions
Dist <- c("Gaussian", "Chisq", "Beta")

# Calculate standardized cumulants
# Those for the normal distribution are rounded to ensure the correct values 
# are obtained.
M1 <- round(calc_theory(Dist = "Gaussian", params = c(0, 1)), 8)
M2 <- calc_theory(Dist = "Chisq", params = 4)
M3 <- calc_theory(Dist = "Beta", params = c(4, 2))
M <- cbind(M1, M2, M3)

# Binary and Ordinal Distributions
marginal <- list(c(0.3, 0.75), c(0.2, 0.5, 0.9))
support <- list() # default support will be generated inside simulation

# Poisson Distributions
lam <- c(1, 5, 10)

# Negative Binomial Distributions
size <- c(3, 6)
prob <- c(0.2, 0.8)

ncat <- length(marginal)
ncont <- ncol(M)
npois <- length(lam)
nnb <- length(size)

# Create correlation matrix from a uniform distribution (0.2, 0.7)
set.seed(seed)
Rey <- diag(1, nrow = (ncat + ncont + npois + nnb))
for (i in 1:nrow(Rey)) {
  for (j in 1:ncol(Rey)) {
    if (i > j) Rey[i, j] <- runif(1, 0.2, 0.7)
    Rey[j, i] <- Rey[i, j]
  }
}

# Check to see if Rey is positive-definite
min(eigen(Rey, symmetric = TRUE)$values) < 0

## ---- warning = FALSE----------------------------------------------------
Lower <- list()

# list of standardized kurtosis values to add in case only invalid power 
#     method pdfs are produced
Skurt <- list(seq(0.5, 2, 0.5), seq(0.02, 0.05, 0.01), seq(0.02, 0.05, 0.01))

start.time <- Sys.time()
for (i in 1:ncol(M)) {
  Lower[[i]] <- calc_lower_skurt(method = "Polynomial", skews = M[3, i], 
                                 fifths = M[5, i], sixths = M[6, i], 
                                 Skurt = Skurt[[i]], seed = 104)
}
stop.time <- Sys.time()
Time <- round(difftime(stop.time, start.time, units = "min"), 3)
cat("Total computation time:", Time, "minutes \n")

# Note the message given for Distribution 1 (Normal).

## ------------------------------------------------------------------------
as.matrix(Lower[[1]]$Min[1, c("skew", "fifth", "sixth", "valid.pdf", 
                              "skurtosis")], 
          nrow = 1, ncol = 5, byrow = TRUE) 

## ------------------------------------------------------------------------
as.matrix(Lower[[2]]$Min[1, c("skew", "fifth", "sixth", "valid.pdf", 
                              "skurtosis")], 
          nrow = 1, ncol = 5, byrow = TRUE) 
Lower[[2]]$SkurtCorr1

## ------------------------------------------------------------------------
as.matrix(Lower[[3]]$Min[1, c("skew", "fifth", "sixth", "valid.pdf", 
                              "skurtosis")], 
          nrow = 1, ncol = 5, byrow = TRUE) 
Lower[[3]]$SkurtCorr1

## ---- warning = FALSE----------------------------------------------------
# Make sure Rey is within upper and lower correlation limits
valid <- valid_corr(k_cat = ncat, k_cont = ncont, k_pois = npois,
                    k_nb = nnb, method = "Polynomial", means =  M[1, ],
                    vars =  (M[2, ])^2, skews = M[3, ], skurts = M[4, ],
                    fifths = M[5, ], sixths = M[6, ], marginal = marginal, 
                    lam = lam, size = size, prob = prob, rho = Rey, 
                    seed = seed)

## ---- warning = FALSE, message = FALSE-----------------------------------
A <- rcorrvar(n = 10000, k_cont = ncont, k_cat = ncat, k_pois = npois,
              k_nb = nnb, method = "Polynomial", means =  M[1, ], 
              vars =  (M[2, ])^2, skews = M[3, ], skurts = M[4, ], 
              fifths = M[5, ], sixths = M[6, ], marginal = marginal,
              lam = lam, size = size, prob = prob, rho = Rey, seed = seed)

## ------------------------------------------------------------------------
Acorr_error = round(A$correlations - Rey, 6)
summary(as.numeric(Acorr_error))

## ---- warning = FALSE, message = FALSE-----------------------------------
B <- rcorrvar(n = 10000, k_cont = ncont, k_cat = ncat, k_pois = npois,
              k_nb = nnb, method = "Polynomial", means =  M[1, ], 
              vars =  (M[2, ])^2, skews = M[3, ], skurts = M[4, ], 
              fifths = M[5, ], sixths = M[6, ], marginal = marginal,
              lam = lam, size = size, prob = prob, rho = Rey, seed = seed, 
              errorloop = TRUE)

## ------------------------------------------------------------------------
Bcorr_error = round(B$correlations - Rey, 6)
summary(as.numeric(Bcorr_error))

## ------------------------------------------------------------------------
knitr::kable(B$summary_ordinal[[1]], caption = "Variable 1")
knitr::kable(B$summary_ordinal[[2]], caption = "Variable 2")

## ------------------------------------------------------------------------
as.matrix(B$summary_Poisson[, c(1, 3:6, 8:9)], nrow = 3, ncol = 7, 
          byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(B$summary_Neg_Bin[, c(1, 3:7, 9:10)], nrow = 2, ncol = 8, 
          byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(round(B$constants, 6), nrow = 3, ncol = 6, byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(round(B$summary_targetcont, 5), nrow = 3, ncol = 7, byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(round(B$summary_continuous[, c("Distribution", "mean", "sd", 
                                         "skew", "skurtosis", "fifth", 
                                         "sixth")], 5), nrow = 3, ncol = 7, 
          byrow = TRUE)

## ------------------------------------------------------------------------
B$valid.pdf

## ---- warning = FALSE, message = FALSE-----------------------------------
as.matrix(t(round(stats_pdf(c = B$constants[1, ], method = "Polynomial", 
                            alpha = 0.025), 4)))

## ---- warning = FALSE, message = FALSE-----------------------------------
as.matrix(t(round(stats_pdf(c = B$constants[2, ], method = "Polynomial", 
                            alpha = 0.025), 4)))

## ---- warning = FALSE, message = FALSE-----------------------------------
as.matrix(t(round(stats_pdf(c = B$constants[3, ], method = "Polynomial", 
                            alpha = 0.025), 4)))

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_cdf(B$continuous_variables[, 2], calc_cprob = TRUE, delta = 10)

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_pdf_theory(B$continuous_variables[, 2], Dist = "Chisq", params = 4)

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_cdf(B$ordinal_variables[, 2])

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_theory(B$Poisson_variables[, 2], cont_var = FALSE, Dist = "Poisson", 
                params = 5)

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_pdf_theory(B$Poisson_variables[, 2], cont_var = FALSE, 
                    Dist = "Poisson", params = 5)

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_theory(B$Neg_Bin_variables[, 1], cont_var = FALSE, 
                Dist = "Negative_Binomial", params = c(3, 0.2))

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_pdf_theory(B$Neg_Bin_variables[, 1], cont_var = FALSE, 
                Dist = "Negative_Binomial", params = c(3, 0.2))

## ---- warning = FALSE----------------------------------------------------
pois_eps <- rep(0.0001, npois)
nb_eps <- rep(0.0001, nnb)

# Make sure Rey is within upper and lower correlation limits
valid2 <- valid_corr2(k_cat = ncat, k_cont = ncont, k_pois = npois,
                      k_nb = nnb, method = "Polynomial", means =  M[1, ],
                      vars =  (M[2, ])^2, skews = M[3, ], skurts = M[4, ],
                      fifths = M[5, ], sixths = M[6, ], marginal = marginal, 
                      lam = lam, pois_eps = pois_eps, size = size, 
                      prob = prob, nb_eps = nb_eps, rho = Rey, seed = seed)

## ---- warning = FALSE, message = FALSE-----------------------------------
C <- rcorrvar2(n = 10000, k_cont = ncont, k_cat = ncat, k_pois = npois,
               k_nb = nnb, method = "Polynomial", means =  M[1, ], 
               vars =  (M[2, ])^2, skews = M[3, ], skurts = M[4, ], 
               fifths = M[5, ], sixths = M[6, ], marginal = marginal, 
               lam = lam, pois_eps = pois_eps, size = size, prob = prob, 
               nb_eps = nb_eps, rho = Rey, seed = seed)

## ------------------------------------------------------------------------
Ccorr_error = round(C$correlations - Rey, 6)
summary(as.numeric(Ccorr_error))

## ---- warning = FALSE, message = FALSE-----------------------------------
D <- rcorrvar2(n = 10000, k_cont = ncont, k_cat = ncat, k_pois = npois,
               k_nb = nnb, method = "Polynomial", means =  M[1, ], 
               vars =  (M[2, ])^2, skews = M[3, ], skurts = M[4, ], 
               fifths = M[5, ], sixths = M[6, ], marginal = marginal, 
               lam = lam, pois_eps = pois_eps, size = size, prob = prob, 
               nb_eps = nb_eps, rho = Rey, seed = seed, errorloop = TRUE)

## ------------------------------------------------------------------------
Dcorr_error = round(D$correlations - Rey, 6)
summary(as.numeric(Dcorr_error))

## ------------------------------------------------------------------------
knitr::kable(D$summary_ordinal[[1]], caption = "Variable 1")
knitr::kable(D$summary_ordinal[[2]], caption = "Variable 2")

## ------------------------------------------------------------------------
as.matrix(D$summary_Poisson[, c(1, 3:6, 8:9)], nrow = 3, ncol = 7, 
          byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(D$summary_Neg_Bin[, c(1, 3:7, 9:10)], nrow = 2, ncol = 8, 
          byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(round(D$summary_targetcont, 5), nrow = 3, ncol = 7, byrow = TRUE)

## ------------------------------------------------------------------------
as.matrix(round(D$summary_continuous[, c("Distribution", "mean", "sd", 
                                         "skew", "skurtosis", "fifth", 
                                         "sixth")], 5), nrow = 3, ncol = 7, 
          byrow = TRUE)

## ------------------------------------------------------------------------
D$valid.pdf

## ---- warning = FALSE, message = FALSE-----------------------------------
as.matrix(t(round(stats_pdf(c = D$constants[1, ], method = "Polynomial", 
                            alpha = 0.025), 4)))

## ---- warning = FALSE, message = FALSE-----------------------------------
as.matrix(t(round(stats_pdf(c = B$constants[2, ], method = "Polynomial", 
                            alpha = 0.025), 4)))

## ---- warning = FALSE, message = FALSE-----------------------------------
as.matrix(t(round(stats_pdf(c = B$constants[3, ], method = "Polynomial", 
                            alpha = 0.025), 4)))

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_theory(D$Poisson_variables[, 2], cont_var = FALSE, Dist = "Poisson", 
                params = 5)

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_pdf_theory(D$Poisson_variables[, 2], cont_var = FALSE, 
                    Dist = "Poisson", params = 5)

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_theory(D$Neg_Bin_variables[, 1], cont_var = FALSE, 
                Dist = "Negative_Binomial", params = c(3, 0.2))

## ---- warning = FALSE, message = FALSE-----------------------------------
plot_sim_pdf_theory(D$Neg_Bin_variables[, 1], cont_var = FALSE, 
                    Dist = "Negative_Binomial", params = c(3, 0.2))