# Determine Correlation Bounds for Ordinal, Continuous, Poisson, and/or Negative Binomial Variables: Correlation Method 2 This function calculates the lower and upper correlation bounds for the given distributions and checks if a given target correlation matrix `rho` is within the bounds. It should be used before simulation with `corrvar2`. However, even if all pairwise correlations fall within the bounds, it is still possible that the desired correlation matrix is not feasible. This is particularly true when ordinal variables ($r \ge 2$ categories) are generated or negative correlations are desired. Therefore, this function should be used as a general check to eliminate pairwise correlations that are obviously not reproducible. It will help prevent errors when executing the simulation. The **ordering** of the variables in `rho` must be 1st ordinal, 2nd continuous non-mixture, 3rd components of continuous mixture, 4th regular Poisson, 5th zero-inflated Poisson, 6th regular NB, and 7th zero-inflated NB. Note that it is possible for `k_cat`, `k_cont`, `k_mix`, `k_pois`, and/or `k_nb` to be 0. The target correlations are specified with respect to the components of the continuous mixture variables. There are no parameter input checks in order to decrease simulation time. All inputs should be checked prior to simulation with `validpar`. Please see the Comparison of Correlation Methods 1 and 2 vignette for the differences between the two correlation methods, and the Variable Types vignette for a detailed explanation of how the correlation boundaries are calculated. ```r validcorr2(n = 10000, k_cat = 0, k_cont = 0, k_mix = 0, k_pois = 0, k_nb = 0, method = c("Fleishman", "Polynomial"), means = NULL, vars = NULL, skews = NULL, skurts = NULL, fifths = NULL, sixths = NULL, Six = list(), mix_pis = list(), mix_mus = list(), mix_sigmas = list(), mix_skews = list(), mix_skurts = list(), mix_fifths = list(), mix_sixths = list(), mix_Six = list(), marginal = list(), lam = NULL, p_zip = 0, size = NULL, prob = NULL, mu = NULL, p_zinb = 0, pois_eps = 0.0001, nb_eps = 0.0001, rho = NULL, seed = 1234, use.nearPD = TRUE, quiet = FALSE) ``` ## Arguments - `n`: the sample size (i.e. the length of each simulated variable; default = 10000) - `k_cat`: the number of ordinal (r >= 2 categories) variables (default = 0) - `k_cont`: the number of continuous non-mixture variables (default = 0) - `k_mix`: the number of continuous mixture variables (default = 0) - `k_pois`: the number of regular Poisson and zero-inflated Poisson variables (default = 0) - `k_nb`: the number of regular Negative Binomial and zero-inflated Negative Binomial variables (default = 0) - `method`: the method used to generate the k_cont non-mixture and k_mix mixture continuous variables. "Fleishman" uses Fleishman's third-order polynomial transformation and "Polynomial" uses Headrick's fifth-order transformation. - `means`: a vector of means for the k_cont non-mixture and k_mix mixture continuous variables (i.e. `rep(0, (k_cont + k_mix))`) - `vars`: a vector of variances for the k_cont non-mixture and k_mix mixture continuous variables (i.e. `rep(1, (k_cont + k_mix))`) - `skews`: a vector of skewness values for the `k_cont` non-mixture continuous variables - `skurts`: a vector of standardized kurtoses (kurtosis - 3, so that normal variables have a value of 0) for the `k_cont` non-mixture continuous variables - `fifths`: a vector of standardized fifth cumulants for the `k_cont` non-mixture continuous variables (not necessary for `method` = "Fleishman") - `sixths`: a vector of standardized sixth cumulants for the `k_cont` non-mixture continuous variables (not necessary for `method` = "Fleishman") - `Six`: a list of vectors of sixth cumulant correction values for the `k_cont` non-mixture continuous variables if no valid PDF constants are found, ex: `Six = list(seq(0.01, 2, 0.01), seq(1, 10, 0.5))`; if no correction is desired for variable $Y_{cont_i}$, set set the i-th list component equal to `NULL`; if no correction is desired for any of the $Y_{cont}$ keep as `Six = list()` (not necessary for `method` = "Fleishman") - `mix_pis`: a list of length `k_mix` with i-th component a vector of mixing probabilities that sum to 1 for component distributions of $Y_{mix_i}$ - `mix_mus`: a list of length `k_mix` with i-th component a vector of means for component distributions of $Y_{mix_i}$ - `mix_sigmas`: a list of length `k_mix` with i-th component a vector of standard deviations for component distributions of $Y_{mix_i}$ - `mix_skews`: a list of length `k_mix` with i-th component a vector of skew values for component distributions of $Y_{mix_i}$ - `mix_skurts`: a list of length `k_mix` with i-th component a vector of standardized kurtoses for component distributions of $Y_{mix_i}$ - `mix_fifths`: a list of length `k_mix` with i-th component a vector of standardized fifth cumulants for component distributions of $Y_{mix_i}$ (not necessary for `method` = "Fleishman") - `mix_sixths`: a list of length `k_mix` with i-th component a vector of standardized sixth cumulants for component distributions of $Y_{mix_i}$ (not necessary for `method` = "Fleishman") - `mix_Six`: a list of length `k_mix` with i-th component a list of vectors of sixth cumulant correction values for component distributions of $Y_{mix_i}$; use `NULL` if no correction is desired for a given component or mixture variable; if no correction is desired for any of the $Y_{mix}$ keep as `mix_Six = list()` (not necessary for `method` = "Fleishman") - `marginal`: a list of length equal to `k_cat`; the i-th element is a vector of the cumulative probabilities defining the marginal distribution of the i-th variable; if the variable can take r values, the vector will contain r - 1 probabilities (the r-th is assumed to be 1); for binary variables, these should be input the same as for ordinal variables with more than 2 categories (i.e. the user-specified probability is the probability of the 1st category, which has the smaller support value) - `lam`: a vector of lambda (> 0) constants for the Poisson variables (see `stats::dpois`); the order should be 1st regular Poisson variables, 2nd zero-inflated Poisson variables - `p_zip`: a vector of probabilities of structural zeros (not including zeros from the Poisson distribution) for the zero-inflated Poisson variables (see `VGAM::dzipois`); if `p_zip` = 0, $Y_{pois}$ has a regular Poisson distribution; if `p_zip` is in (0, 1), $Y_{pois}$ has a zero-inflated Poisson distribution; if `p_zip` is in `(-(exp(lam) - 1)^(-1), 0)`, $Y_{pois}$ has a zero-deflated Poisson distribution and `p_zip` is not a probability; if `p_zip = -(exp(lam) - 1)^(-1)`, $Y_{pois}$ has a positive-Poisson distribution (see `VGAM::dpospois`); if `length(p_zip) < length(lam)`, the missing values are set to 0 (and ordered 1st) - `size`: a vector of size parameters for the Negative Binomial variables (see `stats::dnbinom`); the order should be 1st regular NB variables, 2nd zero-inflated NB variables - `prob`: a vector of success probability parameters for the NB variables; order the same as in `size` - `mu`: a vector of mean parameters for the NB variables (*Note: either `prob` or `mu` should be supplied for all Negative Binomial variables, not a mixture; default = NULL); order the same as in `size`; for zero-inflated NB this refers to the mean of the NB distribution (see `VGAM::dzinegbin`) - `p_zinb`: a vector of probabilities of structural zeros (not including zeros from the NB distribution) for the zero-inflated NB variables (see `VGAM::dzinegbin`); if `p_zinb` = 0, $Y_{nb}$ has a regular NB distribution; if `p_zinb` is in `(-prob^size/(1 - prob^size),` `0)`, $Y_{nb}$ has a zero-deflated NB distribution and `p_zinb` is not a probability; if `p_zinb = -prob^size/(1 - prob^size)`, $Y_{nb}$ has a positive-NB distribution (see `VGAM::dposnegbin`); if `length(p_zinb) < length(size)`, the missing values are set to 0 (and ordered 1st) - `pois_eps`: a vector of length `k_pois` containing total cumulative probability truncation values; if none are provided, the default is 0.0001 for each variable - `nb_eps`: a vector of length `k_nb` containing total cumulative probability truncation values; if none are provided, the default is 0.0001 for each variable - `rho`: the target correlation matrix which must be ordered **1st ordinal, 2nd continuous non-mixture, 3rd components of continuous mixtures, 4th regular Poisson, 5th zero-inflated Poisson, 6th regular NB, 7th zero-inflated NB**; note that `rho` is specified in terms of the components of $Y_{mix}$ - `seed`: the seed value for random number generation (default = 1234) - `use.nearPD`: TRUE to convert `rho` to the nearest positive definite matrix with `Matrix::nearPD` if necessary - `quiet`: if FALSE prints messages, if TRUE suppresses message printing ## Returns A list with components: `rho` the target correlation matrix, which will differ from the supplied matrix (if provided) if it was converted to the nearest positive-definite matrix `L_rho` the lower correlation bound `U_rho` the upper correlation bound If continuous variables are desired, additional components are: `constants` the calculated constants `sixth_correction` a vector of the sixth cumulant correction values `valid.pdf` a vector with i-th component equal to "TRUE" if variable Y_i has a valid power method PDF, else "FALSE" If a target correlation matrix `rho` is provided, each pairwise correlation is checked to see if it is within the lower and upper bounds. If the correlation is outside the bounds, the indices of the variable pair are given. `valid.rho` TRUE if all entries of `rho` are within the bounds, else FALSE ## Reasons for Function Errors 1) The most likely cause for function errors is that no solutions to `fleish` or `poly` converged when using `find_constants`. If this happens, the function will stop. It may help to first use `find_constants` for each continuous variable to determine if a sixth cumulant correction value is needed. If the standardized cumulants are obtained from `calc_theory`, the user may need to use rounded values as inputs (i.e. `skews = round(skews, 8)`). For example, in order to ensure that skew is exactly 0 for symmetric distributions. 2) The kurtosis may be outside the region of possible values. There is an associated lower boundary for kurtosis associated with a given skew (for Fleishman's method) or skew and fifth and sixth cumulants (for Headrick's method). Use `calc_lower_skurt` to determine the boundary for a given set of cumulants. ## Examples ```r validcorr2(n = 1000, k_cat = 1, k_cont = 1, method = "Polynomial", means = 0, vars = 1, skews = 0, skurts = 0, fifths = 0, sixths = 0, marginal = list(c(1/3, 2/3)), rho = matrix(c(1, 0.4, 0.4, 1), 2, 2), quiet = TRUE) ## Not run: # 2 continuous mixture, 1 binary, 1 zero-inflated Poisson, and # 1 zero-inflated NB variable n <- 10000 seed <- 1234 # Mixture variables: Normal mixture with 2 components; # mixture of Logistic(0, 1), Chisq(4), Beta(4, 1.5) # Find cumulants of components of 2nd mixture variable L <- calc_theory("Logistic", c(0, 1)) C <- calc_theory("Chisq", 4) B <- calc_theory("Beta", c(4, 1.5)) skews <- skurts <- fifths <- sixths <- NULL Six <- list() mix_pis <- list(c(0.4, 0.6), c(0.3, 0.2, 0.5)) mix_mus <- list(c(-2, 2), c(L[1], C[1], B[1])) mix_sigmas <- list(c(1, 1), c(L[2], C[2], B[2])) mix_skews <- list(rep(0, 2), c(L[3], C[3], B[3])) mix_skurts <- list(rep(0, 2), c(L[4], C[4], B[4])) mix_fifths <- list(rep(0, 2), c(L[5], C[5], B[5])) mix_sixths <- list(rep(0, 2), c(L[6], C[6], B[6])) mix_Six <- list(list(NULL, NULL), list(1.75, NULL, 0.03)) Nstcum <- calc_mixmoments(mix_pis[[1]], mix_mus[[1]], mix_sigmas[[1]], mix_skews[[1]], mix_skurts[[1]], mix_fifths[[1]], mix_sixths[[1]]) Mstcum <- calc_mixmoments(mix_pis[[2]], mix_mus[[2]], mix_sigmas[[2]], mix_skews[[2]], mix_skurts[[2]], mix_fifths[[2]], mix_sixths[[2]]) means <- c(Nstcum[1], Mstcum[1]) vars <- c(Nstcum[2]^2, Mstcum[2]^2) marginal <- list(0.3) support <- list(c(0, 1)) lam <- 0.5 p_zip <- 0.1 pois_eps <- 0.0001 size <- 2 prob <- 0.75 p_zinb <- 0.2 nb_eps <- 0.0001 k_cat <- k_pois <- k_nb <- 1 k_cont <- 0 k_mix <- 2 Rey <- matrix(0.39, 8, 8) diag(Rey) <- 1 rownames(Rey) <- colnames(Rey) <- c("O1", "M1_1", "M1_2", "M2_1", "M2_2", "M2_3", "P1", "NB1") # set correlation between components of the same mixture variable to 0 Rey["M1_1", "M1_2"] <- Rey["M1_2", "M1_1"] <- 0 Rey["M2_1", "M2_2"] <- Rey["M2_2", "M2_1"] <- Rey["M2_1", "M2_3"] <- 0 Rey["M2_3", "M2_1"] <- Rey["M2_2", "M2_3"] <- Rey["M2_3", "M2_2"] <- 0 # check parameter inputs validpar(k_cat, k_cont, k_mix, k_pois, k_nb, "Polynomial", means, vars, skews, skurts, fifths, sixths, Six, mix_pis, mix_mus, mix_sigmas, mix_skews, mix_skurts, mix_fifths, mix_sixths, mix_Six, marginal, support, lam, p_zip, size, prob, mu = NULL, p_zinb, pois_eps, nb_eps, Rey) # check to make sure Rey is within the feasible correlation boundaries validcorr2(n, k_cat, k_cont, k_mix, k_pois, k_nb, "Polynomial", means, vars, skews, skurts, fifths, sixths, Six, mix_pis, mix_mus, mix_sigmas, mix_skews, mix_skurts, mix_fifths, mix_sixths, mix_Six, marginal, lam, p_zip, size, prob, mu = NULL, p_zinb, pois_eps, nb_eps, Rey, seed) ## End(Not run) ``` ## References Please see references for `SimCorrMix`. ## See Also `find_constants`, `corrvar2`, `validpar`

validcorr2 function