Simulating a Wide Variety of Data Shapes, Sizes and Distributions

Mohammed Saqr

University of Eastern Finland

Sonsoles López-Pernas

University of Eastern Finland

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.

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.

drawn <- simulate_regression(seed = 1)
drawn
#> <simulab_sim:regression> 100 rows x 4 columns
#>    id         x1         x2     outcome
#> 1   1 -0.8356286 -0.9109216  0.08428839
#> 2   2  1.5952808  0.1580288 -0.59203209
#> 3   3  0.3295078 -0.6545846 -3.40256878
#> 4   4 -0.8204684  1.7672873  3.19376689
#> 5   5  0.4874291  0.7167075  0.69788065
#> 6   6  0.7383247  0.9101742  0.76841891
#> 7   7  0.5757814  0.3841854 -0.24200033
#> 8   8 -0.3053884  1.6821761  1.21368455
#> 9   9  1.5117812 -0.6357365 -1.08891255
#> 10 10  0.3898432 -0.4616447 -0.52937755
#> ... 90 more rows
#> 
#> Truth (coefficients):
#>          term coefficient
#> 1 (Intercept)       -0.63
#> 2          x1        0.15
#> 3          x2        0.82
#> 
#> Other tables: effects, predictor_correlation. Read one with as.data.frame(x, what = "effects").

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.

known <- simulate_regression(
  n = 100,
  coefficients = c("(Intercept)" = 1, x1 = 0.5, x2 = -0.3),
  seed = 1
)
known
#> <simulab_sim:regression> 100 rows x 4 columns
#>    id         x1          x2    outcome
#> 1   1 -0.6264538 -0.62036668  1.2822849
#> 2   2  0.1836433  0.04211587  2.7680602
#> 3   3 -0.8356286 -0.91092165  2.4420506
#> 4   4  1.5952808  0.15802877  1.4193240
#> 5   5  0.3295078 -0.65458464 -0.9241063
#> 6   6 -0.8204684  1.76728727  2.5572412
#> 7   7  0.4874291  0.71670748  1.6957685
#> 8   8  0.7383247  0.91017423  1.6374374
#> 9   9  0.5757814  0.38418536  1.1592355
#> 10 10 -0.3053884  1.68217608  0.8527614
#> ... 90 more rows
#> 
#> Truth (coefficients):
#>          term coefficient
#> 1 (Intercept)         1.0
#> 2          x1         0.5
#> 3          x2        -0.3
#> 
#> Other tables: effects, predictor_correlation. Read one with as.data.frame(x, what = "effects").

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.

fit <- lm(outcome ~ x1 + x2, data = known)
validate_recovery(fit, known)
#>          term   estimate truth        bias absolute_error relative_error
#> 1 (Intercept)  1.0253534   1.0  0.02535343     0.02535343     0.02535343
#> 2          x1  0.5211102   0.5  0.02111017     0.02111017     0.04222033
#> 3          x2 -0.3534668  -0.3 -0.05346682     0.05346682     0.17822272
#>   recovered
#> 1      TRUE
#> 2      TRUE
#> 3      TRUE

In this data set the estimates are 1.025, 0.521 and −0.353, and the largest absolute error is 0.053, for x2, the sampling error of one sample of 100.

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.

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
#>          term variable observations         mean        sd    minimum   maximum
#> 1 (Intercept)     bias          500  0.004524666 0.1003488 -0.3130872 0.3106459
#> 2          x1     bias          500  0.001262413 0.1018302 -0.2850844 0.3985527
#> 3          x2     bias          500 -0.001266961 0.1003781 -0.2668912 0.3682355

