---
title: "Simulating a Wide Variety of Data Shapes, Sizes and Distributions"
author:
  - name: "[Mohammed Saqr](https://saqr.me)"
    affiliation: University of Eastern Finland
  - name: "[Sonsoles López-Pernas](https://sonsoles.me)"
    affiliation: University of Eastern Finland
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simulating a Wide Variety of Data Shapes, Sizes and Distributions}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_igraph <- requireNamespace("igraph", quietly = TRUE)
has_tna <- requireNamespace("tna", quietly = TRUE)
has_survival <- requireNamespace("survival", quietly = TRUE)
fmt <- function(x, digits = 2) sub("^-", "\u2212", formatC(x, format = "f", digits = digits))
pct <- function(x, digits = 1) paste0(fmt(100 * x, digits), "%")
join_and <- function(v) sub(", ([^,]*)$", " and \\1", toString(v))
statistic_of <- function(summary, term, variable, statistic) {
  summary[summary$term == term & summary$variable == variable, statistic]
}
coefficient_of <- function(fit, term) unname(stats::coef(fit)[term])
off_diagonal <- function(m, lag) {
  rows <- seq_len(nrow(m) - lag)
  m[cbind(rows, rows + lag)]
}
upper_values <- function(m) m[upper.tri(m)]
pair_value <- function(m, a, b) m[a, b]
aov_sums <- function(fit) summary(fit)[[1]][["Sum Sq"]]
n_distributions <- nrow(simulab::list_distributions())
n_learning_states <- nrow(simulab::learning_states())
n_state_categories <- length(simulab::learning_state_categories())
n_name_regions <- length(simulab::global_name_regions())
```

simulab generates a wide variety of data, in different shapes, sample sizes,
distributions and parameters. The parameters are either drawn at random or
set in advance, and in both cases they are stored with the data, so they are
known, recoverable and testable. It is designed for simulation, for the
evaluation of statistical methods, for training and testing models and
software, and for learning and teaching.

Real data are often scarce, and the scarcity is greatest for learners, who
rarely have access to data sets of the size, structure and design that a
statistical method requires. Real data about people also face privacy and
ethical constraints, and access to them is often subject to consent, ethical
review, data access agreements, and limits on storage and redistribution.
Data generated by simulab from specified parameters describe no real persons,
so they can be shared, published, modified and given to a class without these
restrictions. Their main value lies in the control and breadth they provide. A
real data set has one sample size, one pattern of missing values, one set of
distributions and one structure of relationships. simulab generates data that
cover small and large samples; complete and incomplete data; normal, skewed,
heavy-tailed, bounded, discrete and zero-inflated distributions; and
relationships that are weak, strong, linear, nonlinear or absent. It can also
add outliers, measurement error, unequal variances, dependence and other
conditions that affect statistical methods. This breadth gives learners
examples that show how statistical methods behave under different conditions.

Data of different types and sizes serve learning, teaching and assessment.
Small samples illustrate sampling variation, unstable estimates, low power and
wide confidence intervals, and large samples show greater precision and the
distinction between statistical and practical significance. Data with
different distributions, missing-data mechanisms, outliers and dependence show
why assumptions matter and when alternative methods are needed. Data with known
population parameters let learners compare estimates with true values and
examine bias, variance, coverage, power and error rates. simulab can give each
learner a different data set from the same model, so the answers differ between
learners while the method being assessed is the same. The same models support a
progression from simple examples to complex and realistic cases. Real data
provide context and authenticity, and simulated data provide control,
reproducibility and systematic variation.

Simulated data also serve the evaluation of newly developed models and methods.
A new method requires tests under controlled conditions that cover the range of
situations in which it may be used. Data with known population parameters let
researchers compare model estimates with true values and assess bias, variance,
coverage, power, prediction error, calibration, robustness and computational
performance. Data can also expose the weaknesses of a method under specific
conditions, such as small samples, high dimensionality, missing values, class
imbalance, skewed distributions, outliers or model misspecification. simulab
produces as many data sets as an evaluation requires, in batches from the same
model, each data set reproducible from its own seed. Repeated simulation gives
evidence about the consistency and reliability of a method across many
samples, where one empirical data set gives a single result. Real data then
provide an additional test of performance under practical conditions. A
collection of diverse data sets is a common resource for method development,
validation, comparison, teaching and assessment.

simulab groups its generators into families by the kind of data they produce.
General data are declared variable by variable from named distributions and
formulas, or drawn as correlated variables with normal, non-normal or ordinal
margins. Statistical designs cover group comparisons for t-tests and analysis
of variance, linear regression, prediction data with continuous and
categorical predictors, and Gaussian clusters. Latent variable data come from
factor, latent profile, two-level latent profile and latent class models, and
measurement data from item response models. Longitudinal data come from
multilevel models, growth curves and multivariate vector autoregressions.
Survival data have event times under proportional hazards, with censoring and
competing risks. Sequence data come from Markov chains, hidden Markov models,
mixtures of transition matrices and groups with their own transition matrices,
and behavioral data record timed events of actors in groups and courses.
Network data are random, block, preferential-attachment and small-world graphs,
temporal, bipartite and multiplex networks, and transition networks. Empirical
data are resampled from the rows of an observed data set, and so reproduce its
records, or drawn from its kernel density estimate. Functional data add a
variable generated from a spline curve. Missing values, treatment assignment
and the calibration of parameters to target quantities apply across the
families.

Each result is a tidy data frame, one row per observation, event or edge, and
it can be passed directly to modeling functions such as `lm()` and `glm()`.
A simulated transition network is returned as a native `tna` model instead.
The generating parameters, the design and, where they are defined, the
population quantities are stored with the result as tidy data frames. Matrix
inputs, among them transition matrices, factor loadings and profile means, can
also be given as tidy data frames, and inputs that describe the same
variables, states or profiles are matched by their names.

```{r library}
library(simulab)
```

## Linear regression

In the linear regression model, the outcome of observation $i$ is

$$y_i = \beta_0 + \mathbf{x}_i^\top\boldsymbol\beta + \varepsilon_i,$$

where the predictors $\mathbf{x}_i$ are multivariate normal and the errors
$\varepsilon_i$ are independent and normal with variance $\sigma^2$. Called
with no arguments, `simulate_regression()` draws the parameters itself. The
intercept is drawn from the standard normal distribution and the slopes of two
predictors from the uniform distribution on $[-1, 1]$, and 100 observations
are generated with standard normal predictors and an error standard deviation
of 1. The drawn values are printed under the data as the truth, so the data
have known parameters although none were given.

```{r regression-drawn}
drawn <- simulate_regression(seed = 1)
drawn
```

When the parameters matter to the question, they are given. In
`coefficients`, the names are the terms and the values the generating
coefficients, here an intercept of 1 and slopes of 0.5 for `x1` and −0.3 for
`x2`. The printed truth is then exactly these values.

```{r regression}
known <- simulate_regression(
  n = 100,
  coefficients = c("(Intercept)" = 1, x1 = 0.5, x2 = -0.3),
  seed = 1
)
known
```

A linear model fitted to `known` estimates the coefficients from the data
alone. In `validate_recovery()`, each estimate is set against the value that
generated it. The bias is the estimate minus the generating value, the
absolute error is its size, and `recovered` is `TRUE` when the absolute error
is at most 0.1.

```{r regression-recovery}
fit <- lm(outcome ~ x1 + x2, data = known)
validate_recovery(fit, known)
```

```{r regression-recovery-values, include = FALSE}
single_recovery <- validate_recovery(fit, known)
largest_single <- max(single_recovery$absolute_error)
largest_term <- subset(single_recovery, absolute_error == largest_single)$term
```

In this data set the estimates are
`r join_and(fmt(single_recovery$estimate, 3))`, and the largest absolute
error is `r fmt(largest_single, 3)`, for `` `r largest_term` ``, the sampling
error of one sample of `r nrow(known)`.

## Batches

With `batch = 500`, the same call returns 500 data sets from the same model,
each with its own seed. In `apply_batch()`, the model is fitted to each data
set, and `validate_recovery()` sets every estimate against the truth of its
own data set.

```{r batch}
datasets <- simulate_regression(
  n = 100,
  coefficients = c("(Intercept)" = 1, x1 = 0.5, x2 = -0.3),
  seed = 1,
  batch = 500
)
estimates <- apply_batch(datasets, lm, formula = outcome ~ x1 + x2)
recovery <- validate_recovery(estimates, datasets)
performance <- summarize_simulations(recovery, by = "term", variables = "bias")
performance
```

```{r batch-values, include = FALSE}
largest_bias <- max(abs(performance$mean))
empirical_se <- mean(performance$sd)
```

Over the `r length(datasets)` data sets, the mean bias of each coefficient
estimates the bias of the estimator, here at most `r fmt(largest_bias, 3)` in
absolute value. The standard deviation of the bias is the empirical standard
error, about `r fmt(empirical_se)`, which matches $\sigma/\sqrt{n} = 0.1$ for
`r nrow(known)` observations with $\sigma = 1$.

## Seeds

A call with the same arguments and the same `seed` reproduces the same data.
The random-number state of the session is restored after a seeded call, so
the draws that follow it continue the session's stream.

```{r seeds}
set.seed(2024)
state_before <- .Random.seed
again <- simulate_regression(
  n = 100,
  coefficients = c("(Intercept)" = 1, x1 = 0.5, x2 = -0.3),
  seed = 1
)
identical(again, known)
identical(state_before, .Random.seed)
```

The data of `again` are identical to those of `known`, and the random-number
state of the session is the same before and after the call. With
`seed = NULL`, a call draws from the session's stream and advances it.

## Variables as distributions

A data set can be declared variable by variable, with each variable a call to
a named distribution. The parameters of a distribution may depend on the
variables declared before it. In the declaration below, age is normal,
treatment uptake is binary with log odds linear in age, the number of visits
is Poisson with a log mean linear in age, and the outcome is normal with a
mean that depends on age and treatment.

```{r specification}
specification <- define_variables(
  age     = normal(mean = 50, sd = 10),
  treated = binary(prob = plogis(-2 + 0.04 * age)),
  visits  = poisson(lambda = exp(0.5 + 0.02 * age)),
  outcome = normal(mean = 10 + 0.2 * age + 2 * treated, sd = 2)
)

