simulate_hierarchical_factors: Simulates Hierarchical Latent Factor Data

Description

Simulates data from a hierarchical (second-order) latent factor model based on many manipulable parameters. Variables load onto lower-order factors, and lower-order factors load onto higher-order factors. Parameters do not have default values and must each be set, with the exception of off_disturbances and lower_correlations (and their _range variants), which are mutually exclusive alternatives for specifying the same underlying structure (exactly one must be supplied). See examples to get started

The population correlation matrix (variable-level) implied by the hierarchical factor model is:

$$\boldsymbol{\Sigma} = \boldsymbol{\Lambda}_L \boldsymbol{\Phi}_L \boldsymbol{\Lambda}_L' + \boldsymbol{\Theta}$$

where $\boldsymbol{\Lambda}_L$ is the lower-order loading matrix (variables x lower-order factors; lower_loadings and lower_cross_loadings), $\boldsymbol{\Phi}_L$ is the implied lower-order factor correlation matrix, and $\boldsymbol{\Theta}$ is a diagonal matrix of uniquenesses that ensures unit variance for each variable

$\boldsymbol{\Phi}_L$ is, in turn, implied by the higher-order factor structure:

$$\boldsymbol{\Phi}_L = \boldsymbol{\Lambda}_H \boldsymbol{\Phi}_H \boldsymbol{\Lambda}_H' + \boldsymbol{\Psi}$$

where $\boldsymbol{\Lambda}_H$ is the higher-order loading matrix (lower-order factors x higher-order factors; higher_loadings and higher_cross_loadings), $\boldsymbol{\Phi}_H$ is the higher-order factor correlation matrix (higher_correlations), and $\boldsymbol{\Psi}$ is the disturbance (residual) correlation matrix of the lower-order factors, with its off-diagonal (correlated disturbance) elements set directly via off_disturbances — its diagonal is not controllable: whatever variance $\boldsymbol{\Phi}_H$ and $\boldsymbol{\Lambda}_H$ leave unexplained in each lower-order factor is always rescaled back to unit variance, so only the correlation between disturbances (not their magnitude) can be set

Usage

simulate_hierarchical_factors(
  lower_factors,
  variables,
  variables_range = NULL,
  lower_loadings,
  lower_loadings_range = NULL,
  lower_cross_loadings,
  lower_cross_loadings_range = NULL,
  lower_correlations = NULL,
  lower_correlations_range = NULL,
  higher_factors,
  higher_variables = NULL,
  higher_loadings,
  higher_loadings_range = NULL,
  higher_cross_loadings,
  higher_cross_loadings_range = NULL,
  higher_correlations = NULL,
  higher_correlations_range = NULL,
  off_disturbances = NULL,
  off_disturbances_range = NULL,
  sample_size,
  variable_categories = Inf,
  categorical_limit = 7,
  skew = 0,
  skew_range = NULL
)

Value

Returns a list containing:

data

Simulated data from the specified hierarchical factor model

population_correlation

Population correlation matrix (variable-level)

lower_correlation

Population correlation matrix between lower-order factors, implied by the higher-order factor structure and disturbances (or, when lower_correlations is supplied, equal to that target matrix)

parameters

A list containing the parameters used to generate the data:

lower --- A list of lower-order parameters: factors (number of lower-order factors), variables (variables on each lower-order factor), loadings (loading matrix of variables onto lower-order factors), cross_loadings (cross-loadings used for each lower-order factor), and correlations (population correlation matrix between lower-order factors; same matrix as the top-level lower_correlation)
higher --- A list of higher-order parameters: factors (number of higher-order factors), variables (lower-order factors nested within each higher-order factor), loadings (loading matrix of lower-order factors onto higher-order factors), cross_loadings (cross-loadings used for each higher-order factor), correlations (correlations between higher-order factors), and disturbances (disturbance (residual) correlation matrix of lower-order factors; when off_disturbances is supplied directly, only its off-diagonal elements are meaningful, since its diagonal is always reset to 1; when implied instead, via lower_correlations, the diagonal reflects the residual variance actually solved for)
categories --- Categories for each variable
categorical_limit --- Category limit used to determine continuous variables
skew --- Skew for each variable

Arguments