Over the 500 data sets, the mean bias of each coefficient estimates the bias of the estimator, here at most 0.005 in absolute value. The standard deviation of the bias is the empirical standard error, about 0.10, which matches $\sigma/\sqrt{n} = 0.1$ for 100 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.

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)
#> [1] TRUE
identical(state_before, .Random.seed)
#> [1] TRUE

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.

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)
#> <simulab_sim:study> 6 rows x 5 columns
#>   id      age treated visits  outcome
#> 1  1 63.70958       0      4 23.24307
#> 2  2 44.35302       0      7 18.31476
#> 3  3 53.63128       0      4 17.27679
#> 4  4 56.32863       0      5 17.25232
#> 5  5 54.04268       0      3 18.22492
#> 6  6 48.93875       1      5 22.51943
#> 
#> Truth (definitions):
#>   variable distribution parameter                        value
#> 1      age       normal      mean                           50
#> 2      age       normal        sd                           10
#> 3  treated       binary      prob      plogis(-2 + 0.04 * age)
#> 4   visits      poisson    lambda        exp(0.5 + 0.02 * age)
#> 5  outcome       normal      mean 10 + 0.2 * age + 2 * treated
#> 6  outcome       normal        sd                            2
as.data.frame(study, what = "definitions")
#>   variable distribution parameter                        value
#> 1      age       normal      mean                           50
#> 2      age       normal        sd                           10
#> 3  treated       binary      prob      plogis(-2 + 0.04 * age)
#> 4   visits      poisson    lambda        exp(0.5 + 0.02 * age)
#> 5  outcome       normal      mean 10 + 0.2 * age + 2 * treated
#> 6  outcome       normal        sd                            2

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

uptake_fit <- glm(treated ~ age, family = binomial, data = study)
coef(uptake_fit)
#> (Intercept)         age 
#> -1.78449485  0.03678837
visits_fit <- glm(visits ~ age, family = poisson, data = study)
coef(visits_fit)
#> (Intercept)         age 
#>  0.48526262  0.02077051

The estimated coefficient of age is 0.037 for treatment uptake and 0.021 for visits, against the generating values 0.04 and 0.02. The 79 distributions and their parameters, in the order of a positional call, are listed with list_distributions().

head(list_distributions())
#>    distribution            parameters n_parameters copula
#> 1        anglit       location, scale            2   TRUE
#> 2       arcsine              min, max            2   TRUE
#> 3        benini shape1, shape2, scale            3   TRUE
#> 4          beta        shape1, shape2            2   TRUE
#> 5 beta_binomial  size, shape1, shape2            3  FALSE
#> 6    beta_prime        shape1, shape2            2   TRUE

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.

incomplete <- inject_missingness(
  study,
  mechanism = "MAR",
  proportion = 0.3,
  variables = "outcome",
  predictor = "age",
  seed = 7
)
as.data.frame(incomplete, what = "missingness_summary")
#>   variable realized_proportion target_proportion
#> 1  outcome               0.308               0.3
summarize_simulations(study, variables = "outcome")
#>   variable observations     mean       sd  minimum  maximum
#> 1  outcome         1000 20.96394 3.111055 10.46787 30.26189
summarize_simulations(incomplete, variables = "outcome")
#>   variable observations    mean       sd  minimum  maximum
#> 1  outcome          692 20.3551 3.077688 10.46787 29.75574

The realized share of missing outcomes is 30.8%, against the target of 30%. Missingness rises with age, and the outcome rises with age, so the mean of the observed outcomes, 20.36, lies below the mean before removal, 20.96. 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.

trial <- assign_treatment(study, groups = 2, strata = "treated", name = "arm", seed = 3)
as.data.frame(trial, what = "allocation")
#>   stratum treatment observations
#> 1       0         0          245
#> 2       1         0          256
#> 3       0         1          244
#> 4       1         1          255

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.

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
#>   exposed variable observations     mean        sd  minimum  maximum
#> 1       0      age          613 48.03767  9.494049 19.82067 76.23495
#> 2       1      age          387 52.44100 10.260949 16.28261 84.95304

The mean age is 52.4 years among the exposed and 48.0 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.

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")
#>   grand_mean between_sum_squares within_sum_squares eta_squared
#> 1   11.16667            126.6667                468   0.2130045
group_fit <- aov(outcome ~ group, data = groups)
summary(group_fit)
#>              Df Sum Sq Mean Sq F value   Pr(>F)    
#> group         2  131.3   65.64   20.77 1.91e-08 ***
#> Residuals   117  369.8    3.16                     
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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 126.7, $W$ is 468 and $\eta^2$ is 0.213. 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 131.3 between groups and 369.8 within groups, and a sample $\eta^2$ of 0.262.

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.

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)
#>      t1   t2   t3   t4
#> t1 1.00 0.61 0.34 0.18
#> t2 0.61 1.00 0.55 0.33
#> t3 0.34 0.55 1.00 0.59
#> t4 0.18 0.33 0.59 1.00