study <- simulate_study(n = 1000, specification = specification, seed = 42)
head(study)
as.data.frame(study, what = "definitions")
```

The `definitions` table holds the declared parameter expressions. Logistic
and Poisson regressions fitted to the data estimate the coefficients of age.

```{r specification-recovery}
uptake_fit <- glm(treated ~ age, family = binomial, data = study)
coef(uptake_fit)
visits_fit <- glm(visits ~ age, family = poisson, data = study)
coef(visits_fit)
```

The estimated coefficient of age is
`r fmt(coefficient_of(uptake_fit, "age"), 3)` for treatment uptake and
`r fmt(coefficient_of(visits_fit, "age"), 3)` for visits, against the
generating values 0.04 and 0.02. The `r n_distributions` distributions and
their parameters, in the order of a positional call, are listed with
`list_distributions()`.

```{r distributions}
head(list_distributions())
```

## Missing data

Missing values are classified by what the probability of missingness depends
on. Values are missing completely at random when the probability depends on
no variable, missing at random when it depends only on observed values, and
missing not at random when it depends on the values that are missing. With
`mechanism = "MAR"` in `inject_missingness()`, the probability that a value
of `variables` is missing rises with the rank of the fully observed
`predictor`, and the probabilities average `proportion`. The realized share
varies around `proportion` because the removal is random.

```{r missingness}
incomplete <- inject_missingness(
  study,
  mechanism = "MAR",
  proportion = 0.3,
  variables = "outcome",
  predictor = "age",
  seed = 7
)
as.data.frame(incomplete, what = "missingness_summary")
summarize_simulations(study, variables = "outcome")
summarize_simulations(incomplete, variables = "outcome")
```

```{r missingness-values, include = FALSE}
realized_missing <- as.data.frame(incomplete, what = "missingness_summary")$realized_proportion
mean_before <- summarize_simulations(study, variables = "outcome")$mean
mean_after <- summarize_simulations(incomplete, variables = "outcome")$mean
```

The realized share of missing outcomes is `r pct(realized_missing)`, against
the target of 30%. Missingness rises with age, and the outcome rises with
age, so the mean of the observed outcomes, `r fmt(mean_after)`, lies below
the mean before removal, `r fmt(mean_before)`. The mean of the observed
values is biased under missingness at random whenever the variable that
drives missingness is related to the outcome.

The `missingness` table holds the probability and the indicator of every
targeted cell, and `missingness_summary` holds the target and realized
shares. Missingness rules stated as formulas are declared with
`define_missingness()` and applied with `missingness_matrix()`.

## Treatment assignment

In random assignment, observations are allocated to treatment groups at
random, within each stratum when strata are given. With the default balanced
allocation and equal ratios, the observations of each stratum are divided
between the groups as evenly as possible.

```{r randomization}
trial <- assign_treatment(study, groups = 2, strata = "treated", name = "arm", seed = 3)
as.data.frame(trial, what = "allocation")
```

The assigned group is stored in `arm`, and the `allocation` table holds the
number of observations of each stratum in each group. Random allocation fixes
the assignment mechanism, and covariates can still differ between the groups
of a finite sample.

In observational data, exposure depends on covariates. With `link = "logit"`
in `observe_treatment()`, the formula gives the log odds of exposure, which
here rise with age.

```{r propensity}
observed <- observe_treatment(
  study,
  formulas = "-3 + 0.05 * age",
  link = "logit",
  name = "exposed",
  seed = 3
)
age_by_exposure <- summarize_simulations(observed, by = "exposed", variables = "age")
age_by_exposure
```

```{r propensity-values, include = FALSE}
age_exposed <- subset(age_by_exposure, exposed == 1)$mean
age_unexposed <- subset(age_by_exposure, exposed == 0)$mean
```

The mean age is `r fmt(age_exposed, 1)` years among the exposed and
`r fmt(age_unexposed, 1)` years among the unexposed. Age also raises the
outcome in `study`, so a comparison of the exposed with the unexposed is
confounded by age. The `probabilities` table holds the generating
probabilities of exposure, the reference for evaluating the estimation of
propensity scores. The outcome is left as it was, so exposure has no effect
on it.

## Group designs

In a group design, the outcome of each group is normal with its own mean and
standard deviation. Two-group and one-way designs are generated with
`simulate_ttest()` and `simulate_anova()`, from a sample size, a mean and a
standard deviation for each group.

```{r anova}
groups <- simulate_anova(
  n = c(40, 40, 40),
  means = c(10, 11, 12.5),
  sds = 2,
  labels = c("control", "low", "high"),
  seed = 1
)
as.data.frame(groups, what = "effects")
group_fit <- aov(outcome ~ group, data = groups)
summary(group_fit)
```

```{r anova-values, include = FALSE}
design_effects <- as.data.frame(groups, what = "effects")
sample_sums <- aov_sums(group_fit)
sample_eta <- sample_sums[1] / sum(sample_sums)
```

The `effects` table reports $\eta^2 = B/(B+W)$, computed from the design,
where $B = \sum_g n_g(\mu_g-\bar\mu)^2$ is the between-group sum of squares
and $W = \sum_g (n_g-1)\sigma_g^2$ is the within-group sum of squares. For
this design, $B$ is `r fmt(design_effects$between_sum_squares, 1)`, $W$ is
`r fmt(design_effects$within_sum_squares, 0)` and $\eta^2$ is
`r fmt(design_effects$eta_squared, 3)`. Because $W$ uses $n_g - 1$, this
$\eta^2$ is a quantity of the design, and it differs from a ratio of
population variances weighted by group proportions. The fitted analysis of
variance gives sums of squares of `r fmt(sample_sums[1], 1)` between groups
and `r fmt(sample_sums[2], 1)` within groups, and a sample $\eta^2$ of
`r fmt(sample_eta, 3)`.

## Correlation structures

In a first-order autoregressive structure, the correlation between variables
$i$ and $j$ is $\rho^{|i-j|}$, which decreases with the distance between
their positions when $0 < \rho < 1$. With `structure = "ar1"` in
`simulate_correlated()`, the variables follow this structure.

```{r ar1}
correlated <- simulate_correlated(
  n = 2000,
  means = c(t1 = 0, t2 = 0, t3 = 0, t4 = 0),
  rho = 0.6,
  structure = "ar1",
  seed = 1
)
sample_correlation <- cor(subset(correlated, select = -id))
round(sample_correlation, 2)
```

```{r ar1-values, include = FALSE}
lag_one <- off_diagonal(sample_correlation, 1)
lag_two <- off_diagonal(sample_correlation, 2)
lag_three <- off_diagonal(sample_correlation, 3)
```

With $\rho = 0.6$, the generating correlations are 0.6, `r 0.6^2` and
`r 0.6^3` at distances of one, two and three positions. The sample
correlations are `r fmt(min(lag_one))` to `r fmt(max(lag_one))` at distance
one, `r join_and(fmt(lag_two))` at distance two, and `r fmt(lag_three)` at
distance three. The generating correlation matrix is returned with
`as.data.frame(correlated, what = "correlation")`.

## Non-normal margins with a Gaussian copula

A copula describes the dependence between variables separately from their
marginal distributions. In a Gaussian copula, latent variables
$\mathbf{z} \sim N(\mathbf{0}, \mathbf{R})$ are transformed to uniform
variables with the standard normal distribution function $\Phi$, and each
uniform variable is transformed with the quantile function of its margin.

```{r copula}
copula <- simulate_copula(
  n = 2000,
  specification = define_variables(
    score  = normal(mean = 10, sd = 2),
    visits = poisson(lambda = 4),
    rate   = beta(shape1 = 2, shape2 = 5)
  ),
  rho = 0.6,
  seed = 1
)
head(copula, 3)
rank_correlation <- cor(subset(copula, select = -id), method = "spearman")
round(rank_correlation, 2)
```

```{r copula-values, include = FALSE}
continuous_spearman <- (6 / pi) * asin(0.3)
score_rate <- pair_value(rank_correlation, "score", "rate")
score_visits <- pair_value(rank_correlation, "score", "visits")
visits_rate <- pair_value(rank_correlation, "visits", "rate")
```

The argument `rho = 0.6` sets the correlation of the latent normal variables,
and the `latent_correlation` table holds that matrix. The correlations of the
generated variables depend on their margins as well. For two continuous
margins, a latent correlation of 0.6 gives a Spearman correlation of
$(6/\pi)\arcsin(0.3) =$ `r fmt(continuous_spearman)`, and the sample Spearman
correlation of the normal and beta variables is `r fmt(score_rate)`. The
sample correlations of the Poisson variable with the other two are
`r fmt(score_visits)` and `r fmt(visits_rate)`. The formula holds for
continuous margins, and the ties of a discrete margin change the rank
correlation.

## Ordinal variables

In the latent variable model of an ordinal variable, a latent normal variable
is cut into categories at thresholds. In `simulate_ordinal()`, the thresholds
give each category the probability set in `probabilities`, and correlated
latent variables give correlated ordinal variables. This is the model behind
the polychoric correlation.

```{r ordinal}
ordinal <- simulate_ordinal(
  n = 2000,
  probabilities = c(0.2, 0.3, 0.3, 0.2),
  n_variables = 2,
  rho = 0.7,
  seed = 1
)
head(ordinal, 3)
ordinal_correlation <- cor(subset(ordinal, select = -id))
ordinal_correlation
```

The generating latent correlation is 0.7, and the sample Pearson correlation
of the category codes is `r fmt(upper_values(ordinal_correlation), 3)`. The
first is a correlation of the latent variables and the second a correlation
of their ordinal codes. A polychoric estimator estimates the latent
correlation.

## Common factors

In the common factor model, the vector of indicators is

$$\mathbf{x} = \boldsymbol\tau + \boldsymbol\Lambda\boldsymbol\eta + \boldsymbol\varepsilon,$$

where the factors $\boldsymbol\eta$ have unit variance and correlation matrix
$\boldsymbol\Phi$, and the errors are independent with covariance matrix
$\boldsymbol\Psi$. The covariance of the indicators is then
$\boldsymbol\Lambda\boldsymbol\Phi\boldsymbol\Lambda^\top + \boldsymbol\Psi$.
By default, the uniqueness of indicator $j$ is $1-\sum_k\lambda_{jk}^2$,
which gives every indicator unit variance when each loads on one factor. When
an indicator loads on several correlated factors, its variance under this
default differs from one.

```{r factors}
item_loadings <- matrix(
  c(0.8, 0.7, 0.6, 0.0, 0.0, 0.0,
    0.0, 0.0, 0.0, 0.7, 0.6, 0.5),
  ncol = 2,
  dimnames = list(paste0("item_", 1:6), c("verbal", "numeric"))
)
factor_data <- simulate_factors(
  n = 1000,
  loadings = item_loadings,
  factor_correlation = matrix(c(1, 0.3, 0.3, 1), 2),
  seed = 1
)
head(as.data.frame(factor_data, what = "parameters"), 4)
```

The `parameters` table holds the loadings, uniquenesses and intercepts. A
maximum likelihood factor analysis with an oblique rotation estimates the
loadings.

```{r factors-recovery}
factor_fit <- factanal(subset(factor_data, select = -id), factors = 2, rotation = "promax")
estimated_loadings <- round(unclass(loadings(factor_fit)), 2)
estimated_loadings
```

```{r factors-values, include = FALSE}
primary_loadings <- apply(abs(estimated_loadings), 1, max)
cross_loadings <- apply(abs(estimated_loadings), 1, min)
```

The estimated primary loadings of the six items are
`r join_and(fmt(primary_loadings))`, against the generating values 0.8, 0.7,
0.6, 0.7, 0.6 and 0.5, and the largest cross-loading is
`r fmt(max(cross_loadings))`. The order and signs of estimated factors are
arbitrary, so factor solutions are compared after the factors are matched.

## Item response theory

In the two-parameter logistic model of item response theory, ability is
$\theta \sim N(0,1)$ and the probability of a correct response to item $j$ is

$$P(Y_j=1\mid\theta)=\operatorname{logit}^{-1}\{a_j(\theta-b_j)\},$$

where $b_j$ is the difficulty and $a_j$ the discrimination of the item. In
`simulate_irt()`, `model = "rasch"` fixes every discrimination at 1, `"3pl"`
adds a guessing parameter, and `"graded"` gives ordered responses under the
graded response model.

```{r irt}
irt <- simulate_irt(
  n = 2000,
  discrimination = c(0.8, 1.2, 1.5, 1, 2),
  difficulty = c(-1.5, -0.5, 0, 0.5, 1.5),
  seed = 1
)
as.data.frame(irt, what = "parameters")
item_summary <- summarize_simulations(subset(irt, select = -id))
item_summary
```

```{r irt-values, include = FALSE}
item_parameters <- as.data.frame(irt, what = "parameters")
middle_item <- subset(item_parameters, abs(value) < sqrt(.Machine$double.eps))$item
middle_share <- subset(item_summary, variable == middle_item)$mean
```

The proportion correct falls from `r fmt(max(item_summary$mean))` for the
easiest item to `r fmt(min(item_summary$mean))` for the hardest. Because the
ability distribution is symmetric about 0, the marginal probability of a
correct response to the item of difficulty 0 is 0.5 under this model, and its
sample proportion is `r fmt(middle_share, 3)`. The `abilities` table holds
the generated abilities.

## Latent profiles

A latent profile model is a finite mixture of multivariate normal
distributions. Each observation belongs to a profile with the probability of
that profile, and its indicators are drawn with the means, standard
deviations and correlations of its profile. In `simulate_lpa()`, the profile
means can be given as a tidy data frame with `profile`, `variable` and
`mean`.

```{r lpa}
means <- data.frame(
  profile  = rep(c("low", "average", "high"), each = 3),
  variable = rep(c("effort", "planning", "reflection"), times = 3),
  mean     = c(-1, -1, -0.5,  0, 0, 0,  1.2, 1, 1.5)
)
profiles <- simulate_lpa(n = 600, means = means, proportions = c(0.25, 0.5, 0.25), seed = 1)
head(profiles, 3)
effort_by_profile <- summarize_simulations(profiles, by = "profile", variables = "effort")
effort_by_profile
```

```{r lpa-values, include = FALSE}
effort_low <- subset(effort_by_profile, profile == "low")$mean
effort_average <- subset(effort_by_profile, profile == "average")$mean
effort_high <- subset(effort_by_profile, profile == "high")$mean
```

The mean effort is `r fmt(effort_low)`, `r fmt(effort_average)` and
`r fmt(effort_high)` in the low, average and high profiles, against the
generating means −1, 0 and 1.2. The generated `profile` column is the
reference against which a classification by a fitted mixture model is
judged, after the estimated profiles are matched to the generating ones.

## Latent classes

In a latent class model, the indicators are categorical, and each class has
its own probabilities for the categories of each indicator. The indicators
are independent given class membership, the assumption of local
independence.

```{r lca}
item_probabilities <- data.frame(
  class       = rep(c("engaged", "disengaged"), each = 6),
  indicator   = rep(rep(c("attends", "submits", "posts"), each = 2), times = 2),
  category    = rep(c("no", "yes"), times = 6),
  probability = c(0.1, 0.9,  0.2, 0.8,  0.4, 0.6,
                  0.6, 0.4,  0.7, 0.3,  0.9, 0.1)
)
classes <- simulate_lca(n = 500, probabilities = item_probabilities,
                        proportions = c(0.6, 0.4), seed = 1)
