2 added 3 characters in body

I have two (related) questions regarding the simulation of correlated, non-normal data with a user-specified value of Mardia’s 1970 multivariate kurtosis.

(1) It is very common in my field to use the method developed by Vale and Maurelli (1983) of doing a polynomial transformation of standard, normal random variables ($$Y=a + bX + cX^{2} + dX^{3}$$ where $$X \sim N(0,1)$$) where the coefficients c and d control the (univariate) skewness and kurtosis estimates. The non-normality is then achieved through the 1-dimensional marginal distributions but there is no control over the multivariate skewness/kurtosis estimates (as defined in Mardia (1970). So here comes my first question:

How can I obtain the population values of multivariate kurtosis if I have only the population values of all the univariate marginals?

Here's my attempt at a solution. It's using the semTools package in R (calculates Mardia's measure) and the lavaan package (implements the Vale & Maurelli method)

    library(semTools)
library(lavaan)

### three variables all correlated at 0.5
model <- 'x1 ~~ 0.5*x2
x2 ~~ 0.5*x3
x3 ~~ 0.5*x1'

### vector to store Mardia's kurtosis values
mardia_values <- repdouble(100)

for (i in 1:100){

### generate simulated data with N=100,000 population skewness of 2 and kurtosis of 7
dataz <- simulateData(model, sample.nobs=100000, skewness=c(2,2,2), kurtosis=c(7,7,7))

mardia_values[i] <- as.numeric(mardiaKurtosis(datz)[1])
}

mean(mardia_values)


So essentially I'm approximating it via a simulation (at a large sample size)

Would this method be acceptable? I'm very open to suggestions (particularly if it can be derived analytically

(2) Now, on a related note... does anyone know of any method to simulate data where the value of Mardia Kurtosis can be specified by the user? So far, the only thing I have been able to come up with is, once again, using the Vale and Maurelli (1983) method and sort of just go by trail and error until I get the value that I want. Since Mardia defined these measures in the population, I assume it should be possible to have some method to generate data where the value can be specified in advance, but I do not know how.

Thank you!

References

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519-530.

Vale, C. D., & Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.

I have two (related) questions regarding the simulation of correlated, non-normal data with a user-specified value of Mardia’s 1970 multivariate kurtosis.

(1) It is very common in my field to use the method developed by Vale and Maurelli (1983) of doing a polynomial transformation of standard, normal random variables ($$Y=a + bX + cX^{2} + dX^{3}$$ where $$X \sim N(0,1)$$) where the coefficients c and d control the (univariate) skewness and kurtosis estimates. The non-normality is then achieved through the 1-dimensional marginal distributions but there is no control over the multivariate skewness/kurtosis estimates (as defined in Mardia (1970). So here comes my first question:

How can I obtain the population values of multivariate kurtosis if I have only the population values of all the univariate marginals?

Here's my attempt at a solution. It's using the semTools package in R (calculates Mardia's measure) and the lavaan package (implements the Vale & Maurelli method)

    library(semTools)
library(lavaan)

### three variables all correlated at 0.5
model <- 'x1 ~~ 0.5*x2
x2 ~~ 0.5*x3
x3 ~~ 0.5*x1'

### vector to store Mardia's kurtosis values
mardia_values <- rep(100)

for (i in 1:100){

### generate simulated data with N=100,000 population skewness of 2 and kurtosis of 7
dataz <- simulateData(model, sample.nobs=100000, skewness=c(2,2,2), kurtosis=c(7,7,7))

mardia_values[i] <- as.numeric(mardiaKurtosis(datz)[1])
}

mean(mardia_values)


So essentially I'm approximating it via a simulation (at a large sample size)

Would this method be acceptable? I'm very open to suggestions (particularly if it can be derived analytically

(2) Now, on a related note... does anyone know of any method to simulate data where the value of Mardia Kurtosis can be specified by the user? So far, the only thing I have been able to come up with is, once again, using the Vale and Maurelli (1983) method and sort of just go by trail and error until I get the value that I want. Since Mardia defined these measures in the population, I assume it should be possible to have some method to generate data where the value can be specified in advance, but I do not know how.

Thank you!

References

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519-530.

Vale, C. D., & Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.

I have two (related) questions regarding the simulation of correlated, non-normal data with a user-specified value of Mardia’s 1970 multivariate kurtosis.

(1) It is very common in my field to use the method developed by Vale and Maurelli (1983) of doing a polynomial transformation of standard, normal random variables ($$Y=a + bX + cX^{2} + dX^{3}$$ where $$X \sim N(0,1)$$) where the coefficients c and d control the (univariate) skewness and kurtosis estimates. The non-normality is then achieved through the 1-dimensional marginal distributions but there is no control over the multivariate skewness/kurtosis estimates (as defined in Mardia (1970). So here comes my first question:

How can I obtain the population values of multivariate kurtosis if I have only the population values of all the univariate marginals?

Here's my attempt at a solution. It's using the semTools package in R (calculates Mardia's measure) and the lavaan package (implements the Vale & Maurelli method)

    library(semTools)
library(lavaan)

### three variables all correlated at 0.5
model <- 'x1 ~~ 0.5*x2
x2 ~~ 0.5*x3
x3 ~~ 0.5*x1'

### vector to store Mardia's kurtosis values
mardia_values <- double(100)

for (i in 1:100){

### generate simulated data with N=100,000 population skewness of 2 and kurtosis of 7
dataz <- simulateData(model, sample.nobs=100000, skewness=c(2,2,2), kurtosis=c(7,7,7))

mardia_values[i] <- as.numeric(mardiaKurtosis(datz)[1])
}

mean(mardia_values)


So essentially I'm approximating it via a simulation (at a large sample size)

Would this method be acceptable? I'm very open to suggestions (particularly if it can be derived analytically

(2) Now, on a related note... does anyone know of any method to simulate data where the value of Mardia Kurtosis can be specified by the user? So far, the only thing I have been able to come up with is, once again, using the Vale and Maurelli (1983) method and sort of just go by trail and error until I get the value that I want. Since Mardia defined these measures in the population, I assume it should be possible to have some method to generate data where the value can be specified in advance, but I do not know how.

Thank you!

References

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519-530.

Vale, C. D., & Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.

1

# How to simulate data with values of Mardia's Kurtosis?

I have two (related) questions regarding the simulation of correlated, non-normal data with a user-specified value of Mardia’s 1970 multivariate kurtosis.

(1) It is very common in my field to use the method developed by Vale and Maurelli (1983) of doing a polynomial transformation of standard, normal random variables ($$Y=a + bX + cX^{2} + dX^{3}$$ where $$X \sim N(0,1)$$) where the coefficients c and d control the (univariate) skewness and kurtosis estimates. The non-normality is then achieved through the 1-dimensional marginal distributions but there is no control over the multivariate skewness/kurtosis estimates (as defined in Mardia (1970). So here comes my first question:

How can I obtain the population values of multivariate kurtosis if I have only the population values of all the univariate marginals?

Here's my attempt at a solution. It's using the semTools package in R (calculates Mardia's measure) and the lavaan package (implements the Vale & Maurelli method)

    library(semTools)
library(lavaan)

### three variables all correlated at 0.5
model <- 'x1 ~~ 0.5*x2
x2 ~~ 0.5*x3
x3 ~~ 0.5*x1'

### vector to store Mardia's kurtosis values
mardia_values <- rep(100)

for (i in 1:100){

### generate simulated data with N=100,000 population skewness of 2 and kurtosis of 7
dataz <- simulateData(model, sample.nobs=100000, skewness=c(2,2,2), kurtosis=c(7,7,7))

mardia_values[i] <- as.numeric(mardiaKurtosis(datz)[1])
}

mean(mardia_values)


So essentially I'm approximating it via a simulation (at a large sample size)

Would this method be acceptable? I'm very open to suggestions (particularly if it can be derived analytically

(2) Now, on a related note... does anyone know of any method to simulate data where the value of Mardia Kurtosis can be specified by the user? So far, the only thing I have been able to come up with is, once again, using the Vale and Maurelli (1983) method and sort of just go by trail and error until I get the value that I want. Since Mardia defined these measures in the population, I assume it should be possible to have some method to generate data where the value can be specified in advance, but I do not know how.

Thank you!

References

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57(3), 519-530.

Vale, C. D., & Maurelli, V. A. (1983). Simulating multivariate nonnormal distributions. Psychometrika, 48(3), 465-471.