With $\rho = 0.6$, the generating correlations are 0.6, 0.36 and 0.216 at distances of one, two and three positions. The sample correlations are 0.55 to 0.61 at distance one, 0.34 and 0.33 at distance two, and 0.18 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.

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)
#> <simulab_sim:copula> 3 rows x 4 columns
#>   id    score visits       rate
#> 1  1 7.706109      2 0.07629022
#> 2  2 9.679819      1 0.30106191
#> 3  3 9.711422      7 0.40336193
#> 
#> 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").
rank_correlation <- cor(subset(copula, select = -id), method = "spearman")
round(rank_correlation, 2)
#>        score visits rate
#> score   1.00   0.58 0.57
#> visits  0.58   1.00 0.52
#> rate    0.57   0.52 1.00

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) =$ 0.58, and the sample Spearman correlation of the normal and beta variables is 0.57. The sample correlations of the Poisson variable with the other two are 0.58 and 0.52. 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.

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)
#> <simulab_sim:ordinal> 3 rows x 3 columns
#>   id V1 V2
#> 1  1  1  1
#> 2  2  2  1
#> 3  3  2  4
#> 
#> Truth (probabilities):
#>   variable category probability
#> 1       V1        1         0.2
#> 2       V1        2         0.3
#> 3       V1        3         0.3
#> 4       V1        4         0.2
#> 5       V2        1         0.2
#> 6       V2        2         0.3
#> 7       V2        3         0.3
#> 8       V2        4         0.2
#> 
#> Other tables: correlation. Read one with as.data.frame(x, what = "correlation").
ordinal_correlation <- cor(subset(ordinal, select = -id))
ordinal_correlation
#>           V1        V2
#> V1 1.0000000 0.6583364
#> V2 0.6583364 1.0000000

The generating latent correlation is 0.7, and the sample Pearson correlation of the category codes is 0.658. 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.

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)
#>   variable factor loading uniqueness intercept
#> 1   item_1 verbal     0.8       0.36         0
#> 2   item_2 verbal     0.7       0.51         0
#> 3   item_3 verbal     0.6       0.64         0
#> 4   item_4 verbal     0.0       0.51         0

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

factor_fit <- factanal(subset(factor_data, select = -id), factors = 2, rotation = "promax")
estimated_loadings <- round(unclass(loadings(factor_fit)), 2)
estimated_loadings
#>        Factor1 Factor2
#> item_1    0.78    0.01
#> item_2    0.75   -0.03
#> item_3    0.58    0.00
#> item_4   -0.02    0.72
#> item_5    0.00    0.61
#> item_6    0.03    0.47

The estimated primary loadings of the six items are 0.78, 0.75, 0.58, 0.72, 0.61 and 0.47, against the generating values 0.8, 0.7, 0.6, 0.7, 0.6 and 0.5, and the largest cross-loading is 0.03. 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.

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  parameter category value discrimination guessing
#> 1 item_1 difficulty       NA  -1.5            0.8        0
#> 2 item_2 difficulty       NA  -0.5            1.2        0
#> 3 item_3 difficulty       NA   0.0            1.5        0
#> 4 item_4 difficulty       NA   0.5            1.0        0
#> 5 item_5 difficulty       NA   1.5            2.0        0
item_summary <- summarize_simulations(subset(irt, select = -id))
item_summary
#>   variable observations   mean        sd minimum maximum
#> 1   item_1         2000 0.7185 0.4498432       0       1
#> 2   item_2         2000 0.6160 0.4864795       0       1
#> 3   item_3         2000 0.5120 0.4999810       0       1
#> 4   item_4         2000 0.3905 0.4879844       0       1
#> 5   item_5         2000 0.1320 0.3385754       0       1