head(classes, 3)
xtabs(~ latent_class + attends, data = classes)
```

```{r lca-values, include = FALSE}
engaged_n <- sum(classes$latent_class == "engaged")
engaged_yes <- sum(classes$latent_class == "engaged" & classes$attends == "yes")
disengaged_n <- sum(classes$latent_class == "disengaged")
disengaged_yes <- sum(classes$latent_class == "disengaged" & classes$attends == "yes")
```

The generating probability of attendance is 0.9 in the engaged class and 0.4
in the disengaged class. In this sample, `r engaged_yes` of the
`r engaged_n` engaged participants attend, a share of
`r fmt(engaged_yes / engaged_n)`, and `r disengaged_yes` of the
`r disengaged_n` disengaged participants, a share of
`r fmt(disengaged_yes / disengaged_n)`. The class proportions are
probabilities of membership, so the class sizes vary between samples.

## Two-level latent profiles

In a two-level latent profile model, individuals are nested in clusters, and
each cluster belongs to a latent cluster class that sets the probabilities of
the individual profiles. The profiles have the same means in every cluster
class, and their prevalence differs between cluster classes.

```{r ml-lpa}
nested <- simulate_ml_lpa(
  clusters = 40,
  cluster_size = 25,
  means = means,
  profile_probabilities = data.frame(
    cluster_class = rep(c("mostly_low", "mostly_high"), each = 3),
    profile       = rep(c("low", "average", "high"), times = 2),
    probability   = c(0.6, 0.3, 0.1,  0.1, 0.3, 0.6)
  ),
  seed = 1
)
head(nested, 3)
xtabs(~ cluster_class + profile, data = nested)
```

The `mostly_low` cluster class gives the `low`, `average` and `high` profiles
the probabilities 0.6, 0.3 and 0.1, and the `mostly_high` class gives them
0.1, 0.3 and 0.6. The cross-tabulation counts the individuals of the 40
clusters of 25 by cluster class and profile. The `clusters` table holds the
generated class of each cluster.

## Random intercepts

In the random intercept model, the outcome of observation $i$ in cluster $j$
is

$$y_{ij}=\gamma_0+\gamma_1x_{ij}+u_j+e_{ij},$$

where $u_j\sim N(0,\tau^2)$ and $e_{ij}\sim N(0,\sigma^2)$. Conditional on the
fixed predictors, the intraclass correlation is $\tau^2/(\tau^2+\sigma^2)$,
the correlation between two observations of the same cluster. The
random-intercept variance that gives a target intraclass correlation for a
given within-cluster variance is computed with `calibrate_icc()`.

```{r calibrate-icc}
icc_calibration <- calibrate_icc(icc = 0.2, distribution = "normal", within_variance = 1)
icc_calibration
```

For a target intraclass correlation of 0.2 and a within-cluster variance of
1, the random-intercept variance is
`r icc_calibration$random_effect_variance`, a standard deviation of
`r sqrt(icc_calibration$random_effect_variance)`.

```{r multilevel}
multilevel <- simulate_multilevel(
  clusters = 50,
  cluster_size = 20,
  intercept = 50,
  slopes = c(hours = 2),
  random_intercept_sd = 0.5,
  residual_sd = 1,
  seed = 1
)
head(multilevel, 3)
as.data.frame(multilevel, what = "variance_components")
```

The data have 50 clusters of 20 observations. The `variance_components`
table holds the generating variances, and `random_effects` the realized
effects of the clusters. The intraclass correlation describes the clustering
that remains after the fixed predictor, and the variance explained by `hours`
lies outside its denominator.

## Growth curves

In a linear growth model, the outcome of person $i$ at time $t$ is

$$y_{ti}=(\beta_0+b_{0i})+(\beta_1+b_{1i})t+e_{ti},$$

where $b_{0i}$ and $b_{1i}$ are normal random effects of the person. A
quadratic term adds curvature to the trajectories.

```{r growth}
growth <- simulate_growth(
  n = 200,
  times = 0:4,
  intercept = 10,
  slope = 1.5,
  random_sd = c(2, 0.5),
  seed = 1
)
mean_trajectory <- aggregate(outcome ~ time, data = growth, FUN = mean)
mean_trajectory
```

```{r growth-values, include = FALSE}
mean_first <- subset(mean_trajectory, time == 0)$outcome
mean_last <- subset(mean_trajectory, time == 4)$outcome
```

The generating mean trajectory is $10 + 1.5t$, which is 10 at time 0 and
`r 10 + 1.5 * 4` at time 4. The sample means rise from `r fmt(mean_first)` at
time 0 to `r fmt(mean_last)` at time 4.

## Vector autoregression

In a first-order vector autoregression with a zero intercept, the vector of
variables of person $i$ at occasion $t$ is

$$\mathbf{y}_{it}=\boldsymbol\mu_i+\mathbf{A}(\mathbf{y}_{i,t-1}-\boldsymbol\mu_i)+\boldsymbol\varepsilon_{it},$$

where the person means $\boldsymbol\mu_i$ vary with the given between-person
covariance. In `simulate_longitudinal()`, `transition` holds the lagged
coefficients $\mathbf{A}$. In a tidy data frame, `from` is the variable at
the previous occasion and `to` the variable at the current occasion.

```{r var}
lagged_effects <- data.frame(
  from        = c("mood", "mood", "stress", "stress"),
  to          = c("mood", "stress", "mood", "stress"),
  coefficient = c(0.5, 0, -0.2, 0.4)
)
experience <- simulate_longitudinal(
  n = 50,
  occasions = 40,
  transition = lagged_effects,
  seed = 1
)
head(experience, 3)
as.data.frame(experience, what = "transition")
```

Mood has an autoregressive coefficient of 0.5 and stress one of 0.4. Stress
at the previous occasion lowers mood with a coefficient of −0.2, and the
cross-lagged coefficient of mood on stress is 0. In the `transition` table,
rows are the variables predicted and columns the lagged predictors, as in
$\mathbf{A}$, and the `wide` table holds the repeated measurements of each
person in one row. With `beeps_per_day`, the occasions are grouped into days
of experience-sampling prompts, and the dependence on the previous occasion
is reset at the first prompt of each day.

## Proportional hazards

In the proportional hazards model, the hazard of an event at time $t$ for
covariates $\mathbf{x}$ is

$$h(t\mid\mathbf{x})=h_0(t)\exp(\mathbf{x}^\top\boldsymbol\beta).$$

For a Weibull baseline with cumulative hazard $H_0(t)=\lambda t^k$, the event
time is obtained by inverse transformation of a uniform variable,

$$T=\left\{\frac{-\log U}{\lambda\exp(\mathbf{x}^\top\boldsymbol\beta)}\right\}^{1/k},\qquad U\sim U(0,1).$$

The censoring times are exponential, with a rate calibrated so that the mean
probability of censoring over the generated event times equals `censoring`.

```{r survival}
survival_data <- simulate_proportional_survival(
  n = 2000,
  coefficients = c(treatment = -0.5, age = 0.3),
  baseline = "weibull",
  rate = 0.05,
  shape = 1.5,
  censoring = 0.3,
  seed = 1
)
head(survival_data, 3)
as.data.frame(survival_data, what = "parameters")
as.data.frame(survival_data, what = "diagnostics")
```

```{r survival-values, include = FALSE}
realized_censoring <- as.data.frame(survival_data, what = "diagnostics")$realized_censoring
```

The coefficients are log hazard ratios. Covariates are standard normal by
default, the covariate named `treatment` included, and
`covariate_distribution = "binary"` gives binary covariates. The realized
share of censored times is `r pct(realized_censoring)`, against the target
of 30%. The `parameters` table holds the coefficients and the parameters of
the hazard, and `diagnostics` the target and realized censoring. A Cox model
estimates the coefficients from the observed times and event indicators.

```{r survival-recovery, eval = has_survival}
cox_fit <- survival::coxph(survival::Surv(time, status) ~ treatment + age, data = survival_data)
coef(cox_fit)
```

`r if (has_survival) sprintf("The estimates are %s and %s, against the generating values −0.5 and 0.3.", fmt(coefficient_of(cox_fit, "treatment"), 3), fmt(coefficient_of(cox_fit, "age"), 3))`

## Competing risks

Cause-specific models of event times are declared with `define_survivals()`,
and their latent event times are generated with `simulate_survival()`. In
`combine_competing_risks()`, the earliest time and its event type are kept,
and an event named in `censor` is treated as censoring.

```{r competing-risks}
latent_times <- simulate_survival(
  n = 500,
  specification = define_survivals(
    relapse = hazard(log_rate = -3 + 0.5 * treatment),
    dropout = hazard(log_rate = -3.5)
  ),
  covariates = data.frame(id = 1:500, treatment = rep(0:1, each = 250)),
  seed = 1
)
observed_times <- combine_competing_risks(latent_times, events = c("relapse", "dropout"),
                                          censor = "dropout")
