Analysis of Change: Fit Change Models to One or Many Trajectories

Fits a change model to longitudinal data: any of the package's four nonlinear change models (negative exponential, logistic, Gompertz, Richards; the parameterizations of Kelley, 2005, 2008, exactly as generated by the simulate_longitudinal_*() simulators) or a polynomial of any order (the linear change model of simulate_longitudinal_polynomial). Two estimation methods are offered. The default, method = "two_stage", fits each unit's curve separately, using only that unit's data, and then summarizes the unit-level parameters: their mean, their standard deviation and variance across units (the individual differences), and the standard error of the mean. method = "mixed" fits the proper random-coefficients mixed-effects model simultaneously, with every parameter carrying a random effect, so the reported between-unit spread is a variance component purged of estimation noise. With a single trajectory (\(N = 1\), or id = NULL) the two-stage method reduces to one least squares fit of that unit's change, reported with its standard errors.

Usage

analysis_of_change(
  data,
  id,
  time,
  outcome,
  model = c("negative_exponential", "logistic", "gompertz", "richards", "polynomial"),
  method = c("two_stage", "mixed"),
  order = 1L,
  start = NULL,
  maxiter = 500L
)

Arguments

data: A data.frame in long format: one row per observation, as returned by the simulate_longitudinal_*() simulators.
id: Name of the column identifying units (persons, animals, trees, classrooms). NULL treats all rows as a single trajectory (method = "two_stage" only).
time: Name of the time column.
outcome: Name of the outcome column.
model: Which change model to fit: "negative_exponential" (parameters alpha, zeta, gamma), "logistic" (alpha, beta, gamma, zeta), "gompertz" (alpha, beta, gamma, zeta), "richards" (alpha, beta, gamma, delta, zeta), or "polynomial" (coefficients b0 through bP, with the order \(P\) set by order). See the corresponding simulator's help page for what each parameter means.
method: How the model is estimated. "two_stage" (the default) fits one curve per unit, each from that unit's data alone, and summarizes across units; it is transparent, works at \(N = 1\), and never lets one unit's data influence another's fit. "mixed" estimates the random-coefficients model simultaneously: lme4::lmer() for the polynomial (a linear mixed model) and nlme::nlme() for the nonlinear curves (started at the two-stage estimates), with uncorrelated random effects on every parameter. See Details for how to choose.
order: Polynomial order \(P\) when model = "polynomial" (1 is straight-line change, 2 quadratic, and so on). Default 1. Ignored for the nonlinear models.
start: Optional named numeric vector of starting values for a nonlinear model's parameters, used for every unit (and, under method = "mixed", for the fixed effects). NULL (default) derives data-driven starting values per unit (see Details). Ignored for the polynomial, whose fit is closed form.
maxiter: Maximum number of iterations passed to nls (and, under method = "mixed", to the nonlinear mixed-effects optimizer). Default 500.

Value

A data.frame with one row per model parameter and columns

term: The parameter name.
estimate: The mean of the unit-level estimates (two-stage; with one trajectory, that unit's estimate) or the fixed effect (mixed).
se: The standard error of estimate: the between-unit standard deviation over the square root of the number of fitted units (two-stage), the asymptotic standard error of the single fit (one trajectory), or the fixed-effect standard error (mixed).
sd_units: The between-unit standard deviation of the parameter: the spread of the unit-level estimates (two-stage; NA for a single trajectory) or the random-effect standard deviation (mixed).
var_units: The corresponding variance.

Attributes: "model", "method", "n_units" (units supplied), "n_used" (units whose fit converged; equal to "n_units" under method = "mixed"), "sigma" (the level-one residual standard deviation), and "per_unit_estimates" (a units-by-parameters matrix: the separate unit-level estimates under two-stage, or the unit-level predictions coef() under mixed, which are shrunken toward the fixed effects).

Details

Choosing between the methods. The two-stage (curve-by-curve) route is the transparent classic: every unit's curve is inspectable, no unit's data influence another's fit, and it is the only method available for a single trajectory. Its known cost is that sd_units reflects the spread of estimates, which adds each fit's estimation noise to the true individual differences; with short or noisy trajectories it therefore overstates the population standard deviation. The mixed method estimates that between-unit variation as a variance component, separating it from level-one error, and borrows strength across units, at the price of a harder estimation problem (and, for the nonlinear curves, occasional convergence trouble; the fit is started at the two-stage estimates, and a failure suggests simplifying the model or falling back to two-stage). Under method = "mixed" the random effects are uncorrelated across parameters, matching the simulators' default.

Relation to existing tools. The two-stage method is the idea behind nlme::lmList() and nlme::nlsList(), and the mixed method wraps lme4::lmer() and nlme::nlme(); base R also ships self-starting curves (SSasymp, SSfpl, SSgompertz) in other parameterizations. What this function adds is the package's landmark parameterizations (the intercept-shifting \(\zeta\) floor, the Richards \(\delta\); SSgompertz's \(a \exp(-b_2 b_3^x)\) answers no substantive question directly), the exact match to the simulate_longitudinal_*() simulators so design studies close the loop, one interface across linear and nonlinear change, and the package's tidy summary with failed fits dropped under a single counted warning.