The proportion correct falls from 0.72 for the easiest item to 0.13 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 0.512. 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.

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)
#> <simulab_sim:lpa> 3 rows x 5 columns
#>   id profile    effort  planning reflection
#> 1  1 average 1.1948510 -1.010465  -1.113169
#> 2  2 average 0.4955447  2.401222  -1.405505
#> 3  3    high 1.1180083  1.207912   1.170042
#> 
#> Truth (parameters):
#>   profile   variable mean sd proportion
#> 1     low     effort -1.0  1       0.25
#> 2 average     effort  0.0  1       0.50
#> 3    high     effort  1.2  1       0.25
#> 4     low   planning -1.0  1       0.25
#> 5 average   planning  0.0  1       0.50
#> 6    high   planning  1.0  1       0.25
#> 7     low reflection -0.5  1       0.25
#> 8 average reflection  0.0  1       0.50
#> 9    high reflection  1.5  1       0.25
effort_by_profile <- summarize_simulations(profiles, by = "profile", variables = "effort")
effort_by_profile
#>   profile variable observations         mean       sd   minimum   maximum
#> 1 average   effort          314  0.001224158 1.056226 -2.996949 3.0557424
#> 2    high   effort          140  1.155264445 1.131875 -1.349107 4.3539714
#> 3     low   effort          146 -0.979527540 1.015811 -4.008049 0.9713374

The mean effort is −0.98, 0.00 and 1.16 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.

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)
#> <simulab_sim:lca> 3 rows x 5 columns
#>   id latent_class attends submits posts
#> 1  1      engaged     yes     yes   yes
#> 2  2      engaged     yes     yes    no
#> 3  3      engaged     yes     yes   yes
#> 
#> Truth (parameters):
#>    latent_class indicator category probability proportion
#> 1       engaged   attends       no         0.1        0.6
#> 2    disengaged   attends       no         0.6        0.4
#> 3       engaged   submits       no         0.2        0.6
#> 4    disengaged   submits       no         0.7        0.4
#> 5       engaged     posts       no         0.4        0.6
#> 6    disengaged     posts       no         0.9        0.4
#> 7       engaged   attends      yes         0.9        0.6
#> 8    disengaged   attends      yes         0.4        0.4
#> 9       engaged   submits      yes         0.8        0.6
#> 10   disengaged   submits      yes         0.3        0.4
#> ... 2 more rows
xtabs(~ latent_class + attends, data = classes)
#>             attends
#> latent_class  no yes
#>   disengaged 102  86
#>   engaged     36 276

The generating probability of attendance is 0.9 in the engaged class and 0.4 in the disengaged class. In this sample, 276 of the 312 engaged participants attend, a share of 0.88, and 86 of the 188 disengaged participants, a share of 0.46. 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.

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)
#> <simulab_sim:ml_lpa> 3 rows x 7 columns
#>   id cluster cluster_class profile   effort planning reflection
#> 1  1       1   mostly_high    high 1.760049 1.234904  3.3405227
#> 2  2       1   mostly_high    high 1.342451 1.701258 -0.4621991
#> 3  3       1   mostly_high    high 3.029672 2.326598  0.4338980
#> 
#> Truth (parameters):
#>   profile   variable mean sd proportion
#> 1     low     effort -1.0  1       0.35
#> 2 average     effort  0.0  1       0.30
#> 3    high     effort  1.2  1       0.35
#> 4     low   planning -1.0  1       0.35
#> 5 average   planning  0.0  1       0.30
#> 6    high   planning  1.0  1       0.35
#> 7     low reflection -0.5  1       0.35
#> 8 average reflection  0.0  1       0.30
#> 9    high reflection  1.5  1       0.35
#> 
#> Other tables: profile_probabilities, clusters. Read one with as.data.frame(x, what = "profile_probabilities").
xtabs(~ cluster_class + profile, data = nested)
#>              profile
#> cluster_class average high low
#>   mostly_high     154  314  57
#>   mostly_low      129   58 288

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().

icc_calibration <- calibrate_icc(icc = 0.2, distribution = "normal", within_variance = 1)
icc_calibration
#>   icc distribution random_effect_variance
#> 1 0.2       normal                   0.25