head(observed_times, 3)
as.data.frame(observed_times, what = "events")
xtabs(~ treatment + event_type, data = observed_times)
```

```{r competing-risks-values, include = FALSE}
relapse_treated <- sum(observed_times$treatment == 1 & observed_times$event_type == "relapse")
relapse_untreated <- sum(observed_times$treatment == 0 & observed_times$event_type == "relapse")
```

The `events` table codes dropout as 0 and relapse as 1. Dropout is treated as
censoring, so the outcome is relapse censored by dropout. The treatment
coefficient of 0.5 multiplies the hazard of relapse by
$e^{0.5} =$ `r fmt(exp(0.5))`. Relapse comes first for `r relapse_treated`
of the 250 treated and `r relapse_untreated` of the 250 untreated
participants.

## Sequences of states

In a first-order Markov chain, the next state depends only on the current
state, through the transition probabilities
$p_{ij}=P(X_t=j\mid X_{t-1}=i)$. Called with no arguments,
`simulate_sequences()` draws a transition matrix over five states from a
Dirichlet distribution and generates 100 sequences of 20 states from it, as a
tidy data frame, one row per sequence and position. The drawn matrix is
printed under the data as the truth.

```{r sequences}
sequences <- simulate_sequences(seed = 1)
sequences
```

The estimated transition probability from one state to another is the number
of transitions between them divided by the number of transitions from the
first state, the maximum likelihood estimate for a first-order chain. In
`validate_recovery()`, the estimates are matched to the generating
probabilities by `from` and `to`.

```{r sequences-recovery}
sequence_recovery <- validate_recovery(summarize_transitions(sequences), sequences,
                                       term = "to", estimate = "probability")
