5
$\begingroup$

I am trying to use R to simulate random variations in a real dataset with a known number of categorical and continuous predictor variables, as well as known marginal probabilities for each ordinal/categorical variable. These parameters are as follows:

ncont <- 4  # number of continuous predictor variables
ncat  <- 8  # number of categorical predictor variables

margprobs <- list(c(0.24, 0.3, 0.28),
                  c(0.01, 0.071, 0.17, 0.069, 0.005, 0.05),
                  c(0.07, 0.06, 0.14, 0.3),
                  c(0.05, 0.45, 0.5),
                  c(0.9),
                  c(0.05, 0.8),
                  c(0.8),
                  c(0.5))

I can simulate this dataset if I assume that the categorical variables are independent, and I can create valid correlation matrices for just the continuous variables. However, every time I try to create a simulated variance-covariance or correlation matrix for the categorical variables, it is not valid in the context of the known marginal probabilities.

I have consulted the vignettes for the OrdNor, SimCorrMix, and SimMultiCorrData packages for how to produce a valid correlation or variance-covariance matrix with known marginal probabilities, but I can't figure it out. I can use these packages to find the upper and lower bounds of acceptable correlations in the context of the marginal probabilities, but every correlation matrix I sample passes the validation check these packages provide but ultimately fails when being converted to an intermediate correlation matrix. I provide some example code for this below:

library(SimCorrMix)
margprobs <- lapply(margprobs, FUN=sort)  # must sort margprobs for corrvar() to work
corlims   <- validcorr(1000, k_cat=8, marginal=margprobs)

set.seed(1)
Rey <- diag(8)
for (i in 1:nrow(Rey)) {
  for (j in 1:ncol(Rey)) {
    if (i > j) Rey[i, j] <- runif(1, corlims$L_rho[i, j], corlims$U_rho[i, j])
    Rey[j, i] <- Rey[i, j]
  }
}

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

# Use SimCorrMix to check for valid correlation matrix
validcorr(1000, k_cat=8, marginal=margprobs, rho=Rey)
# TRUE

ords1 <- corrvar(1000, k_cat=8, marginal=margprobs, rho=Rey)
# Error message: It is not possible to find a correlation matrix for MVN 
# ensuring rho for the ordinal variables.  Try the error loop.
# Intermediate correlation matrix is not positive definite. Nearest 
# positive definite matrix is used.

# Note: this error happens whether or not I use Matrix::nearPD or 
# psych::cor.smooth to ensure that the rho matrix is positive definite

Please help me understand how to create a valid ordinal correlation matrix that will work with these simulations.

I want to eventually produce a simulation of continuous and categorical variables based on a valid correlation matrix for all variables, and then I can see how the magnitudes of these correlations affect the ability of models to recover known coefficients (beta) and outcome values (y) that I provide in the simulations.

$\endgroup$
0
$\begingroup$

According to the manual page of corrvar, "marginal" should be a vector of cumulative probabilities excluding the last (which is always = 1). So instead of c(0.24, 0.3, 0.28), for example, try cumsum(c(0.24, 0.3, 0.28)).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.