For a target intraclass correlation of 0.2 and a within-cluster variance of 1, the random-intercept variance is 0.25, a standard deviation of 0.5.

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)
#> <simulab_sim:multilevel> 3 rows x 4 columns
#>   id cluster       hours  outcome
#> 1  1       1 -0.62036668 48.68794
#> 2  2       1  0.04211587 48.63825
#> 3  3       1 -0.91092165 49.35484
#> 
#> Truth (fixed_effects):
#>          term coefficient
#> 1 (Intercept)          50
#> 2       hours           2
#> 
#> Other tables: variance_components, random_effects. Read one with as.data.frame(x, what = "variance_components").
as.data.frame(multilevel, what = "variance_components")
#>          component variance
#> 1 random_intercept     0.25
#> 2     random_slope     0.00
#> 3         residual     1.00

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.

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
#>   time  outcome
#> 1    0 10.02925
#> 2    1 11.49174
#> 3    2 13.11879
#> 4    3 14.55091
#> 5    4 16.18706

The generating mean trajectory is $10 + 1.5t$, which is 10 at time 0 and 16 at time 4. The sample means rise from 10.03 at time 0 to 16.19 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.

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)
#> <simulab_sim:longitudinal_var> 3 rows x 4 columns
#>   id occasion       mood     stress
#> 1  1        1 -1.3004499  0.7078784
#> 2  1        2 -0.4608929  1.8697398
#> 3  1        3 -3.1020560 -1.5373396
#> 
#> Truth (transition):
#>      row column coefficient
#> 1   mood   mood         0.5
#> 2 stress   mood         0.0
#> 3   mood stress        -0.2
#> 4 stress stress         0.4
#> 
#> Other tables: wide, innovation_covariance, between_covariance, person_means. Read one with as.data.frame(x, what = "wide").
as.data.frame(experience, what = "transition")
#>      row column coefficient
#> 1   mood   mood         0.5
#> 2 stress   mood         0.0
#> 3   mood stress        -0.2
#> 4 stress stress         0.4

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.

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)
#> <simulab_sim:proportional_survival> 3 rows x 5 columns
#>   id      time status  treatment        age
#> 1  1 11.534176      1 -0.6264538 -0.8861496
#> 2  2  6.293952      0  0.1836433 -1.9222549
#> 3  3  1.606610      1 -0.8356286  1.6197007
#> 
#> Truth (parameters):
#>             term       value
#> 1      treatment -0.50000000
#> 2            age  0.30000000
#> 3  baseline_rate  0.05000000
#> 4          shape  1.50000000
#> 5 censoring_rate  0.05612131
#> 
#> Other tables: diagnostics. Read one with as.data.frame(x, what = "diagnostics").
as.data.frame(survival_data, what = "parameters")
#>             term       value
#> 1      treatment -0.50000000
#> 2            age  0.30000000
#> 3  baseline_rate  0.05000000
#> 4          shape  1.50000000
#> 5 censoring_rate  0.05612131
as.data.frame(survival_data, what = "diagnostics")
#>   target_censoring realized_censoring baseline
#> 1              0.3             0.3115  weibull

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 31.1%, 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.

cox_fit <- survival::coxph(survival::Surv(time, status) ~ treatment + age, data = survival_data)
coef(cox_fit)
#>  treatment        age 
#> -0.4636587  0.3135970

The estimates are −0.464 and 0.314, against the generating values −0.5 and 0.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.

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)
#> <simulab_sim:competing_risks> 3 rows x 5 columns
#>   id treatment      time event event_type
#> 1  1         0 11.855911     1    relapse
#> 2  2         0  7.503283     1    relapse
#> 3  3         0  8.405052     1    relapse
#> 
#> Other tables: events. Read one with as.data.frame(x, what = "events").
as.data.frame(observed_times, what = "events")
#>   event_type event_code censoring
#> 1    dropout          0      TRUE
#> 2    relapse          1     FALSE
xtabs(~ treatment + event_type, data = observed_times)
#>          event_type
#> treatment dropout relapse
#>         0      97     153
#>         1      67     183

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} =$ 1.65. Relapse comes first for 183 of the 250 treated and 153 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.