head(sequence_recovery)
```

The `r nrow(sequence_recovery)` estimated probabilities differ from the
generating ones by at most `r fmt(max(sequence_recovery$absolute_error), 3)`.
A given transition matrix is passed in `transition`, as a matrix or as a tidy
data frame with `from`, `to` and `probability`.

## Hidden Markov models

In a hidden Markov model, a latent state follows a Markov chain, and at each
occasion the state emits an observed category with the probabilities of that
state. In `simulate_hmm()`, `transition` holds the transition probabilities
of the states and `emission` the probabilities of the categories in each
state, and the result holds both the generated `state` and the emitted
`observation`.

```{r hmm}
hidden <- simulate_hmm(
  n = 200,
  chain_length = 30,
  transition = data.frame(
    from = c("engaged", "engaged", "stuck", "stuck"),
    to   = c("engaged", "stuck", "engaged", "stuck"),
    probability = c(0.9, 0.1, 0.3, 0.7)
  ),
  emission = data.frame(
    state       = rep(c("engaged", "stuck"), each = 3),
    observation = rep(c("submit", "view", "idle"), times = 2),
    probability = c(0.6, 0.3, 0.1,  0.1, 0.3, 0.6)
  ),
  initial = c(engaged = 0.5, stuck = 0.5),
  seed = 1
)
head(hidden, 3)
xtabs(~ state + observation, data = hidden)
```

The engaged state emits `submit` with probability 0.6, and the stuck state
emits `idle` with probability 0.6. The latent sequence is the reference
against which the states reconstructed by a fitted model are judged, after
the estimated states are matched to the generating ones.

## Groups, mixtures and behavioral data

In `simulate_group_sequences()`, sequences are generated for known groups,
each with its own transition matrix. In `simulate_sequence_clusters()`, they
are generated from a mixture of transition matrices, and the generated
membership is stored in `sequence_cluster`.

Behavioral data are generated with `simulate_event_log()` as a tidy data
frame, one row per event, with the group, the actor, the course, the
achievement level, the position of the event in the actor's sequence, the
state and a timestamp. The states are drawn from the catalogue of learning
states listed with `learning_states()`.

```{r behavior}
behavior <- simulate_event_log(groups = 2, actors = 3, n_states = 4,
                               sequence_length = c(5, 8), seed = 1)