lower_factors: Numeric (length = 1). Number of lower-order (first-order) factors
variables: Numeric (length = 1 or lower_factors). Number of variables per lower-order factor. Can be a single value or as many values as there are lower-order factors. Minimum two variables per factor
variables_range: Numeric (length = 2). Range of variables to randomly select from a random uniform distribution. Minimum three variables per factor
lower_loadings: Numeric or matrix (length = 1, lower_factors, or total variables x lower_factors). Loadings of variables onto their lower-order factor, drawn from a random uniform distribution using +/- 0.10 of value input. Can be a single value or as many values as there are lower-order factors (corresponding to the factors); a vector as long as the total number of variables is not applied per variable unless generated automatically via lower_loadings_range. Can also be a loading matrix. Columns must match lower_factors and rows must match total variables (lower_factors x variables) General effect sizes range from small (0.40), moderate (0.55), to large (0.70)
lower_loadings_range: Numeric (length = 2). Range of lower-order loadings to randomly select from a random uniform distribution. General effect sizes range from small (0.40), moderate (0.55), to large (0.70)
lower_cross_loadings: Numeric (length = 1 or lower_factors). Cross-loadings of variables onto non-dominant lower-order factors, drawn from a random normal distribution with a mean of 0 and standard deviation of value input. Can be a single value or as many values as there are lower-order factors (corresponding to the factors). Unlike lower_loadings, a full loading matrix is not supported: if a matrix is supplied, only its first lower_factors elements (taken column-wise) are used as the per-factor standard deviations
lower_cross_loadings_range: Numeric (length = 2). Range of lower-order cross-loadings to randomly select from a random uniform distribution
lower_correlations: Numeric or matrix (length = 1 or lower_factors x lower_factors). Target correlation matrix between lower-order factors, supplied directly as an alternative to off_disturbances. Can be a single value (applied to all off-diagonal correlations) or a full lower_factors x lower_factors correlation matrix. When supplied, off_disturbances (and, if higher_correlations is also omitted, higher_correlations too) is implied rather than set directly (see Details). Mutually exclusive with off_disturbances/off_disturbances_range: exactly one of off_disturbances/off_disturbances_range or lower_correlations/lower_correlations_range must be supplied
lower_correlations_range: Numeric (length = 2). Range of lower-order correlations to randomly select (per pair of lower-order factors) from a random uniform distribution. Somewhat redundant with lower_correlations but more flexible. Mutually exclusive with off_disturbances/off_disturbances_range
higher_factors: Numeric (length = 1). Number of higher-order (second-order) factors. lower_factors must be evenly divisible by higher_factors unless higher_loadings is supplied as a matrix or higher_variables is supplied
higher_variables: Numeric (length = 1 or higher_factors). Number of lower-order factors nested within each higher-order factor. Can be a single value (lower-order factors are split evenly across higher-order factors, equivalent to leaving higher_variables as NULL) or as many values as there are higher-order factors (allowing a different number of lower-order factors to load onto each higher-order factor). Values must sum to lower_factors. Ignored when higher_loadings is supplied as a matrix
higher_loadings: Numeric or matrix (length = 1, higher_factors, or lower_factors x higher_factors). Loadings of lower-order factors onto their higher-order factor, drawn from a random uniform distribution using +/- 0.10 of value input. Can be a single value or as many values as there are higher-order factors (corresponding to the factors); a vector as long as lower_factors is not applied per lower-order factor unless generated automatically via higher_loadings_range. Can also be a loading matrix. Columns must match higher_factors and rows must match lower_factors General effect sizes range from small (0.40), moderate (0.55), to large (0.70)
higher_loadings_range: Numeric (length = 2). Range of higher-order loadings to randomly select from a random uniform distribution. General effect sizes range from small (0.40), moderate (0.55), to large (0.70)
higher_cross_loadings: Numeric (length = 1 or higher_factors). Cross-loadings of lower-order factors onto non-dominant higher-order factors, drawn from a random normal distribution with a mean of 0 and standard deviation of value input. Can be a single value or as many values as there are higher-order factors (corresponding to the factors). Unlike higher_loadings, a full loading matrix is not supported: if a matrix is supplied, only its first higher_factors elements (taken column-wise) are used as the per-factor standard deviations
higher_cross_loadings_range: Numeric (length = 2). Range of higher-order cross-loadings to randomly select from a random uniform distribution
higher_correlations: Numeric (length = 1 or higher_factors x higher_factors). Can be a single value that will be used for all correlations between higher-order factors. Can also be a square matrix (higher_factors x higher_factors). Negative values are currently converted to their absolute value (only positive correlations between higher-order factors are supported). General effect sizes range from orthogonal (0.00), small (0.30), moderate (0.50), to large (0.70). Required when off_disturbances/off_disturbances_range is used. Optional when lower_correlations/lower_correlations_range is used instead, in which case omitting it (and higher_correlations_range) implies it instead of supplying it directly (see Details)
higher_correlations_range: Numeric (length = 2). Range of correlations between higher-order factors to randomly select from a random uniform distribution. General effect sizes range from orthogonal (0.00), small (0.30), moderate (0.50), to large (0.70). Required when off_disturbances/off_disturbances_range is used, unless higher_correlations is supplied instead. Optional when lower_correlations/lower_correlations_range is used (see higher_correlations)
off_disturbances: Numeric or matrix (length = 1 or lower_factors x lower_factors). Off-diagonal (correlated) disturbances between lower-order factors' residuals, left unexplained by the higher-order factor(s); the diagonal is not settable here (see Details for why). A single value is applied uniformly (no jitter) to every off-diagonal element; a full lower_factors x lower_factors matrix is used exactly as supplied instead (its diagonal is ignored and reset to 1). Values must be between -1 and 1 and the resulting matrix (with unit diagonal) must be positive semi-definite. Mutually exclusive with lower_correlations/lower_correlations_range: exactly one of off_disturbances/off_disturbances_range or lower_correlations/lower_correlations_range must be supplied
off_disturbances_range: Numeric (length = 2). Range of off-diagonal (correlated) disturbances to randomly and independently select (per pair of lower-order factors) from a random uniform distribution. Mutually exclusive with lower_correlations/lower_correlations_range
sample_size: Numeric (length = 1). Number of cases to generate from a random multivariate normal distribution using rmvnorm
variable_categories: Numeric (length = 1 or total variables (lower_factors x variables)). Number of categories for each variable. Inf is used for continuous variables; otherwise, values reflect number of categories
categorical_limit: Numeric (length = 1). Values greater than input value are considered continuous. Defaults to 7 meaning that 8 or more categories are considered continuous (i.e., variables are not categorized from continuous to categorical)
skew: Numeric (length = 1 or categorical variables). Skew to be included in categorical variables. It is randomly sampled from provided values. Can be a single value or as many values as there are (total) variables. Current skew implementation is between -10 and 10 in increments of 0.05. Skews that are not in this sequence will be converted to their nearest value in the sequence. Not recommended to use with variables_range
skew_range: Numeric (length = 2). Randomly selects skews within in the range. Somewhat redundant with skew but more flexible. Compatible with variables_range

