simulate_longitudinal_richards: Simulate Data From a Richards Change (Growth) Model

Description

Generates longitudinal data from a random-coefficients Richards change model, the flexible sigmoidal family whose point of inflection is itself a parameter rather than a fixed fraction of total change (Richards, 1959). The deterministic part of the five parameter Richards curve is $$\mu(t) = \frac{\alpha}{\left(1 + \delta \exp(-\gamma (t - \beta))\right)^{1/\delta}} + \zeta,$$ the parameterization of Kelley (2005, 2008). The Richards family subsumes the package's other sigmoidal curves as special cases: at $\delta = 1$ it is exactly the logistic, and in the limit $\delta \rightarrow 0$ it is the Gompertz, so $\delta$ lets the data, or the theory, choose where in its course the process turns from acceleration to deceleration. Guo, Cheng, and Kelley (2016) use this flexibility to model self-replicating malware propagation, where the network structure moves the inflection of the outbreak.

Usage

simulate_longitudinal_richards(
  n,
  target_times = NULL,
  fixed_parameters,
  time_range = NULL,
  occasions = NULL,
  time_distribution = "uniform",
  random_variances = 0,
  random_correlation = NULL,
  error_variance = NULL,
  reliability = NULL,
  error_structure = c("independent", "ar1", "compound_symmetry", "toeplitz"),
  error_correlation = NULL,
  timing_sd = 0
)

Value

A long-format data.frame with columns id,

population, occasion, target_time, time,

true_score, and y, directly usable with

plot_trajectories and nonlinear mixed-model fitters such as nlme::nlme(). Attributes carry the model,

fixed_parameters, random_covariance,

error_variance, error_covariance,

reliability_by_occasion, and schedule

("shared" or "unit_specific").

Arguments

n

A single positive integer, the number of units (persons, animals, trees, classrooms) whose trajectories are drawn, or a vector giving the number of units for each population, one entry per parameter vector in fixed_parameters.

target_times

Numeric vector of the nominal measurement times, one shared schedule for every unit. Give either this or time_range.

fixed_parameters

The population parameters c(alpha = , beta = , gamma = , delta = , zeta = ) (an unnamed length-5 vector is taken in that order), or a list of such vectors, one per population (distinct data generating parameter vectors):

alpha: The total change: the curve travels alpha units from the lower asymptote zeta to the upper asymptote $\alpha + \zeta$ (for $\gamma > 0$).

beta

The point of inflection on the time axis.

gamma

The curvature: how sharply the curve rises through its inflection. Negative values flip the curve to decreasing.

delta

The shape parameter ($\delta > 0$): where the inflection falls on the outcome axis, $y^{*} = \alpha (1 + \delta)^{-1/\delta} + \zeta$. At $\delta = 1$ the inflection sits at half the total change (the logistic); as $\delta \rightarrow 0$ it slides down to $\alpha / e + \zeta$, about 36.8% (the Gompertz); larger $\delta$ pushes it later than halfway.

zeta

The lower asymptote (for $\gamma > 0$). The intercept is $\phi = \alpha (1 + \delta \exp(\gamma \beta))^{-1/\delta} + \zeta$.

time_range

Alternative to target_times: c(lower, upper) bounds from which each unit draws its own measurement times, so no two units share a schedule (e.g., age in weeks at testing rather than a fixed grade). Requires occasions; with time_range, the level-one error must be a single error_variance with the default independent structure, and timing_sd does not apply.

occasions

With time_range: a single positive integer (every unit measured the same number of times) or c(min, max), from which each unit's number of measurement times is drawn uniformly.

time_distribution

Distribution of the unit-specific times over time_range; currently "uniform".

random_variances

Between-unit variances of (alpha, beta, gamma, delta, zeta), a single number recycled to all five or a length-5 vector. Default 0. A named vector over any subset of the parameter names (e.g., c(delta = 0.04)) varies only those named and leaves the rest fixed. A positive variance on delta must be small enough that no unit's drawn $\delta$ falls at or below zero, where the curve is undefined; the function stops with a count if that happens.

random_correlation

Optional 5-by-5 correlation matrix among the random parameters; default uncorrelated.

error_variance, reliability, error_structure, error_correlation, timing_sd