sequences <- simulate_sequences(seed = 1)
sequences
#> <simulab_sim:sequences> 2000 rows x 3 columns
#>    id period     state
#> 1   1      1   Reflect
#> 2   1      2 Elaborate
#> 3   1      3   Analyze
#> 4   1      4 Elaborate
#> 5   1      5     Learn
#> 6   1      6   Reflect
#> 7   1      7 Elaborate
#> 8   1      8     Learn
#> 9   1      9   Reflect
#> 10  1     10   Reflect
#> ... 1990 more rows
#> 
#> Truth (transitions):
#>         from      to probability
#> 1    Reflect Reflect  0.21906379
#> 2      Learn Reflect  0.29843157
#> 3       Plan Reflect  0.30462200
#> 4  Elaborate Reflect  0.31668452
#> 5    Analyze Reflect  0.10757643
#> 6    Reflect   Learn  0.02478853
#> 7      Learn   Learn  0.03837401
#> 8       Plan   Learn  0.08173499
#> 9  Elaborate   Learn  0.14993651
#> 10   Analyze   Learn  0.02386813
#> ... 15 more rows
#> 
#> Other tables: initial_probabilities, wide, settings. Read one with as.data.frame(x, what = "initial_probabilities").

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.

sequence_recovery <- validate_recovery(summarize_transitions(sequences), sequences,
                                       term = "to", estimate = "probability")
head(sequence_recovery)
#>        from count      term   estimate      truth          bias absolute_error
#> 1   Analyze   248   Analyze 0.56750572 0.55310055  0.0144051751   0.0144051751
#> 2   Analyze    54 Elaborate 0.12356979 0.12152212  0.0020476763   0.0020476763
#> 3   Analyze    10     Learn 0.02288330 0.02386813 -0.0009848386   0.0009848386
#> 4   Analyze    88      Plan 0.20137300 0.19393277  0.0074402300   0.0074402300
#> 5   Analyze    37   Reflect 0.08466819 0.10757643 -0.0229082427   0.0229082427
#> 6 Elaborate    18   Analyze 0.05590062 0.07754096 -0.0216403400   0.0216403400
#>   relative_error recovered
#> 1     0.02604441      TRUE
#> 2     0.01685024      TRUE
#> 3    -0.04126165      TRUE
#> 4     0.03836500      TRUE
#> 5    -0.21294852      TRUE
#> 6    -0.27908269      TRUE

The 25 estimated probabilities differ from the generating ones by at most 0.087. 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.

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)
#> <simulab_sim:hmm> 3 rows x 4 columns
#>   id occasion state observation
#> 1  1        1 stuck        idle
#> 2  1        2 stuck        idle
#> 3  1        3 stuck        idle
#> 
#> Truth (transitions):
#>      from      to probability
#> 1 engaged engaged         0.9
#> 2   stuck engaged         0.3
#> 3 engaged   stuck         0.1
#> 4   stuck   stuck         0.7
#> 
#> Other tables: emissions, initial_probabilities. Read one with as.data.frame(x, what = "emissions").
xtabs(~ state + observation, data = hidden)
#>          observation
#> state     idle submit view
#>   engaged  439   2632 1264
#>   stuck   1015    178  472

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().

behavior <- simulate_event_log(groups = 2, actors = 3, n_states = 4,
                               sequence_length = c(5, 8), seed = 1)
head(behavior, 4)
#> <simulab_sim:event_log> 4 rows x 7 columns
#>     group    id   course achievement period state           timestamp
#> 1 Group 1 G1_A1 Course 1      medium      1 Doubt 2020-01-01 00:00:00
#> 2 Group 1 G1_A1 Course 1      medium      2 Doubt 2020-01-01 00:01:10
#> 3 Group 1 G1_A1 Course 1      medium      3 Doubt 2020-01-01 00:05:15
#> 4 Group 1 G1_A1 Course 1      medium      4 Doubt 2020-01-01 00:06:32
#> 
#> Truth (transitions):
#>      group      from        to probability
#> 1  Group 1  Complete  Complete  0.10070835
#> 2  Group 1 Encourage  Complete  0.42327312
#> 3  Group 1  Continue  Complete  0.43868385
#> 4  Group 1     Doubt  Complete  0.28123564
#> 5  Group 1  Complete Encourage  0.79360111
#> 6  Group 1 Encourage Encourage  0.26046548
#> 7  Group 1  Continue Encourage  0.24067761
#> 8  Group 1     Doubt Encourage  0.10338034
#> 9  Group 1  Complete  Continue  0.06142097
#> 10 Group 1 Encourage  Continue  0.03534803
#> ... 22 more rows
#> 
#> Other tables: actors, groups, wide, one_hot. Read one with as.data.frame(x, what = "actors").

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.