Author

Alexander P. Christensen <alexpaulchristensen@gmail.com>

Details

How $\boldsymbol{\Phi}_H$ and $\boldsymbol{\Psi}$ are obtained depends on which parameters are supplied. Exactly one of off_disturbances/off_disturbances_range or lower_correlations/lower_correlations_range must be supplied; higher_correlations/higher_correlations_range is required alongside the former, but optional alongside the latter:

Supplied	Implied
`higher_correlations` + `off_disturbances`	`lower_correlations` (i.e., $\boldsymbol{\Phi}_L$)
`higher_correlations` + `lower_correlations`	`off_disturbances` (i.e., $\boldsymbol{\Psi}$)
`lower_correlations` (alone)	`higher_correlations` and `off_disturbances` (i.e., $\boldsymbol{\Phi}_H$ and $\boldsymbol{\Psi}$)

When $\boldsymbol{\Psi}$ is implied, it is solved for as the residual needed to reproduce the target lower_correlations matrix given $\boldsymbol{\Phi}_H$ and $\boldsymbol{\Lambda}_H$, and may end up with a non-unit diagonal and/or off-diagonal values (i.e., correlated disturbances between lower-order factors are permitted); unlike a directly supplied off_disturbances, this implied $\boldsymbol{\Psi}$ keeps whatever diagonal falls out of the subtraction, since it is not passed back through the unit-variance rescaling that a supplied off_disturbances is. When $\boldsymbol{\Phi}_H$ is also implied (i.e., higher_correlations is omitted), it is obtained first, from the target lower_correlations matrix and $\boldsymbol{\Lambda}_H$, via a generalized-inverse (least-squares) regression:

$$\boldsymbol{\Phi}_H = (\boldsymbol{\Lambda}_H'\boldsymbol{\Lambda}_H)^{-1}\boldsymbol{\Lambda}_H' \boldsymbol{\Phi}_L \boldsymbol{\Lambda}_H(\boldsymbol{\Lambda}_H'\boldsymbol{\Lambda}_H)^{-1}$$

This raw least-squares result is not guaranteed to be a valid correlation matrix (unit diagonal, [-1, 1] off-diagonals), so it is standardized with cov2cor, after which its diagonal is reset to 1 and off-diagonal elements are taken in absolute value (consistent with how a directly supplied higher_correlations is handled); $\boldsymbol{\Psi}$ is then solved for as described above

In every case, $\boldsymbol{\Psi}$ (whether fixed or implied) must form a valid (positive semi-definite) covariance matrix; if a target correlation matrix or fixed higher-order structure is incompatible (e.g., it implies negative unique variances), values are redrawn when higher_loadings_range, higher_correlations_range, and/or higher_cross_loadings_range are used, or an error is raised when the higher-order structure is entirely fixed

Examples

Run this code

# Generate hierarchical factor data
# 4 lower-order factors (2 per higher-order factor), 2 higher-order factors
hierarchical <- simulate_hierarchical_factors(
  lower_factors = 4, # lower-order factors = 4
  variables = 6, # variables per lower-order factor = 6
  lower_loadings = 0.55, # lower-order loadings = 0.45 to 0.65
  lower_cross_loadings = 0.05, # lower-order cross-loadings N(0, 0.05)
  higher_factors = 2, # higher-order factors = 2
  higher_loadings = 0.60, # higher-order loadings = 0.50 to 0.70
  higher_cross_loadings = 0.05, # higher-order cross-loadings N(0, 0.05)
  higher_correlations = 0.30, # correlation between higher-order factors = 0.30
  off_disturbances = 0.10, # off-diagonal (correlated) disturbances = 0.10
  sample_size = 1000 # number of cases = 1000
)

# Different number of lower-order factors per higher-order factor
# 6 lower-order factors (2 on higher-order factor A, 4 on higher-order factor B)
hierarchical_uneven <- simulate_hierarchical_factors(
  lower_factors = 6, variables = 6,
  lower_loadings = 0.55, lower_cross_loadings = 0.05,
  higher_factors = 2, higher_variables = c(2, 4), # A = 2 factors, B = 4 factors
  higher_loadings = 0.60, higher_cross_loadings = 0.05,
  higher_correlations = 0.30, off_disturbances = 0.10,
  sample_size = 1000
)

# Randomly vary higher-order loadings
hierarchical_loadings <- simulate_hierarchical_factors(
  lower_factors = 4, variables = 6,
  lower_loadings = 0.55, lower_cross_loadings = 0.05,
  higher_factors = 2,
  higher_loadings_range = c(0.40, 0.80), # higher-order loadings = 0.40 to 0.80
  higher_cross_loadings = 0.05,
  higher_correlations = 0.30,
  off_disturbances_range = c(0.00, 0.20), # off-diagonal disturbances = 0.00 to 0.20
  sample_size = 1000
)

# Generate dichotomous data
hierarchical_dichotomous <- simulate_hierarchical_factors(
  lower_factors = 4, variables = 6,
  lower_loadings = 0.55, lower_cross_loadings = 0.05,
  higher_factors = 2, higher_loadings = 0.60,
  higher_cross_loadings = 0.05, higher_correlations = 0.30,
  off_disturbances = 0.10, sample_size = 1000,
  variable_categories = 2 # dichotomous data
)

# Supply a target lower-order correlation matrix directly instead of `off_disturbances`;
# the disturbances are implied (solved for) and may end up correlated
hierarchical_lower_correlations <- simulate_hierarchical_factors(
  lower_factors = 4, variables = 6,
  lower_loadings = 0.55, lower_cross_loadings = 0.05,
  lower_correlations = 0.40, # target correlation of 0.40 between all lower-order factors
  higher_factors = 2, higher_loadings = 0.60,
  higher_cross_loadings = 0.05, higher_correlations = 0.30,
  sample_size = 1000
)

# Omit `higher_correlations` when `lower_correlations` is supplied: it is implied
# (rather than the disturbances) from the target lower-order correlations and `higher_loadings`
hierarchical_implied_correlations <- simulate_hierarchical_factors(
  lower_factors = 4, variables = 6,
  lower_loadings = 0.55, lower_cross_loadings = 0.05,
  lower_correlations = 0.40, # target correlation of 0.40 between all lower-order factors
  higher_factors = 2, higher_loadings = 0.60,
  higher_cross_loadings = 0.05,
  sample_size = 1000
)

Run the code above in your browser using DataLab