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

## ----library------------------------------------------------------------------
library(simulab)

## ----regression-drawn---------------------------------------------------------
drawn <- simulate_regression(seed = 1)
drawn

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

## ----regression-recovery------------------------------------------------------
fit <- lm(outcome ~ x1 + x2, data = known)
validate_recovery(fit, known)

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

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

## ----batch-values, include = FALSE--------------------------------------------
largest_bias <- max(abs(performance$mean))
empirical_se <- mean(performance$sd)

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

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

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

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

## ----distributions------------------------------------------------------------
head(list_distributions())

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

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

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

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

## ----propensity-values, include = FALSE---------------------------------------
age_exposed <- subset(age_by_exposure, exposed == 1)$mean
age_unexposed <- subset(age_by_exposure, exposed == 0)$mean

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

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

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

## ----ar1-values, include = FALSE----------------------------------------------
lag_one <- off_diagonal(sample_correlation, 1)
lag_two <- off_diagonal(sample_correlation, 2)
lag_three <- off_diagonal(sample_correlation, 3)

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

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

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

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

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

## ----factors-values, include = FALSE------------------------------------------
primary_loadings <- apply(abs(estimated_loadings), 1, max)
cross_loadings <- apply(abs(estimated_loadings), 1, min)

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

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

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

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

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

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

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

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

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

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

## ----growth-values, include = FALSE-------------------------------------------
mean_first <- subset(mean_trajectory, time == 0)$outcome
mean_last <- subset(mean_trajectory, time == 4)$outcome

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

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

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

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

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

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

## ----sequences----------------------------------------------------------------
sequences <- simulate_sequences(seed = 1)
sequences

## ----sequences-recovery-------------------------------------------------------
sequence_recovery <- validate_recovery(summarize_transitions(sequences), sequences,
                                       term = "to", estimate = "probability")
head(sequence_recovery)

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

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

## ----tna, eval = has_tna------------------------------------------------------
network <- simulate_tna(seed = 1)
network

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

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

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

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

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

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

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

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