Fit a longitudinal joint regression model — gamlss_longitudinal • gamlss.longitudinal

This function fits a longitudinal model to a dataset with gamlss margins and copula fit to dependence. Any linear or factor covariates can be fit to any parameters of the copula or margin distributions. The model is fit using RS() optimisation and the joint likelihood by default. Select use dlcopdpar=FALSE to fit separately optimised models for the margin and copula likelihoods which can be quicker with a slight loss to overall fit.

Usage

gamlss_longitudinal(
  dataset,
  margin_dist,
  copula_dist,
  time_var = NA,
  subject_var = NA,
  mu.formula = ("response ~ 1"),
  sigma.formula = ("~ 1"),
  nu.formula = ("~ 1"),
  tau.formula = ("~ 1"),
  theta.formula = ("~ 1"),
  zeta.formula = ("~ 1"),
  include_dlcopdpar = TRUE,
  check_dlcopdpar_gradient = FALSE,
  inner_stop_crit = NA,
  outer_stop_crit = NA,
  start_step_size = 0.5,
  step_adjustment = NA,
  max_steps = 5,
  start_from = NA,
  warm_start_joint = TRUE,
  warm_start_joint_iter = 5,
  verbose = 1,
  plot_results = FALSE,
  true_val = NA,
  method = "RS",
  max_outer_iter = 100,
  max_inner_iter = 100,
  max_negative_outer_streak = 10,
  max_elapsed_sec = Inf,
  use_backtracking = TRUE,
  backtracking_max_halves = 50,
  cg_max_stall = 5,
  cg_max_delta = 0.5,
  cg_armijo_c1 = 1e-04,
  cg_grad_tol = NA,
  cg_step_tol = NA,
  cg_update_lambda = TRUE,
  cg_lambda_update_every = 10,
  cg_max_lambda_updates = NA,
  cg_raw_loglik_drop_tol = 10,
  cg_line_search = "best",
  cg_max_line_search_evals = 60,
  cg_gradient_method = "forward",
  discrete_score_method = c("analytical", "finite"),
  cg_zeta_hessian = "analytical",
  cg_hessian_method = c("analytical", "finite", "auto"),
  compute_vcov = TRUE,
  vcov_method = c("analytical", "numderiv"),
  vcov_numderiv = FALSE,
  use_Rcpp = FALSE,
  lambda_start = NA,
  lambda_penalty_K = 2,
  rs_update_lambda = TRUE,
  rs_smooth_trust_radius = Inf
)

Arguments

