Spline Simulation with `simglm`
In simglm: Simulate Models Based on the Generalized Linear Model

library(knitr)
library(dplyr)
library(simglm)
library(splines)
knit_print.data.frame = function(x, ...) {
  res = paste(c('', '', kable(x, output = FALSE)), collapse = '\n')
  asis_output(res)
}

Spline Simulation

Spline terms can be specified directly in the formula argument for fixed-effect simulation. The supported spline basis functions are ns() and bs(), both from the splines package. These terms are expanded into generated basis columns during simulation while the original variable is kept in the simulated data so the same formula can be used for model fitting.

Natural Cubic Spline

The example below simulates two fixed variables, x1 and x2, and uses a natural cubic spline for x2. The formula uses regular R formula syntax.

set.seed(321)

sim_arguments <- list(
  formula = y ~ 1 + x1 + ns(x2, df = 4),
  fixed = list(
    x1 = list(var_type = 'continuous', mean = 0, sd = 1),
    x2 = list(var_type = 'continuous', mean = 0, sd = 1)
  ),
  sample_size = 10
)

simulate_fixed(data = NULL, sim_arguments)

The generated design matrix contains one column for the intercept, one column for x1, four spline basis columns for x2, and the original x2 variable. The spline basis columns are named x2_1, x2_2, x2_3, and x2_4.

Generating an Outcome

Regression weights can be specified with a named list. This is useful for spline terms because one formula term expands to multiple basis columns. The name used for the spline term should match the formula term.

set.seed(321)

sim_arguments <- list(
  formula = y ~ 1 + x1 + ns(x2, df = 4),
  fixed = list(
    x1 = list(var_type = 'continuous', mean = 0, sd = 1),
    x2 = list(var_type = 'continuous', mean = 0, sd = 1)
  ),
  error = list(variance = 1),
  sample_size = 100,
  reg_weights = list(
    `(Intercept)` = 0,
    x1 = 0.25,
    `ns(x2, df = 4)` = c(0.75, 0.5, 0, 0)
  )
)

sim_data <- simulate_fixed(data = NULL, sim_arguments) |>
  simulate_error(sim_arguments) |>
  generate_response(sim_arguments)

head(sim_data)

B-Spline Basis

The same pattern can be used for a B-spline basis with bs().

set.seed(321)

sim_arguments_bs <- list(
  formula = y ~ 1 + x1 + bs(x2, df = 4),
  fixed = list(
    x1 = list(var_type = 'continuous', mean = 0, sd = 1),
    x2 = list(var_type = 'continuous', mean = 0, sd = 1)
  ),
  error = list(variance = 1),
  sample_size = 100,
  reg_weights = list(
    `(Intercept)` = 0,
    x1 = 0.25,
    `bs(x2, df = 4)` = c(0.75, 0.5, 0, 0)
  )
)

bs_data <- simulate_fixed(data = NULL, sim_arguments_bs) |>
  simulate_error(sim_arguments_bs) |>
  generate_response(sim_arguments_bs)

head(bs_data)

Model Fitting

Spline models can be fit with the existing model_fit argument. The original variable is retained in the simulated data, so the fitted model can use the same formula syntax with ns() or bs().

set.seed(321)

sim_arguments_fit <- list(
  formula = y ~ 1 + x1 + ns(x2, df = 4),
  fixed = list(
    x1 = list(var_type = 'continuous', mean = 0, sd = 1),
    x2 = list(var_type = 'continuous', mean = 0, sd = 1)
  ),
  error = list(variance = 1),
  sample_size = 100,
  reg_weights = list(
    `(Intercept)` = 0,
    x1 = 0.25,
    `ns(x2, df = 4)` = c(0.75, 0.5, 0, 0)
  ),
  model_fit = list(model_function = 'lm')
)

fit <- simulate_fixed(data = NULL, sim_arguments_fit) |>
  simulate_error(sim_arguments_fit) |>
  generate_response(sim_arguments_fit) |>
  model_fit(sim_arguments_fit)

broom::tidy(fit)

Power Simulation

Spline models can also be used with the existing power simulation workflow. Each spline basis coefficient is returned as its own model term in the fitted model, for example ns(x2, df = 4)1 through ns(x2, df = 4)4. Alternative power thresholds should be named by the fitted model terms returned by broom::tidy().

set.seed(321)

sim_arguments_power <- list(
  formula = y ~ 1 + x1 + ns(x2, df = 4),
  fixed = list(
    x1 = list(var_type = 'continuous', mean = 0, sd = 1),
    x2 = list(var_type = 'continuous', mean = 0, sd = 1)
  ),
  error = list(variance = 1),
  sample_size = 80,
  reg_weights = list(
    `(Intercept)` = 0,
    x1 = 0.25,
    `ns(x2, df = 4)` = c(0.75, 0.5, 0, 0)
  ),
  model_fit = list(model_function = 'lm'),
  power = list(
    thresholds = list(
      x1 = 0.25,
      `ns(x2, df = 4)1` = c(0.5, 0.75),
      `ns(x2, df = 4)2` = 0.5
    )
  ),
  replications = 10,
  extract_coefficients = TRUE
)

power_out <- replicate_simulation(sim_arguments_power)

compute_statistics(
  power_out,
  sim_arguments_power,
  type_1_error = FALSE,
  alternative_power = TRUE
)

Random effects can be included with spline fixed effects using the existing random-effect syntax. For example, y ~ 1 + ns(time, df = 4) + (1 | id) adds a random intercept while keeping the spline as a fixed-effect term.