head(behavior, 4)
```

The example has 2 groups of 3 actors, 4 states, and sequences of 5 to 8
events for each actor.

## Transition network analysis

A transition network represents the transitions between states as a weighted
directed network, in which the weight of an edge is the probability of moving
from its source state to its target state. Called with no arguments,
`simulate_tna()` draws a transition matrix over five states, labeled with
metacognitive and cognitive states from the package's catalogue of learning
states, generates 100 sequences of 20 states from it, and fits one network to
them with the `tna` package. The result is the native `tna` model, so the functions of the `tna`
package apply to it directly. Its print shows the fitted network, followed by
the transition matrix and initial probabilities that generated the sequences.

```{r tna, eval = has_tna}
network <- simulate_tna(seed = 1)
network
```

A given matrix can come from a network estimated on real data. Here it is the
network the `tna` package estimates from its `group_regulation` data, the
regulation states of students in collaborative learning. Its transition
matrix and initial probabilities generate 2,000 sequences, and the network is
fitted again from them.

```{r tna-given, eval = has_tna}
regulation <- tna::tna(tna::group_regulation)
known_network <- simulate_tna(n = 2000, transition = regulation$weights,
                              initial = regulation$inits, seed = 1)
known_network
```

The fitted network and its truth sit in the same result, so the result is
compared with itself, edge by edge.

```{r tna-recovery, eval = has_tna}
network_recovery <- validate_recovery(known_network, known_network,
                                      term = "to", estimate = "weight")
