Skip to contents

Generates random data \((Y, X_1, \ldots, X_p)\) jointly multivariate normal with user-specified marginal means, marginal SDs, and full correlation structure. Useful as a backbone for sensitivity analyses, Monte Carlo studies of regression sample size methods, and pedagogical demonstrations.

Usage

simulate_regression_data(
  N,
  p,
  rho_YX,
  rho_XX = NULL,
  mu_Y = 0,
  mu_X = 0,
  sigma_Y = 1,
  sigma_X = 1,
  seed = NULL,
  column_names = NULL
)

Arguments

N

The total sample size (a positive integer \(\ge p + 2\)).

p

The number of predictor variables.

rho_YX

A numeric vector of length p giving the population correlations between \(Y\) and each predictor \(X_j\).

rho_XX

A \(p \times p\) symmetric correlation matrix for the predictors. Defaults to the identity matrix (orthogonal predictors).

mu_Y

The population mean of \(Y\) (default 0).

mu_X

A numeric vector of length p giving the population mean of each predictor (default 0, recycled across predictors).

sigma_Y

The population standard deviation of \(Y\) (default 1, in which case the simulated \(Y\) is on the standardized scale).

sigma_X

A single number or a numeric vector of length p giving the population standard deviation of each predictor (default 1, standardized).

seed

Optional integer random seed for reproducibility (default NULL).

column_names

Optional character vector of length p + 1 giving column names for the returned data.frame; defaults to c("y", "x1", "x2", ..., "xp").

Value

A data.frame with N rows and p + 1 columns: the outcome \(Y\) (first column) followed by predictors \(X_1, \ldots, X_p\).

Details

Internally the joint correlation matrix is assembled as $$R = \begin{pmatrix} 1 & \rho_{YX}^\top \\ \rho_{YX} & R_{XX} \end{pmatrix},$$ converted to a covariance matrix via the supplied SDs, and N draws are taken using mvrnorm. The resulting \(Y\) and predictors satisfy the requested marginal means and standard deviations and (in expectation) the requested correlation structure.

Author

Ken Kelley kkelley@nd.edu

Examples

# Five orthogonal predictors, each correlating .30 with Y.
set.seed(113)
d <- simulate_regression_data(
  N      = 200,
  p      = 5,
  rho_YX = rep(0.30, 5)
)
summary(lm(y ~ ., data = d))$r.squared   # ~ 5 * 0.30^2 = 0.45
#> [1] 0.5074721

# Predictors with shared structure (exchangeable correlation matrix).
rho_XX <- matrix(0.5, nrow = 5, ncol = 5); diag(rho_XX) <- 1
simulate_regression_data(
  N      = 300,
  p      = 5,
  rho_YX = c(.50, .40, .30, .20, .10),
  rho_XX = rho_XX,
  seed   = 113
)[1:3, ]
#>            y         x1         x2         x3         x4         x5
#> 1 -0.3960700 -0.3724134 -0.5445050  0.4563194  0.1107002  0.1054664
#> 2 -0.1355309 -1.4695662 -0.4833627 -2.2009360 -0.7832249 -0.7055538
#> 3 -0.6852586 -0.9169909 -0.3397163 -0.5326804  0.1738402 -1.0443534