The level-one error and assessment-time machinery, with the same meaning as in simulate_longitudinal_polynomial: specify exactly one of error_variance or reliability (solved here through the first-order delta method true-score variance; because the curve is nonlinear in its parameters, that approximation, and the reliability_by_occasion attribute with it, can drift from the realized variance ratio when the random variances are large relative to the mean curve); error_structure and error_correlation set the across-occasion error correlation; timing_sd jitters the actual assessment times.

Author

Ken Kelley [email protected]

Details

The logistic and the Gompertz fix where the inflection falls as a fraction of total change (50% and 36.8%); choosing between them is choosing that fraction by assumption. The Richards curve makes the fraction estimable through $\delta$, at the price of one more parameter and a harder estimation problem, since $\delta$ and $\gamma$ carry overlapping information in finite samples (Richards, 1959; Kelley, 2005). Simulating from the Richards family at several $\delta$ values is the natural way to study whether a design can tell those shapes apart.

References

Guo, H., Cheng, H. K., & Kelley, K. (2016). Impact of network structure on malware propagation: A growth curve perspective. Journal of Management Information Systems, 33(1), 296--325.

Kelley, K. (2005). Estimating nonlinear change models in heterogeneous populations when class membership is unknown: Defining and developing the latent classification differential change model (Doctoral dissertation). University of Notre Dame.

Kelley, K. (2008). Nonlinear change models in populations with unobserved heterogeneity. Methodology, 4(3), 97--112.

Richards, F. J. (1959). A flexible growth function for empirical use. Journal of Experimental Botany, 10(2), 290--301.

Examples

Run this code

# The family in one panel: five Richards curves sharing alpha, beta,
# gamma, and zeta and differing only in the shape parameter delta.
# delta near 0 is the Gompertz, delta = 1 is the logistic, and
# larger delta pushes the inflection later than halfway. Each curve
# is its own population of size one, so a single call draws the whole panel.
panel <- simulate_longitudinal_richards(
  n = 1, target_times = seq(0, 12, by = 0.1),
  fixed_parameters = list(
    c(alpha = 80, beta = 6, gamma = 0.9, delta = 0.02, zeta = 10),
    c(alpha = 80, beta = 6, gamma = 0.9, delta = 0.25, zeta = 10),
    c(alpha = 80, beta = 6, gamma = 0.9, delta = 1.00, zeta = 10),
    c(alpha = 80, beta = 6, gamma = 0.9, delta = 3.00, zeta = 10),
    c(alpha = 80, beta = 6, gamma = 0.9, delta = 8.00, zeta = 10)),
  error_variance = 0
)
plot_trajectories(panel, id = "id", time = "time",
                  outcome = "true_score", group = "population")

# Individual differences in a single parameter: only the shape
# varies (a named entry leaves the other variances at zero), so the
# curves agree on floor, ceiling, timing, and curvature yet turn
# from acceleration to deceleration at different heights.
set.seed(113)
d_delta <- simulate_longitudinal_richards(
  n = 25, target_times = seq(0, 12, by = 0.5),
  fixed_parameters = c(alpha = 80, beta = 6, gamma = 0.9,
                       delta = 1, zeta = 10),
  random_variances = c(delta = 0.04), error_variance = 0
)
plot_trajectories(d_delta, id = "id", time = "time",
                  outcome = "true_score")

# Individual differences in every parameter, plus level-one error:
# a late-inflecting outbreak-style curve, with delta = 3 placing the
# inflection at about 63% of total change, (1 + 3)^(-1/3).
set.seed(113)
d <- simulate_longitudinal_richards(
  n = 30, target_times = 0:12,
  fixed_parameters = c(alpha = 80, beta = 6, gamma = 0.9,
                       delta = 3, zeta = 10),
  random_variances = c(alpha = 36, beta = 1, gamma = 0.01,
                       delta = 0.04, zeta = 9),
  error_variance = 16
)
plot_trajectories(d, id = "id", time = "time", outcome = "y")

# delta = 1 reproduces the logistic exactly: with no randomness and
# no error, the two simulators return identical true scores.
d_r <- simulate_longitudinal_richards(
  n = 1, target_times = 0:5,
  fixed_parameters = c(alpha = 80, beta = 3, gamma = 1,
                       delta = 1, zeta = 10),
  error_variance = 0)
d_l <- simulate_longitudinal_logistic(
  n = 1, target_times = 0:5,
  fixed_parameters = c(alpha = 80, beta = 3, gamma = 1, zeta = 10),
  error_variance = 0)
all.equal(d_r$true_score, d_l$true_score)

Run the code above in your browser using DataLab