summarize_simulations(network_recovery, variables = "absolute_error")
```

`r if (has_tna) sprintf("Over the %d edges, the mean absolute error is %s and the largest is %s.", nrow(network_recovery), fmt(mean(network_recovery$absolute_error), 3), fmt(max(network_recovery$absolute_error), 3))`
In `simulate_group_tna()`, one network is simulated and fitted for each of
several groups.

## Random graphs and blocks

In the Bernoulli model, each pair of nodes is connected independently with
probability $p$. Pairs are ordered in a directed network and unordered in an
undirected one, and self-loops are excluded by default. When the probability
is a matrix indexed by node type, the model is a stochastic block model, in
which the probability of a tie depends on the types of both nodes.

```{r block}
node_types <- rep(c("teacher", "student"), times = c(3, 12))
tie_probability <- matrix(
  c(0.1, 0.5,
    0.5, 0.15),
  nrow = 2,
  dimnames = list(c("teacher", "student"), c("teacher", "student"))
)
classroom <- simulate_network(
  nodes = 15,
  probability = tie_probability,
  node_type = node_types,
  directed = FALSE,
  seed = 1
)
head(classroom, 3)
as.data.frame(classroom, what = "settings")
```

```{r block-values, include = FALSE}
classroom_edges <- as.data.frame(classroom, what = "settings")$edges
```

The network has 3 teachers and 12 students. The probability of a tie is 0.1
between teachers, 0.15 between students and 0.5 between a teacher and a
student. The result is a tidy data frame, one row per edge. The `nodes` table
holds the type of each node, `adjacency` the adjacency matrix as a tidy data
frame, and `settings` the realized number of edges, here
`r classroom_edges`.

## Preferential attachment and small worlds

In preferential attachment, a network grows by adding nodes, and each new
node connects to existing nodes with probabilities that rise with their
degree. A small-world network starts from ties between neighboring nodes and
rewires a share of them at random. Both models are generated through
`igraph`. Centrality measures of the nodes are returned as a tidy data frame
with `network_centrality()`.

```{r attachment, eval = has_igraph}
scale_free <- simulate_network(nodes = 200, model = "barabasi_albert",
                               attachment = 2, directed = FALSE, seed = 1)