dataset: Long-format data frame containing the response, subject, time, and covariate columns.
margin_dist: Marginal distribution specified as a gamlss family object, e.g. GA(), NO(), PO(), NBI(), etc.
copula_dist: Copula distribution code, one of "N", "C", "F", "G", "J", or "t".
time_var: Name of the time variable in dataset.
subject_var: Name of the subject identifier in dataset.
mu.formula: Formula for the mean parameter of the marginal distribution
sigma.formula: Formula for the sigma parameter of the marginal distribution
nu.formula: Formula for the nu parameter of the marginal distribution
tau.formula: Formula for the tau parameter of the marginal distribution
theta.formula: Formula for the theta parameter of the copula distribution
zeta.formula: Formula for the zeta parameter of the copula distribution
include_dlcopdpar: Include the derivative of the copula likelihood with respect to the margin parameters in the joint likelihood.
check_dlcopdpar_gradient: If TRUE, run an optional finite-difference diagnostic for the margin score contribution when include_dlcopdpar = TRUE.
inner_stop_crit: Stopping criterion for the inner loop. If NA or NULL, an automatic data-adaptive value is used.
outer_stop_crit: Stopping criterion for the outer loop. If NA or NULL, an automatic data-adaptive value is used.
start_step_size: Initial step size for the backfitting algorithm
step_adjustment: Step size adjustment factor
max_steps: Maximum number of times for reducing the step size
start_from: Starting values for the parameters if needed
warm_start_joint: Logical; if TRUE (default), RS joint fits started without explicit start_from first run a short separate RS stabilisation phase and use those coefficients as the joint starting values.
warm_start_joint_iter: Integer; number of separate RS outer iterations used for the default joint warm start.
verbose: Level of output to the console 3 = ALL, 0 = Minimal
plot_results: Plot the results of the optimisation
true_val: True values for the parameters if known for plotting
method: Optimisation method to use, RS() is the default
max_outer_iter: Maximum number of outer iterations
max_inner_iter: Maximum number of inner iterations
max_negative_outer_streak: Maximum number of consecutive negative outer log-likelihood changes allowed before stopping.
max_elapsed_sec: Optional maximum elapsed fitting time in seconds. If finite, the optimiser stops with an error once this budget is exceeded.
use_backtracking: Logical; if TRUE (default), apply step-halving backtracking to reject downhill inner updates.
backtracking_max_halves: Integer; maximum number of consecutive step halvings attempted after a rejected update before taking no step.
cg_max_stall: Integer; for method = "CG" only. Maximum number of consecutive outer iterations where no improving step is found before CG stops.
cg_max_delta: Numeric; for method = "CG" only. Maximum absolute coefficient step size used to limit Newton/trust-region updates.
cg_armijo_c1: Numeric; for method = "CG" only. Minimum improvement threshold used by the line-search acceptance rule.
cg_grad_tol: Numeric; for method = "CG" only. Penalized-gradient infinity-norm convergence tolerance. If NA, selected from outer_stop_crit.
cg_step_tol: Numeric; for method = "CG" only. Accepted-step L2 convergence tolerance. If NA, selected from outer_stop_crit.
cg_update_lambda: Logical; for method = "CG" only. If TRUE, update smoother penalties during CG iterations.
cg_lambda_update_every: Integer; for method = "CG" only. When cg_update_lambda = TRUE, update each smoother's lambda every this many outer iterations. Use 1 to update every CG iteration.
cg_max_lambda_updates: Integer; for method = "CG" only. Maximum number of smoother penalty update rounds. Use NA for no cap.
cg_raw_loglik_drop_tol: Numeric; for method = "CG" only. Stop CG as not converged if the raw joint log-likelihood drops this far below the best raw joint log-likelihood seen after at least one lambda update. Use NA to disable.
cg_line_search: Character; for method = "CG" only. "best" evaluates candidate steps up to cg_max_line_search_evals before taking the largest improvement, while "first" accepts the first improving candidate step.
cg_max_line_search_evals: Integer; for method = "CG" only. Optional cap on the number of candidate likelihood evaluations per outer iteration.
cg_gradient_method: Character; for method = "CG" only. "analytical" uses the same score components as RS, "forward" uses one-sided finite differences, and "central" uses two-sided finite differences.
discrete_score_method: Character. For discrete margins using exact rectangle likelihoods, choose "analytical" for vectorised rectangle-score assembly or "finite" for slow row-wise finite-difference scores.
cg_zeta_hessian: Character; for method = "CG" only. "analytical" uses the analytical Hessian for the zeta block, while "finite" replaces the zeta-zeta block with central finite differences of the raw joint log-likelihood.
cg_hessian_method: Character; for method = "CG" only. "analytical" uses the semi-analytical Hessian for Newton steps, "finite" uses a full finite-difference Hessian, and "auto" tries analytical then falls back to finite differences when needed.
compute_vcov: Logical; if TRUE (default), compute and store the model variance-covariance output at the end of fitting.
vcov_method: Character; fit-time vcov method when compute_vcov = TRUE. One of "analytical" or "numderiv". Analytical vcov falls back to the numerical reference path if the analytical Hessian cannot be inverted.
vcov_numderiv: Logical; passed to vcov.gamlss.longitudinal() when compute_vcov = TRUE.
use_Rcpp: Use Rcpp for matrix operations
lambda_start: Optional starting value for smooth-term penalties.
lambda_penalty_K: Penalty strength used when updating smooth-term smoothing parameters.
rs_update_lambda: Logical; for method = "RS" only. If TRUE, update smoothing parameters by the RS GAIC step; if FALSE, keep lambda_start fixed.
rs_smooth_trust_radius: Numeric; for method = "RS" only. Optional L2 trust radius applied separately to each smooth coefficient block after the RS weighted least-squares proposal. Use Inf to disable.

Details

Formula inputs are converted to fixed-effect and smooth model matrices by create_model_matrices(). The user-supplied time_var is preserved for formulas as an internal time_covariate column, while an internal numeric time index is used for ordering longitudinal margins and adjacent copula pairs. Numeric and integer time inputs are used directly, numeric-like character time inputs are converted with as.numeric(), and categorical time should be supplied as a factor. Ordered factors are treated with treatment contrasts rather than polynomial contrasts so factor-level effects are interpretable as level comparisons.

The fitting data are converted to a plain data frame. Structurally missing subject-time combinations are expanded to explicit rows with missing responses. Response NA values may be present unless an entire margin or an adjacent copula pair has no complete responses, in which case fitting stops. Submitted predictor values must be observed and finite. Response NA values represent missing outcomes, but response NaN, Inf, and -Inf values stop before fitting. The package does not impute missing values; structurally inserted rows are represented in model.frame(fit, type = "expanded"), while observed response rows are available from model.frame(fit, type = "observed").