Simulate correlated general-distribution variables with a Gaussian copula

Usage

simulate_copula(
  n,
  specification,
  rho = 0,
  tau = NULL,
  structure = c("independent", "exchangeable", "ar1", "custom"),
  correlation = NULL,
  id = "id",
  seed = NULL,
  envir = parent.frame(),
  batch = NULL
)

Arguments

n: Number of observations.
specification: Variable definitions from define_variables(). In the distribution-call form every marginal must carry a quantile function; list_distributions(copula = TRUE) reports which do. The formula/variance column form covers twelve marginals: formula is the mean (through link) and variance carries the second parameter, which is a variance for normal, a precision for beta, a dispersion for gamma and negative_binomial, and a size for binomial. It is unused by binary, exponential, poisson and no_zero_poisson; for uniform and uniform_integer formula instead carries both limits, and for categorical formula carries the category probabilities and variance their labels. A marginal outside those twelve raises simulab_no_quantile; state it as a distribution call instead.
rho: Optional latent Pearson correlation, zero by default.
tau: Optional Kendall correlation, overriding rho. It is converted to the latent Pearson correlation sin(pi * tau / 2), so a continuous margin pair attains Kendall's tau of tau.
structure: Correlation structure: one of "independent", "exchangeable", "ar1", or "custom". If left unset it is chosen from the other arguments: "custom" when correlation is supplied, "exchangeable" when a non-zero rho (or tau) is supplied, and "independent" otherwise. Passing "independent" together with a non-zero rho is a contradiction and raises an error.
correlation: Custom latent correlation matrix or tidy table.
id: Identifier name.
seed: Optional random seed.
envir: Formula evaluation environment.
batch: Optional single positive whole number. When given, the simulator runs batch times and returns a plain list of batch results, each exactly what the same call without batch returns. With a seed, every dataset gets its own seed drawn from seed, so the whole batch is reproducible; without one the datasets are consecutive draws from the session's random-number stream. The default NULL returns a single result.

Value

A simulab_sim base data.frame with one row per observation and columns id (renamed by id) and one column per variable. Two tidy tables come from as.data.frame(): what = "definitions" restates the specification, and what = "latent_correlation" has one row per correlation-matrix cell with columns row, column and correlation. That matrix is the correlation of the latent Gaussian variables, not of the generated margins, which are pulled apart by their own quantile transforms.

Conditions

simulab_no_quantile when a distribution call names a marginal with no quantile function, such as rice() or mixture().

Examples

# Distribution calls. Each marginal keeps its own distribution while the
# variables are correlated through a Gaussian copula.
result <- simulate_copula(
  n = 200,
  specification = define_variables(
    score  = normal(mean = 10, sd = 2),
    visits = poisson(lambda = 4),
    rate   = beta(shape1 = 2, shape2 = 5)
  ),
  rho = 0.5, seed = 1
)
head(result)
#> <simulab_sim:copula> 6 rows x 4 columns
#>   id     score visits      rate
#> 1  1  9.518238      5 0.4463908
#> 2  2 12.036044      9 0.6890614
#> 3  3  8.887993      6 0.2012495
#> 4  4 12.667842      4 0.2524626
#> 5  5  9.347847      0 0.1384510
#> 6  6  9.453232      9 0.2715592
#> 
#> Truth (definitions):
#>   variable distribution parameter value
#> 1    score       normal      mean    10
#> 2    score       normal        sd     2
#> 3   visits      poisson    lambda     4
#> 4     rate         beta    shape1     2
#> 5     rate         beta    shape2     5
#> 
#> Other tables: latent_correlation. Read one with as.data.frame(x, what = "latent_correlation").

specification <- define_variables(
  define_variable("normal_var", formula = "0", variance = "1", distribution = "normal"),
  define_variable("count_var", formula = "3", variance = "0", distribution = "poisson")
)
result <- simulate_copula(n = 200, specification = specification, rho = 0.5, seed = 1)
head(result)
#> <simulab_sim:copula> 6 rows x 3 columns
#>   id normal_var count_var
#> 1  1 -0.4991469         3
#> 2  2  0.6144984         6
#> 3  3 -0.3965160         5
#> 4  4  1.4552777         3
#> 5  5 -0.2731824         0
#> 6  6 -0.1460692         7
#> 
#> Truth (definitions):
#>     variable distribution formula variance     link
#> 1 normal_var       normal       0        1 identity
#> 2  count_var      poisson       3        0 identity
#> 
#> Other tables: latent_correlation. Read one with as.data.frame(x, what = "latent_correlation").
components(result)
#>                table rows columns
#> 1               data  200       3
#> 2        definitions    2       5
#> 3 latent_correlation    4       3