Nothing
Simulating a Wide Variety of Data Shapes, Sizes and Distributions
In simulab: Simulates a Wide Variety of Data Shapes, Sizes and Distributions
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.
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
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
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)
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.
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
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.
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.
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.
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().
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.
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")
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.
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.
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
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.
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)
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.
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)
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.
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)
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.
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.
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.
factor_fit <- factanal(subset(factor_data, select = -id), factors = 2, rotation = "promax")
estimated_loadings <- round(unclass(loadings(factor_fit)), 2)
estimated_loadings
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.
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
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.
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
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.
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)
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.
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().
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).
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.
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
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.
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.
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")
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.
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.
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)
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.
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.
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.
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().
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.
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.
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.
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.
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")
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().
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.
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().
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)
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 |
Try the simulab package in your browser
Any scripts or data that you put into this service are public.
simulab documentation built on Oct. 15, 2026, 5:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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.
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
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
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)
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.
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
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.
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.
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.
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().
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.
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")
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.
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.
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
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.
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)
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.
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)
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.
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)
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.
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.
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.
factor_fit <- factanal(subset(factor_data, select = -id), factors = 2, rotation = "promax") estimated_loadings <- round(unclass(loadings(factor_fit)), 2) estimated_loadings
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.
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
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.
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
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.
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)
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.
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().
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).
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.
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
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.
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.
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")
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.
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.
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)
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.
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.
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.
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().
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.
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.
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.
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.
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")
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().
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.
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().
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)
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 |
Try the simulab package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.