network <- simulate_tna(seed = 1)
network
#> State Labels : 
#> 
#>    Analyze, Elaborate, Learn, Plan, Reflect 
#> 
#> Transition Probability Matrix :
#> 
#>              Analyze  Elaborate      Learn      Plan    Reflect
#> Analyze   0.56750572 0.12356979 0.02288330 0.2013730 0.08466819
#> Elaborate 0.05590062 0.10869565 0.13664596 0.3881988 0.31055901
#> Learn     0.04237288 0.44067797 0.04237288 0.2627119 0.21186441
#> Plan      0.22473868 0.04181185 0.07665505 0.3745645 0.28222997
#> Reflect   0.13363029 0.32962138 0.01781737 0.3006682 0.21826281
#> 
#> Initial Probabilities : 
#> 
#>   Analyze Elaborate     Learn      Plan   Reflect 
#>      0.02      0.31      0.09      0.12      0.46 
#> 
#> Generating Transition Probability Matrix (truth) :
#> 
#>           Analyze Elaborate Learn  Plan Reflect
#> Analyze     0.553     0.122 0.024 0.194   0.108
#> Elaborate   0.078     0.070 0.150 0.386   0.317
#> Learn       0.036     0.381 0.038 0.246   0.298
#> Plan        0.221     0.034 0.082 0.358   0.305
#> Reflect     0.113     0.316 0.025 0.327   0.219
#> 
#> Generating Initial Probabilities (truth) :
#> 
#>   Analyze Elaborate     Learn      Plan   Reflect 
#>     0.045     0.334     0.134     0.079     0.407 
#> 
#> Other tables: edges, sequences, wide, model_info. Read one with as.data.frame(x, what = "edges").

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.

regulation <- tna::tna(tna::group_regulation)
known_network <- simulate_tna(n = 2000, transition = regulation$weights,
                              initial = regulation$inits, seed = 1)
known_network
#> State Labels : 
#> 
#>    adapt, cohesion, consensus, coregulate, discuss, emotion, monitor, plan, synthesis 
#> 
#> Transition Probability Matrix :
#> 
#>                   adapt   cohesion  consensus coregulate    discuss    emotion
#> adapt      0.0000000000 0.29629630 0.45634921 0.02777778 0.06349206 0.10846561
#> cohesion   0.0027386541 0.02895149 0.49647887 0.11306729 0.06964006 0.11893584
#> consensus  0.0048441449 0.01337405 0.08203454 0.19323926 0.18723673 0.07424179
#> coregulate 0.0211012199 0.03956479 0.13748764 0.01978239 0.27266733 0.16419387
#> discuss    0.0697995853 0.04302004 0.32066344 0.08517623 0.19765031 0.10711818
#> emotion    0.0016806723 0.32845138 0.32653061 0.03385354 0.09339736 0.07515006
#> monitor    0.0075037519 0.06003002 0.17408704 0.05802901 0.36868434 0.07903952
#> plan       0.0008699435 0.02718573 0.29284472 0.01729013 0.06383210 0.14582427
#> synthesis  0.2176063304 0.03659743 0.47477745 0.05341246 0.05539070 0.08209693
#>               monitor       plan   synthesis
#> adapt      0.02777778 0.01984127 0.000000000
#> cohesion   0.03129890 0.13536776 0.003521127
#> consensus  0.04454507 0.39311289 0.007371525
#> coregulate 0.09033960 0.23738872 0.017474448
#> discuss    0.02384243 0.01209399 0.140635798
#> emotion    0.03577431 0.10180072 0.003361345
#> monitor    0.01850925 0.21860930 0.015507754
#> plan       0.07394519 0.37635929 0.001848630
#> synthesis  0.01582591 0.06429278 0.000000000
#> 
#> Initial Probabilities : 
#> 
#>      adapt   cohesion  consensus coregulate    discuss    emotion    monitor 
#>     0.0095     0.0650     0.2355     0.0165     0.1710     0.1355     0.1380 
#>       plan  synthesis 
#>     0.2085     0.0205 
#> 
#> Generating Transition Probability Matrix (truth) :
#> 
#>            adapt cohesion consensus coregulate discuss emotion monitor  plan
#> adapt      0.000    0.273     0.477      0.022   0.059   0.120   0.033 0.016
#> cohesion   0.003    0.027     0.498      0.119   0.060   0.116   0.033 0.141
#> consensus  0.005    0.015     0.082      0.188   0.188   0.073   0.047 0.396
#> coregulate 0.016    0.036     0.135      0.023   0.274   0.172   0.086 0.239
#> discuss    0.071    0.048     0.321      0.084   0.195   0.106   0.022 0.012
#> emotion    0.002    0.325     0.320      0.034   0.102   0.077   0.036 0.100
#> monitor    0.011    0.056     0.159      0.058   0.375   0.091   0.018 0.216
#> plan       0.001    0.025     0.290      0.017   0.068   0.147   0.076 0.374
#> synthesis  0.235    0.034     0.466      0.044   0.063   0.071   0.012 0.075
#>            synthesis
#> adapt          0.000
#> cohesion       0.004
#> consensus      0.008
#> coregulate     0.019
#> discuss        0.141
#> emotion        0.003
#> monitor        0.016
#> plan           0.002
#> synthesis      0.000
#> 
#> Generating Initial Probabilities (truth) :
#> 
#>      adapt   cohesion  consensus coregulate    discuss    emotion    monitor 
#>      0.012      0.060      0.214      0.019      0.176      0.152      0.144 
#>       plan  synthesis 
#>      0.204      0.020 
#> 
#> Other tables: edges, sequences, wide, model_info. Read one with as.data.frame(x, what = "edges").

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