Units whose two-stage fit does not converge are dropped with a single warning reporting how many, and the effective count is the "n_used" attribute. If no unit's fit converges the function stops. The polynomial fit is closed form and does not fail on any trajectory with at least \(P + 1\) occasions.

Starting values. Unless start is supplied, each unit's nonlinear starting values are derived from that unit's data: the floor and the span from early and late observations, the inflection time from where the trajectory crosses the middle of its range, the curvature from the time the trajectory needs to travel the central half of its range, and, for the Richards model, a logistic start (delta = 1). Cleanly measured trajectories rarely need more; hard cases (very short series, strong decreasing curves, near-flat change) may need an explicit start.

References

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.

Pinheiro, J. C., & Bates, D. M. (2000). Mixed-effects models in S and S-PLUS. Springer.

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

Author

Ken Kelley kkelley@nd.edu

Examples

# Simulate a Gompertz population with individual differences, then
# recover both the mean curve and the spread of its parameters.
set.seed(113)
d <- simulate_longitudinal_gompertz(
  n = 40, target_times = 0:10,
  fixed_parameters = c(alpha = 75, beta = 3, gamma = 0.55, zeta = 10),
  random_variances = c(alpha = 25, beta = 0.4, gamma = 0.005, zeta = 4),
  error_variance = 4
)
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "gompertz")
#>  term  estimate se     sd_units var_units
#>  alpha 75       1.06   6.71     45       
#>  beta  3.08     0.107  0.674    0.454    
#>  gamma 0.557    0.0138 0.0875   0.00766  
#>  zeta  10       0.529  3.35     11.2     

# The same data as a straight-line (order 1) polynomial: the linear
# model has nothing to say about floors, ceilings, or timing.
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "polynomial", order = 1)
#>  term estimate se   sd_units var_units
#>  b0   12.9     1.24 7.86     61.8     
#>  b1   8.13     0.14 0.887    0.787    

# A single unit's trajectory: N = 1 is one fit, with asymptotic
# standard errors in place of between-unit spread.
one <- d[d$id == levels(d$id)[1], ]
analysis_of_change(one, id = NULL, time = "time", outcome = "y",
                   model = "gompertz")
#>  term  estimate se     sd_units var_units
#>  alpha 73.8     2.65   <NA>     <NA>     
#>  beta  3.11     0.116  <NA>     <NA>     
#>  gamma 0.72     0.0732 <NA>     <NA>     
#>  zeta  9.82     1.85   <NA>     <NA>     

# The proper mixed-effects model, fit simultaneously across units:
# nlme::nlme() for a nonlinear curve, lme4::lmer() for the polynomial
# (the latter requires lme4 to be installed).
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "gompertz", method = "mixed")
#>  term  estimate se     sd_units var_units
#>  alpha 75.1     1.03   5.64     31.8     
#>  beta  3.06     0.113  0.698    0.487    
#>  gamma 0.549    0.0131 0.0646   0.00417  
#>  zeta  9.8      0.448  1.87     3.51     
analysis_of_change(d, id = "id", time = "time", outcome = "y",
                   model = "polynomial", order = 1,
                   method = "mixed")
#>  term estimate se    sd_units var_units
#>  b0   12.9     1.19  5.97     35.6     
#>  b1   8.13     0.135 0.339    0.115