degree <- network_centrality(scale_free, measures = "degree")
degree_summary <- summarize_simulations(degree, variables = "value")
degree_summary
```

`r if (has_igraph) sprintf("With %s, each new node adds two edges to existing nodes. The degree ranges from %s to %s, with a mean of %s.", "\x60attachment = 2\x60", degree_summary$minimum, degree_summary$maximum, fmt(degree_summary$mean))`

## Temporal, bipartite and multiplex networks

In the temporal network model, time is divided into periods, and each pair of
nodes is active or inactive in each period. The probability that a pair is
active in the first period is `initial_probability`, the probability that an
inactive pair becomes active is `formation_probability`, and the probability
that an active pair becomes inactive is `dissolution_probability`.

```{r temporal}
temporal <- simulate_temporal_network(
  nodes = 10,
  periods = 12,
  initial_probability = 0.1,
  formation_probability = 0.05,
  dissolution_probability = 0.3,
  seed = 1
)
head(temporal, 4)
components(temporal)
```

The result is a tidy data frame of spells, one row per period of continuous
activity of a pair, with `from`, `to`, `onset` and `terminus`. A spell is
active from `onset` to `terminus - 1`. A spell still active in the last
period has a `terminus` one past the last period and `censored = TRUE`. The
`events` table holds the formations and dissolutions, and `snapshots` the
active pairs in each period.

A bipartite network has ties only between two disjoint sets of nodes, such as
students and the events they attend, and is generated with
`simulate_bipartite_network()`. A multiplex network has several layers of
ties over the same nodes, each layer with its own probability, and is
generated with `simulate_multiplex_network()`.

```{r multiplex}
layers <- simulate_multiplex_network(
  nodes = 10,
  layers = c("advice", "friendship"),
  probability = c(0.15, 0.3),
  seed = 1
)
as.data.frame(layers, what = "layers")
xtabs(~ layer, data = layers)
```

```{r multiplex-values, include = FALSE}
ordered_pairs <- 10 * 9
advice_edges <- sum(layers$layer == "advice")
friendship_edges <- sum(layers$layer == "friendship")
```

The network has 10 nodes in two directed layers without self-loops. Each
layer has `r ordered_pairs` ordered pairs, so the expected number of edges is
`r ordered_pairs * 0.15` for advice and `r ordered_pairs * 0.3` for
friendship. The realized numbers are `r advice_edges` and
`r friendship_edges`. The `layers` table holds the settings of each layer.

## Simulation functions

Every function below takes `seed`, and every single-dataset function also
takes `batch`. The catalogue, with the shape of each result, is listed with
`list_simulators()`.

| Function | Data generated |
|---|---|
| `simulate_study()` | Variables declared one by one from named distributions and formulas |
| `simulate_correlation()` | Normal variables with a given correlation structure, with defaults for every argument but `n` |
| `simulate_correlated()` | Multivariate normal variables with given means, standard deviations and correlations |
| `simulate_copula()` | Variables with any marginal distributions, correlated through a Gaussian copula |
| `simulate_ordinal()` | Correlated ordinal variables cut from latent normal variables |
| `simulate_ttest()` | Two groups with normal outcomes |
| `simulate_anova()` | Several groups with normal outcomes |
| `simulate_regression()` | A linear regression with correlated normal predictors |
| `simulate_prediction()` | A prediction design with continuous and categorical predictors |
| `simulate_clusters()` | Gaussian clusters with known membership |
| `simulate_lpa()` | A latent profile model |
| `simulate_ml_lpa()` | A two-level latent profile model, individuals in clusters |
| `simulate_lca()` | A latent class model with categorical indicators |
| `simulate_factors()` | A common factor model |
| `simulate_irt()` | Item responses under the Rasch, two- and three-parameter logistic and graded response models |
| `simulate_multilevel()` | A two-level model with random intercepts and slopes |
| `simulate_growth()` | Growth trajectories with random intercepts, slopes and curvature |
| `simulate_longitudinal()` | A multivariate vector autoregression, with optional days of experience-sampling prompts |
| `simulate_survival()` | Event times from cause-specific hazards, for competing risks |
| `simulate_proportional_survival()` | Proportional hazards event times with calibrated censoring |
| `simulate_markov()` | Markov chains from a given transition matrix |
| `generate_transition_system()` | A transition matrix and initial probabilities drawn from a Dirichlet distribution |
| `simulate_sequences()` | Sequences of states, with a drawn or given transition matrix |
| `simulate_sequence_clusters()` | Sequences from a mixture of transition matrices, with cluster membership |
| `simulate_hmm()` | A hidden Markov model, with hidden states and emitted observations |
| `simulate_group_sequences()` | Sequences of actors in groups, each group with its own transition matrix |
| `simulate_event_log()` | Behavioral events of actors in groups and courses, with timestamps and achievement levels |
| `simulate_until_event()` | Records of each unit through its nth event |
| `simulate_tna()` | One transition network, fitted to sequences from a drawn or given matrix |
| `simulate_group_tna()` | One transition network per group |
| `simulate_tna_network()` | A transition network with nodes in groups and given within- and between-group tie probabilities |
| `simulate_network()` | Random, block, preferential-attachment, small-world, regular, geometric and forest-fire networks |
| `simulate_edge_list()` | A tidy edge list with optional node types and edge classes |
| `simulate_network_matrix()` | Adjacency, transition, frequency or co-occurrence matrices |
| `simulate_temporal_network()` | Edges that form and dissolve between periods |
| `simulate_bipartite_network()` | Ties between two disjoint sets of nodes |
| `simulate_multiplex_network()` | Several layers of ties over the same nodes |
| `simulate_synthetic()` | Rows resampled from an observed data set |
| `simulate_density()` | Values drawn from the kernel density estimate of an observed variable |
| `simulate_spline()` | A variable generated from a spline of a predictor |
| `simulate_sequence_batches()` | Repeated sequence data sets, stacked |
| `simulate_tna_batches()` | Repeated fitted transition networks, stacked |
| `simulate_network_batches()` | Repeated networks, stacked |
| `simulate_scenarios()` | Any simulator run over a grid of design values |
| `simulate_data()` | Any simulator, named as a string |