network_recovery <- validate_recovery(known_network, known_network,
                                      term = "to", estimate = "weight")
summarize_simulations(network_recovery, variables = "absolute_error")
#>         variable observations        mean          sd minimum    maximum
#> 1 absolute_error           81 0.003937092 0.004688553       0 0.02321182

Over the 81 edges, the mean absolute error is 0.004 and the largest is 0.023. 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.

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)
#> <simulab_sim:network> 3 rows x 3 columns
#>   from to weight
#> 1    1  4      1
#> 2    3  4      1
#> 3    1  5      1
#> 
#> Truth (settings):
#>       model directed loops weight nodes edges
#> 1 bernoulli    FALSE FALSE binary    15    26
#> 
#> Other tables: nodes, adjacency. Read one with as.data.frame(x, what = "nodes").
as.data.frame(classroom, what = "settings")
#>       model directed loops weight nodes edges
#> 1 bernoulli    FALSE FALSE binary    15    26

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 26.

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().

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
#>   variable observations mean       sd minimum maximum
#> 1    value          200 3.97 3.471217       2      21

With attachment = 2, each new node adds two edges to existing nodes. The degree ranges from 2 to 21, with a mean of 3.97.

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.

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)
#> <simulab_sim:temporal_network> 4 rows x 6 columns
#>     from     to onset terminus weight censored
#> 1 Node 2 Node 1     5        7      1    FALSE
#> 2 Node 5 Node 1     1        3      1    FALSE
#> 3 Node 5 Node 1     7       10      1    FALSE
#> 4 Node 8 Node 1     1        3      1    FALSE
#> 
#> Truth (settings):
#>   periods directed loops initial_probability formation_probability
#> 1      12     TRUE FALSE                 0.1                  0.05
#>   dissolution_probability weight
#> 1                     0.3 binary
#> 
#> Other tables: events, snapshots, nodes. Read one with as.data.frame(x, what = "events").
components(temporal)
#>       table rows columns
#> 1      data   49       6
#> 2    events   83       5
#> 3 snapshots  128       4
#> 4     nodes   10       2
#> 5  settings    1       7

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().

layers <- simulate_multiplex_network(
  nodes = 10,
  layers = c("advice", "friendship"),
  probability = c(0.15, 0.3),
  seed = 1
)
as.data.frame(layers, what = "layers")
#>        layer probability
#> 1     advice        0.15
#> 2 friendship        0.30
xtabs(~ layer, data = layers)
#> layer
#>     advice friendship 
#>         12         30

The network has 10 nodes in two directed layers without self-loops. Each layer has 90 ordered pairs, so the expected number of edges is 13.5 for advice and 27 for friendship. The realized numbers are 12 and 30. 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