Package {antedep}

Convert Bell parameter to Bell mean

Description

Convert Bell parameter to Bell mean

Usage

.bell_mean_from_theta(theta)

Convert Bell mean to Bell parameter

Description

Convert Bell mean to Bell parameter

Usage

.bell_theta_from_mean(mu)

Bell parameter for a given mean, to machine precision

Description

Bell parameter for a given mean, to machine precision

Usage

.bell_theta_newton(mu, start = NULL)

Blend old/new CAT parameters and renormalize probabilities

Description

Blend old/new CAT parameters and renormalize probabilities

Usage

.blend_cat_params_em(
  old_marginal,
  old_transition,
  new_marginal,
  new_transition,
  step,
  order,
  n_time,
  n_blocks,
  homogeneous
)

Arguments

Value

List with blended marginal and transition.

Build AD covariance matrix from parameters

Description

Build AD covariance matrix from parameters

Usage

.build_gau_covariance(order, phi, sigma, n_time)

Arguments

Value

Covariance matrix (n_time x n_time)

Convert probability row to unconstrained logits

Description

Convert probability row to unconstrained logits

Usage

.cat_prob_to_theta(prob, c, eps = 1e-08)

Arguments

Value

Numeric vector of length c-1

Convert unconstrained logits to probability row

Description

Convert unconstrained logits to probability row

Usage

.cat_theta_to_prob(theta_row, c)

Arguments

Value

Probability vector of length c

Upper-tail probability of a chi-bar-square mixture

Description

The j = 0 component is a point mass at zero, so its upper-tail probability P(chi^2_0 >= stat) is 1 when stat <= 0 and 0 when stat > 0.

Usage

.chi_bar_pvalue(stat, r, weights = NULL)

Compute CIs for marginal parameters

Description

Compute CIs for marginal parameters

Usage

.compute_marginal_ci(marginal, cell_counts, p, c, n_time, n_subjects, z, level)

Arguments

Value

Data frame with CIs

Compute CIs for transition parameters

Description

Compute CIs for transition parameters

Usage

.compute_transition_ci(transition, cell_counts, p, c, n_time, z, level)

Arguments

Value

List of data frames with CIs for each time point

Construct AD covariance matrix from fitted model

Description

Construct AD covariance matrix from fitted model

Usage

.construct_gau_covariance(fit, n_time)

Arguments

Value

Covariance matrix.

Count cells for contingency table

Description

Counts occurrences of each combination of categories at specified time indices.

Usage

.count_cells_cat(y, time_indices, n_categories, subject_mask = NULL)

Arguments

Value

Array of counts with dimensions c x c x ... (length = length(time_indices))

Count cells efficiently using table()

Description

Alternative implementation using R's table function.

Usage

.count_cells_table_cat(y, time_indices, n_categories, subject_mask = NULL)

Arguments

Value

Array of counts

Count free parameters for AD(p) categorical model

Description

Count free parameters for AD(p) categorical model

Usage

.count_params_cat(
  order,
  n_categories,
  n_time,
  n_blocks = 1,
  homogeneous = TRUE
)

Arguments

Value

Number of free parameters

Count parameters for AD model

Description

Count parameters for AD model

Usage

.count_params_gau(order, n_time, n_blocks = 1L, estimate_mu = FALSE)

Arguments

Value

Number of free parameters.

Count parameters for homogeneity test

Description

Count parameters for homogeneity test

Usage

.count_params_homogeneity(
  fit,
  test,
  n_blocks,
  order,
  n_time,
  thinning,
  innovation,
  is_null
)

Arguments

Value

Number of free parameters

Free-parameter count of an unconstrained INAD model

Description

Matches .count_params_inad_fit(): order 2 has 2 * n_time - 3 thinning parameters (\alpha_{2,1}, then two per time point), and negative binomial innovations carry one size parameter per time point, as fitted by fit_inad().

Usage

.count_params_inad(order, n_time, n_blocks, thinning, innovation)

Convert counts to probabilities (with safe division)

Description

Convert counts to probabilities (with safe division)

Usage

.counts_to_probs_cat(counts, margin = NULL)

Arguments

Value

Array of probabilities

Compute transition probabilities from counts (simpler version)

Description

Given counts array with last dimension being Y_k and earlier dimensions being conditioning variables, compute P(Y_k | conditioning).

Usage

.counts_to_transition_probs(counts)

Arguments

Value

Array of transition probabilities (same dimensions)

Create default marginal parameters (uniform)

Description

Create default marginal parameters (uniform)

Usage

.default_marginal_cat(p, c, n_time, n_blocks, homogeneous)

Arguments

Value

List of marginal parameters

Create default transition parameters (uniform)

Description

Create default transition parameters (uniform)

Usage

.default_transition_cat(p, c, n_time, n_blocks, homogeneous)

Arguments

Value

List of transition parameters

E-step: Compute expected sufficient statistics

Description

Computes conditional expectations for missing values and second moments.

Usage

.em_e_step_gau(y, order, mu, phi, sigma, blocks, tau)

EM E-step for INAD models

Description

For every subject and time, computes the exact conditional distribution of the latent thinned counts given the observed data, and accumulates (i) the conditional means of the thinned counts, (ii) the conditional pmf of the innovation, aggregated by block and time, and (iii) the observed-data log-likelihood, which is the sum of the log normalizing constants.

Usage

.em_inad_estep(y, order, thinning, innovation, par, blocks, n_blocks)

Innovation pmf parameter by block and time

Description

Returns a n_blocks by n_time matrix of the parameter used by .innov_vec(): the mean for Poisson and negative binomial innovations, and the Bell parameter for Bell innovations.

Usage

.em_inad_lam(theta, tau, n_blocks, innovation)

EM M-step for INAD models

Description

Thinning coefficients use the closed-form update shared by all three operators. Innovation parameters maximize the expected complete-data innovation log-likelihood, evaluated from the conditional innovation pmf.

Usage

.em_inad_mstep(
  y,
  order,
  thinning,
  innovation,
  es,
  par,
  blocks,
  n_blocks,
  nb_inno_size_ub = 50,
  n_cm = 3
)

M-step: Update parameters from sufficient statistics

Description

M-step: Update parameters from sufficient statistics

Usage

.em_m_step_gau(suff_stats, order, blocks, estimate_mu, n_subjects, n_time)

Convert EM expected counts to fit_cat-style cell_counts

Description

Convert EM expected counts to fit_cat-style cell_counts

Usage

.expected_counts_to_cell_counts_cat_em(
  counts,
  order,
  n_time,
  n_blocks,
  homogeneous
)

Arguments

Value

A list shaped like fit_cat(...) cell_counts.

Extract conditional probabilities from transition array

Description

Extract conditional probabilities from transition array

Usage

.extract_conditional_probs(trans_array, history)

Arguments

Value

Vector of probabilities for current time point

Fit CAT model for a single population

Description

Fit CAT model for a single population

Usage

.fit_cat_single_pop(y, p, c, n_time, subject_mask = NULL)

Arguments

Value

List with marginal, transition, cell_counts, log_l

Fit CAT model with missing data via observed-data likelihood optimization

Description

Fit CAT model with missing data via observed-data likelihood optimization

Usage

.fit_cat_single_pop_marginalize(y, p, c, n_time)

Arguments

Value

List with marginal, transition, cell_counts, log_l, convergence

Fit stationary categorical AD model

Description

Internal function to fit a model where both marginal and transition probabilities are constant over time.

Usage

.fit_cat_stationary(y, p, c, blocks = NULL, homogeneous = TRUE)

Arguments

Value

A cat_fit-like object with stationary constraints

Fit stationary model for single population

Description

Fit stationary model for single population

Usage

.fit_cat_stationary_single(y, p, c, n_time, subject_mask = NULL)

Arguments

Value

List with marginal, transition, log_l, n_params

Fit time-invariant categorical AD model

Description

Internal function to fit a model where transition probabilities are constant across time (for k > p).

Usage

.fit_cat_timeinvariant(y, p, c, blocks = NULL, homogeneous = TRUE)

Arguments

Value

A cat_fit-like object with time-invariant transitions

Fit time-invariant model for single population

Description

Fit time-invariant model for single population

Usage

.fit_cat_timeinvariant_single(y, p, c, n_time, subject_mask = NULL)

Arguments

Value

List with marginal, transition, log_l, n_params

Fit an INAD model under time-invariance constraints

Description

Maximizes the observed-data log-likelihood jointly over all parameters left free by the constraint. With complete data it first takes damped Newton steps with the exact score and information from pooled starting values (settings$path = "newton"); if the first-order conditions do not hold at the result, it continues with blockwise BOBYQA updates and a joint BOBYQA refinement, repeated until the log-likelihood stabilizes (settings$path = "search").

Usage

.fit_constrained(
  y,
  order,
  thinning,
  innovation,
  blocks,
  cinfo,
  init_from = NULL,
  na_action = "fail",
  nb_inno_size_ub = 50,
  max_rounds = 25,
  tol = 1e-08,
  verbose = FALSE,
  ...
)

Value

An inad_fit object whose settings$constraint_df records the number of equality constraints, so that n_params, aic, and bic count only the free parameters.

Constrained fit by Newton steps from pooled starting values

Description

Replaces each constrained group of the unconstrained estimates by a pooled value (thinning coefficients weighted by their predecessor totals, with the innovation means shifted so that the fitted means are unchanged), then takes Newton steps with the exact score and information.

Usage

.fit_constrained_newton(
  y,
  order,
  thinning,
  innovation,
  blocks_id,
  B,
  cinfo,
  a1,
  a2,
  theta,
  kappa,
  tau,
  nb_inno_size_ub,
  ub_a1,
  ub_a2
)

Value

NULL if the starting values are infeasible, otherwise a list with the fit and kkt_ok.

EM algorithm for AD with missing data

Description

Implements the EM algorithm for fitting AD(p) models with missing data under MAR (Missing At Random) assumption.

Usage

.fit_gau_em(
  y,
  order,
  blocks = NULL,
  estimate_mu = TRUE,
  max_iter = 100,
  tol = 1e-06,
  verbose = FALSE,
  ...
)

Arguments

Value

List with fitted model (same structure as fit_gau, plus EM info)

Fit fully heterogeneous INAD model

Description

Internal function that fits separate INAD models for each group and combines the log-likelihoods.

Usage

.fit_inad_heterogeneous(
  y,
  blocks,
  n_blocks,
  group_indices,
  order,
  thinning,
  innovation,
  ...
)

Arguments

Value

A list mimicking fit_inad output with combined log-likelihood

Get all category combinations of given length

Description

Get all category combinations of given length

Usage

.get_combinations_cat(n_categories, length)

Arguments

Value

Matrix where each row is a combination (c^length rows, length columns)

Extract history for a subject at time k

Description

Extract history for a subject at time k

Usage

.get_history(y_row, k, order)

Arguments

Value

Integer vector of length min(k-1, order) with category history

Compute group-specific means for two-sample case

Description

Compute group-specific means for two-sample case

Usage

.group_means_gau(y, blocks)

Arguments

Value

List with mu1 and mu2 vectors.

Convert history to array index

Description

Convert history to array index

Usage

.history_to_index(history, n_categories)

Arguments

Value

Integer index (1-based) into flattened array

Thinning coefficients that are not identifiable

Description

Thinning coefficients that are not identifiable

Usage

.inad_alpha_unidentified(y, lag)

Value

Logical vector of length ncol(y): TRUE where the counts at lag lag before the occasion are all zero (or all missing).

Upper limits for thinning coefficients

Description

Binomial thinning requires \alpha \le 1; Poisson and negative binomial thinning require only \alpha \ge 0. For the latter, Fisher's identity gives \sum_s E(Z_{si,j} \mid y_s) = \alpha_{i,i-j} \sum_s y_{s,i-j} at any stationary point with \alpha_{i,i-j} > 0, and Z_{si,j} \le y_{si}, so every stationary point satisfies \alpha_{i,i-j} \le \sum_s y_{si} / \sum_s y_{s,i-j}. Twice that ratio is used as a box limit for the optimizers: it is never active at a maximum, and it keeps the box finite for derivative-free methods. With common = TRUE the limit is for one coefficient shared by all occasions with the given lag, from the sums over those occasions. A coefficient whose predecessor counts are all zero is not identifiable; its limit is 1e-8, and fitting functions set it to 0.

Usage

.inad_alpha_upper(y, lag, thinning, common = FALSE)

Arguments

Value

A vector of length ncol(y) (0 at occasions without the lag), or a scalar when common = TRUE.

Build transition objects for truncated-state missing-data likelihood

Description

Build transition objects for truncated-state missing-data likelihood

Usage

.inad_build_transitions(
  order,
  N,
  K,
  alpha1,
  alpha2,
  thinning,
  innovation,
  lambda,
  nb_inno_size
)

Wald intervals from an observed information matrix

Description

Wald intervals from an observed information matrix

Usage

.inad_ci_from_info(tab, info, level = 0.95, log_scale = TRUE, ridge = 0)

Truncated discrete convolution

Description

Truncated discrete convolution

Usage

.inad_conv_trunc(p, q, K)

Line search along one free parameter, starting from a bound

Description

Maximizes the observed-data log-likelihood over the parameter in row k of the free-parameter table, all other parameters fixed, between its current value and lim. The search starts with a step of size t0 towards lim, shortens it until the log-likelihood increases, lengthens it while the log-likelihood keeps increasing, and finishes with stats::optimize on the bracket found.

Usage

.inad_coord_search(loglik, fit, fp, k, lim, t0, N, ll0 = loglik(fit))

Value

A list with the best value par, the log-likelihood increase gain (0 if no increase was found), and the corresponding fit.

Effective innovation mean under block effects

Description

Effective innovation mean under block effects

Usage

.inad_effective_innovation_mean(theta, tau, innovation)

Effective innovation distribution parameter under block effects

Description

For Poisson and negative binomial innovations, block effects are additive on the mean. For Bell innovations, block effects are additive on the innovation mean theta * exp(theta), then mapped back to the Bell parameter.

Usage

.inad_effective_innovation_param(theta, tau, innovation)

Row of the free-parameter table for a local (time or block) parameter

Description

Row of the free-parameter table for a local (time or block) parameter

Usage

.inad_free_index(tab, group, index)

Free parameters of a fitted INAD model

Description

Free parameters of a fitted INAD model

Usage

.inad_free_params(fit, N, B)

Value

A list with tab, a data frame with one row per free parameter (columns component, group, param, est, index), and the unpacked parameter vectors.

First-order conditions for a fit-like list (see kkt_inad)

Description

First-order conditions for a fit-like list (see kkt_inad)

Usage

.inad_kkt_core(y, fit, blocks, tol = 1e-06)

Observed-data log-likelihood of a fit-like list, as a function of the fit

Description

Observed-data log-likelihood of a fit-like list, as a function of the fit

Usage

.inad_loglik_fun(y, fit, blocks = NULL)

Build truncated thinning pmf for one previous count

Description

Build truncated thinning pmf for one previous count

Usage

.inad_make_thin_pmf(prev, alpha, thinning, K)

Newton polish of an INAD fit using the exact score and information

Description

Takes Newton-Raphson steps on the free parameters, using the Louis score and observed information, with step halving so that the log-likelihood never decreases and parameters stay in their parameter space. Used to finish derivative-free fits to full precision.

Usage

.inad_newton_polish(
  y,
  fit,
  blocks = NULL,
  nb_inno_size_ub = 50,
  max_iter = 25,
  tol = 1e-10
)

Value

A list with the polished fit (parameters only) and its log_l.

Numerical observed-data Hessian with fixed absolute steps

Description

Positive parameters are differentiated on the log scale and block effects on their natural scale. Steps are absolute (step, refined by Richardson extrapolation), because steps proportional to the parameter value become too small for log-scale parameters near zero.

Usage

.inad_numeric_info(y, fit, blocks = NULL, step = 0.001)

Value

The observed information for the free parameters on their natural scale.

Retrieve/calculate order-2 transition distribution for one previous state pair

Description

Retrieve/calculate order-2 transition distribution for one previous state pair

Usage

.inad_order2_dist(prev2, prev1, info, K)

Finish a fit by Newton steps with the exact score and information

Description

Accepts the polished parameters only if they do not decrease the log-likelihood, and refreshes the derived quantities of the fit.

Usage

.inad_polish_fit(fit, y, nb_inno_size_ub = 50)

Write free-parameter values back into a fit-like list

Description

Write free-parameter values back into a fit-like list

Usage

.inad_set_free(fit, fp, values, N)

Heuristic state-space bound for INAD missing-data recursion

Description

Heuristic state-space bound for INAD missing-data recursion

Usage

.inad_state_max(y, order, alpha1, alpha2, theta, tau, innovation = "pois")

Subject-level observed-data likelihood for INAD(0)

Description

Subject-level observed-data likelihood for INAD(0)

Usage

.inad_subject_ll_order0(y_s, innov_list, K)

Subject-level observed-data likelihood for INAD(1)

Description

Subject-level observed-data likelihood for INAD(1)

Usage

.inad_subject_ll_order1(y_s, init_pmf, trans1, K)

Subject-level observed-data likelihood for INAD(2)

Description

Subject-level observed-data likelihood for INAD(2)

Usage

.inad_subject_ll_order2(y_s, init_pmf, trans_t2, order2_info, K)

Initialize parameters for AD EM algorithm

Description

Initialize parameters for AD EM algorithm

Usage

.initialize_gau_em(y, order, blocks, estimate_mu)

Initialize from marginal statistics (ignoring dependence)

Description

Initialize from marginal statistics (ignoring dependence)

Usage

.initialize_gau_marginal(y, order, blocks)

Innovation probability vector

Description

Computes the probability mass for innovation counts.

Usage

.innov_vec(u_vals, lam, innovation, sz)

Arguments

Value

Numeric vector of probabilities.

Compute innovation variance estimates under AD(p)

Description

The MLE of the innovation variance at time i is RSS_i(p) / N.

Usage

.innovation_var_gau(y, p, mu = NULL)

Arguments

Value

Numeric vector of innovation variance estimates.

Compute log-likelihood from cell counts (faster for large datasets)

Description

Instead of looping over subjects, this version uses aggregated cell counts. Equivalent to logL_cat but more efficient.

Usage

.logL_from_counts_cat(y, order, marginal, transition, n_categories)

Arguments

Value

Scalar log-likelihood

Compute observed-data log-likelihood for AD with missing values

Description

Uses multivariate normal marginalization to compute the likelihood of observed values, marginalizing over missing values.

Usage

.logL_gau_missing(y, order, mu, phi, sigma, blocks = NULL, tau = 0)

Arguments

Value

Observed-data log-likelihood (scalar)

Observed-data INAD likelihood with missing values under MAR

Description

Uses forward recursion over a truncated count state space.

Usage

.logL_inad_missing(
  y,
  order,
  thinning,
  innovation,
  alpha,
  theta,
  nb_inno_size = NULL,
  blocks = NULL,
  tau = 0
)

Compute log-likelihood contribution from one subject

Description

Compute log-likelihood contribution from one subject

Usage

.logL_subject_cat(y_s, p, c, marginal, transition)

Arguments

Value

Scalar log-likelihood contribution

Observed-data log-likelihood for order-1 CAT model

Description

Observed-data log-likelihood for order-1 CAT model

Usage

.logL_subject_cat_marginalize_p1(y_s, c, marginal, transition)

Arguments

Value

Scalar log-likelihood contribution

Observed-data log-likelihood for order-2 CAT model

Description

Observed-data log-likelihood for order-2 CAT model

Usage

.logL_subject_cat_marginalize_p2(y_s, c, marginal, transition)

Arguments

Value

Scalar log-likelihood contribution

Compute observed-data log-likelihood contribution from one subject

Description

Compute observed-data log-likelihood contribution from one subject

Usage

.logL_subject_cat_observed(y_s, p, c, marginal, transition, na_action)

Arguments

Value

Scalar log-likelihood contribution

Log-sum-exp for numerical stability

Description

Log-sum-exp for numerical stability

Usage

.log_sum_exp(x)

Compute log-likelihood contribution from counts and probabilities

Description

Computes sum of N * log(pi) over all cells.

Usage

.loglik_contribution(counts, probs)

Arguments

Value

Scalar log-likelihood contribution

Stable log-sum-exp

Description

Stable log-sum-exp

Usage

.logsumexp_cat(x)

Arguments

Value

Scalar log(sum(exp(x)))

Observed score and information for all free parameters, by Louis' identity

Description

Observed score and information for all free parameters, by Louis' identity

Usage

.louis_full_info(y, fit, blocks = NULL)

One-sided observed-data derivatives for a thinning coefficient on a bound

Description

At alpha = 0 the thinned count Z is identically 0, and at binomial alpha = 1 it equals yprev; the support of the complete-data distribution changes there, so the conditional mean of the complete-data score is 0/0 and Fisher's identity cannot be evaluated at the bound. The observed-data likelihood is f(y) = \sum_z p_\alpha(z) h(y - z), where h is the pmf of the other components of y. With n = yprev and H_k = h(y - k), the derivatives of p_\alpha at \alpha = 0 are -n, n, 0, \ldots for z = 0, 1, 2, \ldots (all three thinning operators), so the right derivative is n (H_1 / H_0 - 1). The second derivatives are c, -2c, c for z = 0, 1, 2, with c = n(n - 1), n^2, and n(n + 1) for binomial, Poisson, and negative binomial thinning. Binomial thinning at \alpha = 1 is the mirror image in 1 - \alpha, with z = n, n - 1, n - 2, giving the left derivative -n (H_{n-1} / H_n - 1).

Usage

.louis_thin_boundary(z, lp_other, yprev, alpha, thinning)

Arguments

Value

A named vector: the one-sided derivative s of log f(y) with respect to alpha, and b, minus the one-sided second derivative.

Complete-data score and negative Hessian of one thinned count

Description

Valid for 0 < alpha (and alpha < 1 for binomial thinning). At alpha = 0, or binomial alpha = 1, the thinned count is degenerate and the conditional mean of this score is 0/0; see .louis_thin_boundary.

Usage

.louis_thin_terms(z, yprev, alpha, thinning)

Value

A list with the score s and the negative second derivative b of log P(alpha odot yprev = z) with respect to alpha.

Louis score and information for the parameters entering time i

Description

Local parameters are a1 (i >= 2), a2 (order 2, i >= 3), theta, kappa (negative binomial innovations), and the block effects tau2, ..., tauB. Returns the observed score (the conditional mean of the complete-data score) and the observed information E(B \mid y) - \mathrm{Cov}(S \mid y), both summed over subjects.

Usage

.louis_time_terms(
  y,
  i,
  ord,
  thinning,
  innovation,
  a1,
  a2,
  theta_i,
  kappa_i,
  tau,
  blocks
)

Compute MLE of mean vector under AD(p)

Description

Under AD(p), the MLE of mu is the sample mean vector when there are no covariates.

Usage

.mle_mean_gau(y)

Arguments

Value

Numeric vector of mean estimates.

Pack CAT parameters into unconstrained vector

Description

Pack CAT parameters into unconstrained vector

Usage

.pack_cat_params(marginal, transition, p, c, n_time)

Arguments

Value

Numeric parameter vector

Parse a time-invariance constraint specification

Description

Parse a time-invariance constraint specification

Usage

.parse_constraint(constrain, order, n_time, innovation = "pois")

Compute intervenor-adjusted sample partial correlation

Description

Computes the sample partial correlation between Y_i and Y_j, adjusted for the intervenors Y_(j+1), ..., Y_(i-1), from the residuals of the fitted mean structure.

Usage

.partial_corr_gau(y, i, j, mu = NULL)

Arguments

Details

The intervenor-adjusted partial correlation r_(i,j|(j+1:i-1)) is computed from residuals of regressing Y_i and Y_j on the intervenors.

Under AD(p), the partial correlation r_(i,i-k|(i-k+1:i-1)) should be approximately zero for k > p.

Value

Sample partial correlation (scalar).

Compute matrix of intervenor-adjusted partial correlations

Description

Compute matrix of intervenor-adjusted partial correlations

Usage

.partial_corr_matrix_gau(y, mu = NULL)

Arguments

Value

Matrix where entry (i,j) for i > j is r_(i,j|(j+1:i-1)).

Compute partial residuals for PRISM plot

Description

Compute partial residuals for PRISM plot

Usage

.prism_residuals(y, i, j)

Arguments

Value

List with resid_i and resid_j (partial residuals).

Modified psi function for test statistic correction

Description

Computes the psi function used in the modified likelihood ratio test (Kenward, 1987) for better small-sample approximation.

Usage

.psi_kenward(x, y)

Arguments

Details

psi(x, y) is defined as:

• 0 if x = 0

• (2y - 1) / (2y(y - 1)) if x = 1

• 2 * sum from l=1 to x/2 of (y + 2l - 2)^(-1) if x > 0 and even

• psi(1, y) + psi(x - 1, y + 1) if x > 1 and odd

Value

Value of psi(x, y).

Compute residual sum of squares from AD regression

Description

For time point i, regress Y_i on its p_i predecessors (Y_(i-1), ..., Y_(i-p_i)) and return the residual sum of squares.

Usage

.rss_gau(y, i, p, mu = NULL, include_intercept = TRUE)

Arguments

Details

When mu is provided, the regression is: (Y_i - mu_i) ~ (Y_(i-1) - mu_(i-1)) + ... + (Y_(i-p) - mu_(i-p)) with no intercept.

When mu is NULL and include_intercept is TRUE: Y_i ~ 1 + Y_(i-1) + ... + Y_(i-p)

For i = 1 or p = 0, the RSS is simply the sum of squared deviations from the mean (or from mu_1 if provided).

Value

Residual sum of squares (scalar).

Compute RSS for two-sample case (pooled and separate)

Description

Compute RSS for two-sample case (pooled and separate)

Usage

.rss_two_sample_gau(y, blocks, i, p, pooled = TRUE)

Arguments

Value

RSS value.

Compute RSS vector for all time points under AD(p)

Description

Compute RSS vector for all time points under AD(p)

Usage

.rss_vector_gau(y, p, mu = NULL)

Arguments

Value

Numeric vector of length n_time with RSS for each time point.

Safe log function (0 * log(0) = 0)

Description

Safe log function (0 * log(0) = 0)

Usage

.safe_log(x)

Arguments

Value

Log of x, with 0 for x = 0

Safeguard M-step update via step-halving

Description

Safeguard M-step update via step-halving

Usage

.safeguard_update_cat_em(
  y,
  order,
  old_marginal,
  old_transition,
  new_marginal,
  new_transition,
  blocks_id,
  homogeneous,
  c,
  logL_old,
  n_time,
  n_blocks
)

Arguments

Value

List with safeguarded parameters, acceptance flag, and step size.

Simulate one subject's trajectory

Description

Simulate one subject's trajectory

Usage

.simulate_subject_cat(n_time, p, c, marginal, transition)

Arguments

Value

Integer vector of length n_time

Is a thinning coefficient on a bound where the thinned count is degenerate?

Description

Is a thinning coefficient on a bound where the thinned count is degenerate?

Usage

.thin_at_bound(alpha, thinning)

Thinning probability vector

Description

Computes the probability mass for thinned counts.

Usage

.thin_vec(k_vals, yprev, a, thinning)

Arguments

Value

Numeric vector of probabilities.

Build uniform CAT parameter values

Description

Build uniform CAT parameter values

Usage

.uniform_cat_params(p, c, n_time)

Arguments

Value

List with marginal and transition

Unpack alpha parameters

Description

Converts alpha from various input formats to a standardized list.

Usage

.unpack_alpha(order, alpha, N)

Arguments

Value

List with elements a1 and a2.

Unpack unconstrained CAT parameter vector

Description

Unpack unconstrained CAT parameter vector

Usage

.unpack_cat_params(theta, p, c, n_time)

Arguments

Value

List with marginal and transition parameter lists

Validate blocks parameter

Description

Validate blocks parameter

Usage

.validate_blocks_cat(blocks, n_subjects)

Arguments

Value

Validated blocks vector (integer, 1-indexed)

Validate transition probability parameters

Description

Checks that transition probabilities are valid (non-negative, sum to 1).

Usage

.validate_params_cat(marginal, transition, p, c, n_time, tol = 1e-08)

Arguments

Value

TRUE if valid, otherwise throws error

Validate categorical data matrix

Description

Validate categorical data matrix

Usage

.validate_y_cat(y, n_categories = NULL, allow_na = FALSE)

Arguments

Value

List with validated y and n_categories

The Bell distribution

Description

Density, distribution function, quantile function and random generation for the Bell distribution with parameter theta.

Usage

dbell(x, theta, log = FALSE)

pbell(x, theta)

rbell(n, theta, max_z = 100L)

qbell(p, theta, max_z = 100L)

Arguments

Details

Let B_x denote the xth Bell number. The Bell distribution has probability mass function

P(X = x) = \theta^x \exp(-\exp(\theta) + 1) \frac{B_x}{x!},

for nonnegative integers x and \theta \ge 0.

For \theta > 0, the Bell mean is E[X] = \theta e^\theta. At \theta = 0, the distribution is degenerate at 0.

The functions follow the standard naming used in base R: dbell for the density, pbell for the distribution function, qbell for the quantile function and rbell for random generation.

Value

For dbell, a numeric vector of probabilities. For pbell, a numeric vector of cumulative probabilities. For qbell, an integer vector of quantiles. For rbell, an integer vector of random values.

Examples

dbell(0:5, theta = 1)
pbell(0:5, theta = 1)
qbell(c(0.25, 0.5, 0.9), theta = 1)
set.seed(1)
rbell(10, theta = 1)

Bell numbers

Description

Compute Bell numbers up to order z using a truncated series representation.

Usage

Bz(z)

Arguments

Details

This function is for internal use. The results are stored in .Bell_num and used by the Bell distribution functions.

Value

numeric vector of length z + 1 containing B_0, ..., B_z.

Likelihood-Ratio Comparisons for Nested AD Fits

Description

Compares one or more nested antedependence model fits in the order supplied.

Usage

## S3 method for class 'gau_fit'
anova(object, ..., test = "Chisq")

## S3 method for class 'cat_fit'
anova(object, ..., test = "Chisq")

## S3 method for class 'inad_fit'
anova(object, ..., test = "Chisq")

Arguments

Details

For INAD fits, a comparison between models of different orders is a boundary problem (the added thinning coefficients are zero under the smaller model), so no chi-square p-value is reported for it; see test_order_inad.

Value

An anova data frame.

Convert Antedependence Simulation Output to a Time-Series Object

Description

Convert Antedependence Simulation Output to a Time-Series Object

Usage

## S3 method for class 'gau_sim'
as.ts(x, ...)

## S3 method for class 'cat_sim'
as.ts(x, ...)

## S3 method for class 'inad_sim'
as.ts(x, ...)

Arguments

Value

A ts or mts object.

Bayesian and Akaike Information Criteria for Antedependence Fits

Description

Compute BIC or AIC for a fitted Gaussian, categorical, or INAD antedependence model. Since logLik() methods are registered for all three fit classes, stats::AIC() and stats::BIC() generics also work directly and are preferred in new code. The family-specific helpers are retained as thin compatibility wrappers for existing code.

Usage

bic_cat(fit, n_subjects = NULL)

aic_cat(fit)

bic_gau(fit, n_subjects = NULL)

aic_gau(fit)

bic_inad(fit, n_subjects = NULL)

aic_inad(fit)

Arguments

Details

The BIC is computed as:

BIC = -2 \times \ell + k \times \log(N)

where \ell is the log-likelihood, k is the number of free parameters, and N is the number of subjects.

The AIC is computed as:

AIC = -2 \times \ell + 2k

Value

A numeric scalar BIC or AIC value.

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 30, n_time = 5, order = 1, phi = 0.3)
fit <- fit_gau(y, order = 1)
bic_gau(fit, n_subjects = nrow(y))
aic_gau(fit)

BIC-Based Order Selection for Categorical AD Models

Description

Fits AD models of increasing orders and selects the best by BIC.

Usage

bic_order_cat(
  y,
  max_order = 2,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  criterion = "bic"
)

Arguments

Value

A list containing:

Examples

y <- simulate_cat(100, 5, order = 1, n_categories = 2)
result <- bic_order_cat(y, max_order = 2)
print(result)
print(result$best_order)

BIC-Based Order Selection for Gaussian AD Models

Description

Fits AD models of increasing orders and selects the best by BIC.

Usage

bic_order_gau(y, max_order = 2L, ...)

Arguments

Value

A list with class gau_bic_order containing:

fits

List of fitted models

bic

BIC values for each order

best_order

Order with lowest BIC

table

Summary table

See Also

bic_order_cat, bic_order_inad, fit_gau

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 80, n_time = 6, order = 1, phi = 0.4)
ord <- bic_order_gau(y, max_order = 2)
ord$best_order
print(ord)

BIC-Based Order Selection for INAD Models

Description

Fits INAD models across candidate orders and reports BIC-based selection.

Usage

bic_order_inad(
  y,
  max_order = 2,
  thinning = "binom",
  innovation = "pois",
  blocks = NULL,
  ...
)

Arguments

Value

A list with class "bic_order_inad" containing:

fits

List of fitted INAD models by candidate order

bic

Named numeric vector of BIC values by order

best_order

Order with minimum BIC

table

Data frame with order, logLik, n_params, and BIC

settings

Input and derived settings used for selection

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 30,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
ord <- bic_order_inad(y, max_order = 2, thinning = "binom", innovation = "pois", max_iter = 10)
ord$best_order

Morphine Bolus Analgesia Counts

Description

Morphine bolus self administration counts for two treatment groups recorded at 12 four hour time points. The data are stored in matrix form to facilitate use with antedependence models.

Usage

bolus_inad

Format

A list with four components:

y

integer matrix of dimension N by n_time containing all subjects and time points

y_2mg

integer matrix with rows corresponding to the 2 mg treatment group

y_1mg

integer matrix with rows corresponding to the 1 mg treatment group

blocks

integer vector of length N giving the block or treatment group indicator, 1 for 2 mg and 2 for 1 mg

Source

Dataset bolus from the cold package, converted to matrix form and grouped by treatment.

Cattle Growth Data (Treatments A and B)

Description

Longitudinal cattle growth measurements for two treatment groups from Zimmerman and Nunez-Anton antedependence book companion data. This dataset is continuous-response data suitable for Gaussian AD modeling.

Usage

cattle_growth

Format

A list with five components:

y

numeric matrix of dimension N by n_time containing all subjects

y_A

numeric matrix for Treatment A subjects

y_B

numeric matrix for Treatment B subjects

blocks

integer vector of length N (1 = Treatment A, 2 = Treatment B)

time

integer vector of measurement occasions

Source

https://homepage.divms.uiowa.edu/~dzimmer/Data-for-AD/cattle_growth_data_Treatment%20A.txt and https://homepage.divms.uiowa.edu/~dzimmer/Data-for-AD/cattle_growth_data_Treatment%20B.txt

Confidence Intervals for Fitted Categorical AD Models

Description

Computes Wald-based confidence intervals for the transition probability parameters of a fitted categorical antedependence model.

Usage

ci_cat(fit, y = NULL, level = 0.95, parameters = "all")

Arguments

Details

Confidence intervals are computed using the Wald method based on the asymptotic normality of maximum likelihood estimators.

For a probability estimate \hat{\pi} based on count N, the standard error is:

SE(\hat{\pi}) = \sqrt{\frac{\hat{\pi}(1-\hat{\pi})}{N}}

For conditional probabilities \hat{\pi}_{j|i} based on conditioning count N_i, the standard error is:

SE(\hat{\pi}_{j|i}) = \sqrt{\frac{\hat{\pi}_{j|i}(1-\hat{\pi}_{j|i})}{N_i}}

The confidence interval is then:

\hat{\pi} \pm z_{\alpha/2} \times SE(\hat{\pi})

Note: CIs are truncated to the interval from 0 to 1 when they exceed these bounds.

Missing-data fits with na_action = "marginalize" are not currently supported because observed cell counts are not stored for that path.

Value

A list of class "cat_ci" containing:

References

Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley.

See Also

fit_cat

Examples

# Fit a model
set.seed(123)
y <- simulate_cat(200, 5, order = 1, n_categories = 2)
fit <- fit_cat(y, order = 1)

# Compute confidence intervals
ci <- ci_cat(fit)
print(ci)

# Just marginal CIs
ci_marg <- ci_cat(fit, parameters = "marginal")

Confidence Intervals for Fitted Gaussian AD Models

Description

Computes approximate Wald confidence intervals for selected parameters from a fitted Gaussian AD model.

Usage

ci_gau(fit, level = 0.95, parameters = "all")

Arguments

Details

This helper currently supports complete-data Gaussian AD fits.

Standard errors are based on large-sample approximations:

• SE(\hat{\mu}_t) \approx \hat{\sigma}_t / \sqrt{n}

• SE(\hat{\sigma}_t) \approx \hat{\sigma}_t / \sqrt{2n}

• SE(\hat{\phi}) \approx \sqrt{(1-\hat{\phi}^2)/n} for free \phi entries

Value

An object of class gau_ci, a list with elements settings, level, mu, phi, and sigma. Each non-NULL element is a data frame with columns param, est, se, lower, upper, and level.

See Also

fit_gau, ci_cat, ci_inad

Examples

y <- simulate_gau(n_subjects = 80, n_time = 6, order = 1, phi = 0.4)
fit <- fit_gau(y, order = 1)
ci <- ci_gau(fit)
ci$mu
ci$phi
ci$sigma

Confidence Intervals for Fitted INAD Models

Description

Computes Wald confidence intervals for the free parameters of a fitted INAD model from the observed information.

Usage

ci_inad(
  y,
  fit,
  blocks = NULL,
  level = 0.95,
  idx_time = NULL,
  ridge = 0,
  profile_maxeval = 2500,
  profile_xtol_rel = 1e-06,
  method = c("louis", "hessian"),
  log_scale = TRUE,
  profile_tau = TRUE,
  hessian_step = 0.001
)

Arguments

Details

With method = "louis" (the default), the observed information is computed exactly by Louis' identity: for every time point the latent thinned counts have a finite conditional distribution given the data, so the conditional moments of the complete-data score and Hessian are finite sums. The time-point contributions are combined over all free parameters jointly, including block effects and parameters shared across time by fit_inad(constrain = ...). With method = "hessian", the observed information is instead the negative Hessian of the observed-data log-likelihood computed numerically (Richardson extrapolation with fixed absolute steps; positive parameters on the log scale); the two methods agree up to numerical error, except that numerical differentiation is unreliable for estimates very close to zero, for which the Louis calculation should be used. Estimates on the boundary of the parameter space (a thinning coefficient at zero or, for binomial thinning, at one, or a negative binomial innovation size at its cap nb_inno_size_ub) are held fixed and receive no Wald interval. The information is computed from exact finite sums (up to floating-point arithmetic); the resulting Wald intervals are large-sample approximations.

By default, intervals for the positive parameters (alpha, theta, nb_inno_size) are formed on the log scale and transformed back; intervals for block effects use their natural scale. For unconstrained fits with block effects, profile-likelihood intervals for tau are also computed (element tau_profile).

Value

An object of class inad_ci, a list with elements settings, level, method, alpha, theta, nb_inno_size, tau (Wald), tau_profile (when computed), and vcov, the estimated covariance matrix of the free parameters. Each interval element is a data frame with columns param, est, se, lower, upper, width, level, and scale.

References

Louis, T.A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society, Series B, 44, 226-233.

Examples

data("bolus_inad", package = "antedep")
y <- bolus_inad$y[, 1:4]  # first four occasions, to keep the example short
blocks <- bolus_inad$blocks
fit <- fit_inad(y, order = 1, thinning = "binom", innovation = "pois", blocks = blocks)
ci <- ci_inad(y, fit, blocks = blocks, profile_tau = FALSE)
ci$alpha
ci$theta
ci$tau

Cochlear Implant Speech Recognition Data

Description

Longitudinal speech recognition outcomes for two groups (A/B), including incomplete records, from Zimmerman and Nunez-Anton antedependence book companion data. This dataset is continuous-response data suitable for Gaussian AD modeling.

Usage

cochlear_implant

Format

A list with six components:

y

numeric matrix of dimension N by n_time containing all subjects

y_A

numeric matrix for Group A subjects

y_B

numeric matrix for Group B subjects

blocks

integer vector of length N (1 = Group A, 2 = Group B)

group

character vector of group labels ("A"/"B")

time

integer vector of measurement occasions

Source

https://homepage.divms.uiowa.edu/~dzimmer/Data-for-AD/speech_recognition_data.txt

Extract Model Coefficients from Antedependence Fits

Description

By default, returns a named numeric vector so standard tools can consume fitted objects. Use type = "list" to inspect the family-specific structured representation.

Usage

## S3 method for class 'gau_fit'
coef(object, type = c("vector", "list"), ...)

## S3 method for class 'cat_fit'
coef(object, type = c("vector", "list"), ...)

## S3 method for class 'inad_fit'
coef(object, type = c("vector", "list"), ...)

Arguments

Value

A named numeric vector or a family-specific named list.

Confidence Intervals for Antedependence Model Fits

Description

Returns a standard two-column matrix with lower and upper columns. The structured helpers ci_gau(), ci_cat(), and ci_inad() remain available for detailed inspection.

Usage

## S3 method for class 'gau_fit'
confint(object, parm = NULL, level = 0.95, ...)

## S3 method for class 'cat_fit'
confint(object, parm = NULL, level = 0.95, y = NULL, ...)

## S3 method for class 'inad_fit'
confint(object, parm = NULL, level = 0.95, y, ...)

Arguments

Value

A matrix with columns lower and upper.

Deviance for Antedependence Model Fits

Description

Deviance for Antedependence Model Fits

Usage

## S3 method for class 'ad_fit'
deviance(object, ...)

## S3 method for class 'gau_fit'
deviance(object, ...)

## S3 method for class 'cat_fit'
deviance(object, ...)

## S3 method for class 'inad_fit'
deviance(object, ...)

Arguments

Value

A scalar deviance value.

EM Algorithm for Categorical AD Model Estimation

Description

Fits categorical antedependence models with missing outcomes using the Expectation-Maximization (EM) algorithm for orders 0 and 1.

Usage

em_cat(
  y,
  order = 1,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  max_iter = 100,
  tol = 1e-06,
  epsilon = 1e-08,
  safeguard = TRUE,
  verbose = FALSE
)

Arguments

Details

For complete data (no missing values), this function defers to fit_cat with closed-form MLEs.

For missing data and orders 0/1, each EM iteration computes expected sufficient statistics with a forward-backward E-step, then updates probabilities by normalized expected counts in the M-step. If safeguard = TRUE, a step-halving line search is applied to the M-step update whenever the observed-data likelihood decreases.

A final E-step is run before returning so that log_l/AIC/BIC and expected cell counts correspond exactly to the returned parameter values.

Value

A cat_fit object with fields matching fit_cat. In EM mode, cell_counts stores expected counts from the final E-step, with settings$cell_counts_type = "expected".

See Also

fit_cat, logL_cat

Examples

set.seed(1)
y <- simulate_cat(n_subjects = 40, n_time = 5, order = 1, n_categories = 3)
y[sample(length(y), 10)] <- NA
fit <- em_cat(y, order = 1, n_categories = 3, max_iter = 20, tol = 1e-5)
fit$settings$na_action

EM Algorithm for Gaussian AD Model Estimation

Description

Convenience wrapper around fit_gau with na_action = "em" to provide a parallel entry point to em_inad.

Usage

em_gau(
  y,
  order = 1,
  blocks = NULL,
  estimate_mu = TRUE,
  max_iter = 100,
  tol = 1e-06,
  verbose = FALSE,
  ...
)

Arguments

Details

This is an alias-style helper for users who prefer explicit em_* entry points across model families.

Value

An gau_fit object as returned by fit_gau.

See Also

fit_gau, em_inad, em_cat, fit_cat

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 35, n_time = 5, order = 1, phi = 0.3)
y[sample(length(y), 8)] <- NA
fit <- em_gau(y, order = 1, max_iter = 20, tol = 1e-5)
fit$settings$na_action

EM Algorithm for INAD Model Estimation

Description

Fits INAD models using the Expectation-Maximization algorithm, treating the thinned counts as missing data. This is an alternative to direct likelihood optimization by fit_inad.

Usage

em_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  blocks = NULL,
  max_iter = 2000,
  tol = 1e-09,
  alpha_init = NULL,
  theta_init = NULL,
  tau_init = NULL,
  nb_inno_size = NULL,
  nb_inno_size_ub = 50,
  safeguard = TRUE,
  newton_polish = FALSE,
  verbose = FALSE
)

Arguments

Details

For Gaussian and CAT EM entry points, see em_gau and em_cat. For CAT specifically, fit_cat() supports na_action = "em" for orders 0/1 and na_action = "marginalize" for order 2 missing-data fits.

The complete data for subject s at time i are the thinned counts Z_{si,j} = \alpha_{i,i-j} \odot Y_{s,i-j} and the innovation w_{si} = y_{si} - \sum_j Z_{si,j}. Given the observed data, the latent vector at time i has finite support and a conditional distribution that depends only on y_{si} and its p predecessors. The E-step computes this conditional distribution exactly, so the expected complete-data log-likelihood is evaluated without approximation for every thinning-innovation combination.

For binomial, Poisson, and negative binomial thinning, the parameter-dependent part of the thinning log-likelihood is linear in the thinned count, and the M-step for the thinning coefficient has the closed form \hat\alpha_{i,i-j} = \sum_s E(Z_{si,j} \mid y_s) / \sum_s y_{s,i-j}. For Poisson and Bell innovations the innovation log-likelihood is also linear in w_{si}. For negative binomial innovations it contains \log\Gamma(w_{si} + \kappa_i), which is not linear in w_{si}; its conditional expectation is computed from the conditional distribution of w_{si}, and the size parameters \kappa_i are updated by maximizing the resulting expected log-likelihood. Without block effects the innovation updates are exact maximizations. With block effects the innovation means, sizes, and block effects are updated by conditional maximization (an ECM algorithm), which also never decreases the observed-data log-likelihood.

Value

A list with class "inad_fit" containing estimated parameters, the maximized log-likelihood, n_params, aic, bic, the number of iterations, and the log-likelihood trace.

References

Dempster, A.P., Laird, N.M. and Rubin, D.B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B, 39, 1-38.

Meng, X.-L. and Rubin, D.B. (1993). Maximum likelihood estimation via the ECM algorithm: a general framework. Biometrika, 80, 267-278.

See Also

em_gau, em_cat, fit_inad, fit_cat

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 50,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.25,
  theta = 2
)
fit <- em_inad(y, order = 1, thinning = "binom", innovation = "pois", max_iter = 20, tol = 1e-6)
fit$log_l

Unified Antedependence Model Fitting Interface

Description

Fits one of the package's antedependence model families through a common entry point. The family-specific functions remain available and are still the clearest choice when using family-specific arguments.

Usage

fit_ad(y, family = c("gaussian", "categorical", "inad"), ...)

Arguments

Value

A fitted model object inheriting from ad_fit.

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 30, n_time = 5, order = 1, phi = 0.3)
fit <- fit_ad(y, family = "gaussian", order = 1)
BIC(fit)

Fit Categorical Antedependence Model by Maximum Likelihood

Description

Computes maximum likelihood estimates for the parameters of an AD(p) model for categorical longitudinal data. The model is parameterized by transition probabilities, and MLEs are obtained in closed form.

Usage

fit_cat(
  y,
  order = 1,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  na_action = c("fail", "complete", "marginalize", "em"),
  em_max_iter = 100,
  em_tol = 1e-06,
  em_epsilon = 1e-08,
  em_safeguard = TRUE,
  em_verbose = FALSE
)

Arguments

Details

For AD(p), the model decomposes as:

P(Y_1, \ldots, Y_n) = P(Y_1, \ldots, Y_p) \times \prod_{k=p+1}^{n} P(Y_k | Y_{k-p}, \ldots, Y_{k-1})

MLEs are computed as empirical proportions:

• Marginal/joint probabilities: count / N

• Transition probabilities: conditional count / marginal count

Empty cells receive probability 0 (if denominator is also 0).

When na_action = "em", fit_cat() dispatches to em_cat. In that case, em_safeguard controls step-halving protection against likelihood-decreasing updates, and returned log_l/AIC/BIC/cell_counts are synchronized via a final E-step under the returned parameters. For order = 2, na_action = "em" is not available and errors explicitly; use na_action = "marginalize".

Value

A list of class "cat_fit" containing:

References

Xie, Y. and Zimmerman, D. L. (2013). Antedependence models for nonstationary categorical longitudinal data with ignorable missingness: likelihood-based inference. Statistics in Medicine, 32, 3274-3289.

Examples

# Simulate binary AD(1) data
set.seed(123)
y <- simulate_cat(n_subjects = 100, n_time = 5, order = 1, n_categories = 2)

# Fit model
fit <- fit_cat(y, order = 1)
print(fit)

# Compare orders
fit0 <- fit_cat(y, order = 0)
fit1 <- fit_cat(y, order = 1)
fit2 <- fit_cat(y, order = 2)
c(AIC_0 = fit0$aic, AIC_1 = fit1$aic, AIC_2 = fit2$aic)

# EM fit with missing data
y_miss <- y
y_miss[sample(length(y_miss), size = round(0.15 * length(y_miss)))] <- NA
fit_em <- fit_cat(
  y_miss,
  order = 1,
  na_action = "em",
  em_max_iter = 80,
  em_tol = 1e-6
)
fit_em$settings$n_iter
fit_em$settings$cell_counts_type

Fit Gaussian Antedependence Model by Maximum Likelihood

Description

Fits an AD(0), AD(1), or AD(2) model for Gaussian longitudinal data by maximum likelihood. Missing values can be handled by complete-case deletion or by EM (see em_gau for an explicit EM wrapper).

Usage

fit_gau(
  y,
  order = 1,
  blocks = NULL,
  na_action = c("fail", "complete", "em"),
  estimate_mu = TRUE,
  em_max_iter = 100,
  em_tol = 1e-06,
  em_verbose = FALSE,
  ...
)

Arguments

Details

For missing data with na_action = "em", AD orders 0 and 1 are the primary production path. AD order 2 is available, but the current EM implementation uses simplified second-order updates and should be treated as provisional for high-stakes inference.

For observed-data likelihood evaluation under MAR without fitting, use logL_gau with na_action = "marginalize". In contrast, fit_gau handles missingness via complete-case fitting (na_action = "complete") or EM (na_action = "em").

Value

A list with components including mu, phi, sigma, tau, log_l, n_obs, n_missing.

See Also

em_gau, fit_cat, fit_inad

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 30, n_time = 5, order = 1, phi = 0.3)
fit <- fit_gau(y, order = 1)
fit$log_l

y_miss <- y
y_miss[1, 2] <- NA
fit_em <- fit_gau(y_miss, order = 1, na_action = "em", em_max_iter = 20)
fit_em$settings$na_action

Fit INAD Antedependence Model by Maximum Likelihood

Description

Fits INAD models by maximum likelihood.

Usage

fit_inad(
  y,
  order = 1,
  thinning = c("binom", "pois", "nbinom"),
  innovation = c("pois", "bell", "nbinom"),
  blocks = NULL,
  max_iter = 50,
  tol = 1e-06,
  verbose = FALSE,
  init_alpha = NULL,
  init_theta = NULL,
  init_tau = 0.4,
  init_nb_inno_size = 1,
  nb_inno_size_ub = 50,
  na_action = c("fail", "complete", "marginalize"),
  constrain = NULL,
  newton_polish = TRUE
)

Arguments

Details

No fixed effect: time-separable optimization using logL_inad_i, with the thinning coefficients, innovation parameter, and negative binomial size of each time point optimized jointly from several starting values.

Fixed effect: block coordinate descent using nloptr BOBYQA, updating tau, alpha, theta, and nb_inno_size if needed.

Parameter spaces: thinning coefficients satisfy 0 \le \alpha \le 1 for binomial thinning and \alpha \ge 0 for Poisson and negative binomial thinning, as in em_inad. Internally the optimizers use an upper limit of twice \sum_s y_{si} / \sum_s y_{s,i-j} for the latter, which cannot be active at a maximum because every stationary point satisfies \alpha_{i,i-j} \le \sum_s y_{si} / \sum_s y_{s,i-j}. A coefficient whose predecessor counts are all zero is not identifiable and is set to 0. The negative binomial innovation sizes are bounded above by the numerical cap nb_inno_size_ub; kkt_inad reports estimates on that cap and checks the first-order conditions of a fit.

Value

A list of class "inad_fit" containing:

alpha

Estimated antedependence parameter(s)

theta

Estimated innovation parameter(s)

tau

Estimated block effects (if applicable)

nb_inno_size

Estimated innovation NB size parameter(s), when innovation = "nbinom"

log_l

Maximized log-likelihood

aic

Akaike information criterion

bic

Bayesian information criterion

n_params

Number of free parameters

convergence

Convergence code

settings

Model and fitting settings

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 60,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
fit <- fit_inad(y, order = 1, thinning = "binom", innovation = "pois", max_iter = 20)
fit$log_l

Fitted Values and Residuals for Gaussian AD Fits

Description

Fitted Values and Residuals for Gaussian AD Fits

Usage

## S3 method for class 'gau_fit'
fitted(object, y, ...)

## S3 method for class 'gau_fit'
residuals(object, y, ...)

Arguments

Value

A numeric matrix of fitted values or residuals.

First-Order Optimality Check for a Fitted INAD Model

Description

Evaluates the observed-data score of a fitted INAD model at its estimates, computed exactly by Louis' identity (Fisher's identity for the score), and checks the first-order (Karush-Kuhn-Tucker) conditions for a maximum over the parameter space.

Usage

kkt_inad(fit, y = NULL, blocks = NULL, tol = 1e-06)

Arguments

Details

The conditions are a zero score for every interior estimate, a nonpositive one-sided score for a thinning coefficient or an innovation parameter at zero (an innovation mean on the boundary), a nonnegative one-sided score for a binomial thinning coefficient at one, and a nonnegative score for a negative binomial innovation size at its numerical cap (nb_inno_size_ub). At \alpha = 0, and at binomial \alpha = 1, the thinned count is degenerate and the score is the exact one-sided derivative of the observed-data log-likelihood.

Interior stationarity is checked with the observed information I of the interior parameters, which must be positive definite. The Newton decrement s^\top I^{-1} s, where s is the score, is about twice the log-likelihood increase that a further Newton step would give, and must not exceed tol; each interior parameter must also satisfy s_j^2 / I_{jj} \le tol, about twice the increase available by moving it alone (this is implied by the decrement condition). A boundary estimate whose score points into the parameter space passes only if the largest log-likelihood increase available by moving it alone into the interior, found by a line search on the observed-data log-likelihood, is at most tol / 2. Thinning coefficients whose predecessor counts are all zero are not identifiable; they are reported with status "not identifiable" and are not checked.

Satisfying these conditions shows that the fit is a stationary point that meets the boundary conditions, to the stated tolerance. It does not show that the fit is the global maximum.

Value

An object of class "kkt_inad": a list with table (one row per free parameter, with columns param, est, score, status, check, the quantity compared with tol, and ok), info_pd, min_info_eigen (the smallest eigenvalue of the interior observed information), newton_decrement, max_abs_interior_score, n_lower, n_upper, n_size_cap, and kkt_ok.

See Also

fit_inad, em_inad, ci_inad

Examples

set.seed(1)
y <- simulate_inad(n_subjects = 60, n_time = 5, order = 1, thinning = "binom",
                   innovation = "pois", alpha = 0.3, theta = 2)
fit <- em_inad(y, order = 1, thinning = "binom", innovation = "pois",
               newton_polish = TRUE)
kkt_inad(fit)

Labor Force Longitudinal Categorical Data (Table 1)

Description

Five-year employment-status sequences reconstructed from Table 1 in the labor-force example used in Xie and Zimmerman score/Wald antedependence testing work. Category coding is 1 = employed, 2 = unemployed.

Usage

labor_force_cat

Format

A list with five components:

y

integer matrix of dimension N by 5 containing expanded subject-level sequences

counts

data frame with columns Y1, Y2, Y3, Y4, Y5, Count

n_categories

number of categories (2)

time

integer vector of calendar years (1967:1971)

status_labels

character vector c("employed", "unemployed")

Source

Table 1 (labor-force example) from: Xie, Y. and Zimmerman, D. L. (2013). Score and Wald tests for antedependence in categorical longitudinal data.

Log-Likelihood for Categorical AD Models

Description

Evaluates the log-likelihood of an AD(p) model for categorical longitudinal data at given parameter values.

Usage

logL_cat(
  y,
  order,
  marginal,
  transition = NULL,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  na_action = c("fail", "complete", "marginalize")
)

Arguments

Details

The log-likelihood for AD(p) decomposes into contributions from initial time points and transition time points.

For order 0 (independence), the log-likelihood is the sum of log marginal probabilities at each time point.

Parameter structure for marginal:

• Order 0: List with elements t1, t2, ..., tn, each a vector of length c

• Order 1: List with element t1 (vector of length c)

• Order 2: List with t1 (vector), t2_given_1to1 (c x c matrix)

Parameter structure for transition:

• Order 0: Not used (NULL or empty list)

• Order 1: List with elements t2, t3, ..., tn, each c x c matrix

• Order 2: List with elements t3, t4, ..., tn, each c x c x c array

Value

Scalar log-likelihood value.

References

Xie, Y. and Zimmerman, D. L. (2013). Antedependence models for nonstationary categorical longitudinal data with ignorable missingness: likelihood-based inference. Statistics in Medicine, 32, 3274-3289.

Examples

set.seed(1)
y <- simulate_cat(n_subjects = 40, n_time = 5, order = 1, n_categories = 3)
fit <- fit_cat(y, order = 1, n_categories = 3)
logL_cat(
  y = y,
  order = 1,
  marginal = fit$marginal,
  transition = fit$transition,
  n_categories = 3
)

Log-Likelihood for Gaussian AD Models

Description

Computes the log-likelihood for Gaussian antedependence models of order 0, 1, or 2. Supports missing data under MAR assumption via na_action parameter.

Usage

logL_gau(
  y,
  order,
  mu = NULL,
  phi = NULL,
  sigma = NULL,
  blocks = NULL,
  tau = 0,
  na_action = c("fail", "complete", "marginalize")
)

Arguments

Details

For complete data (no NA), all three na_action options give the same result.

For missing data:

• marginalize: Uses MVN marginalization to compute P(Y_obs). This is the correct observed-data likelihood for MAR missing data.

• complete: Removes subjects with any missing values. May lose information.

• fail: Stops with error. Useful to ensure no missing data present.

Value

Scalar log-likelihood value.

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 30, n_time = 5, order = 1, phi = 0.3)
fit <- fit_gau(y, order = 1)
logL_gau(y, order = 1, mu = fit$mu, phi = fit$phi, sigma = fit$sigma)

Log-Likelihood for INAD Models (with Missing Data Support)

Description

If blocks is NULL, this computes the log likelihood as the sum of per time contributions from logL_inad_i for computational convenience.

Usage

logL_inad(
  y,
  order = 1,
  thinning = c("binom", "pois", "nbinom"),
  innovation = c("pois", "bell", "nbinom"),
  alpha,
  theta,
  nb_inno_size = NULL,
  blocks = NULL,
  tau = 0,
  na_action = c("fail", "complete", "marginalize")
)

Arguments

Value

A scalar log likelihood.

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 60,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
fit <- fit_inad(y, order = 1, thinning = "binom", innovation = "pois", max_iter = 20)
logL_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = fit$alpha,
  theta = fit$theta
)

INAD Log-Likelihood Contribution at Time i

Description

Returns the time i contribution, summed over subjects, for the no fixed effect model.

Usage

logL_inad_i(
  y,
  i,
  order = 1,
  thinning = c("binom", "pois", "nbinom"),
  innovation = c("pois", "bell", "nbinom"),
  alpha,
  theta,
  nb_inno_size = NULL
)

Arguments

Value

A scalar log likelihood contribution for time i.

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 50,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
fit <- fit_inad(y, order = 1, thinning = "binom", innovation = "pois", max_iter = 20)
logL_inad_i(
  y = y,
  i = 3,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = fit$alpha,
  theta = fit$theta
)

Log-Likelihood for Antedependence Model Fits

Description

Extracts the maximized log-likelihood from a fitted antedependence model and returns a "logLik" object compatible with AIC and BIC.

Usage

## S3 method for class 'ad_fit'
logLik(object, ...)

## S3 method for class 'gau_fit'
logLik(object, ...)

## S3 method for class 'cat_fit'
logLik(object, ...)

## S3 method for class 'inad_fit'
logLik(object, ...)

Arguments

Value

A scalar of class "logLik" with attributes df and nobs.

Number of Subjects for Antedependence Model Fits

Description

Number of Subjects for Antedependence Model Fits

Usage

## S3 method for class 'ad_fit'
nobs(object, ...)

## S3 method for class 'gau_fit'
nobs(object, ...)

## S3 method for class 'cat_fit'
nobs(object, ...)

## S3 method for class 'inad_fit'
nobs(object, ...)

Arguments

Value

Integer scalar: number of subjects.

Compute Intervenor-Adjusted Partial Correlation Matrix

Description

Computes the partial correlation between Y[i] and Y[j] adjusting for the "intervenor" variables ⁠Y[i+1], ..., Y[j-1]⁠. Under an antedependence model of order p, partial correlations for |i-j| > p should be approximately zero.

Usage

partial_corr(y, test = FALSE, n_digits = 3)

Arguments

Details

The intervenor-adjusted partial correlation between Y[i] and Y[j] (i < j) is computed as the correlation between the residuals from regressing Y[i] and Y[j] on the intervenor set ⁠Y[i+1], ..., Y[j-1]⁠.

For adjacent time points (|i-j| = 1), the partial correlation equals the ordinary correlation since there are no intervenors.

The diagonal of both returned matrices contains variances (not correlations). This keeps scale information available alongside correlation structure.

The significance test uses an approximate threshold of 2/sqrt(n_eff), which corresponds roughly to a 95% confidence bound under normality. This is a rough screening tool, not a formal hypothesis test.

Value

A list with components:

References

Zimmerman, D. L. and Nunez-Anton, V. (2009). Antedependence Models for Longitudinal Data. CRC Press.

See Also

plot_prism for visual diagnostics

Examples

data("bolus_inad")
pc <- partial_corr(bolus_inad$y, test = TRUE)

# View partial correlations (lower triangle)
pc$partial_correlation

# Extract variances from the diagonal
variances <- diag(pc$partial_correlation)

# Check which are "significant" (rough screen for AD order)
pc$significant

Plot marginal probabilities of a categorical AD fit

Description

Displays a stacked bar chart of the estimated marginal category probabilities at each time point.

Usage

## S3 method for class 'cat_fit'
plot(x, ...)

Arguments

Value

x, invisibly.

Examples

set.seed(1)
y <- simulate_cat(n_subjects = 80, n_time = 5, order = 1, n_categories = 2)
fit <- fit_cat(y, order = 1, n_categories = 2)
plot(fit)

Plot estimated parameters of a Gaussian AD fit

Description

Produces a multi-panel base-graphics plot showing the estimated mean (\mu), innovation standard deviation (\sigma), and (for AD(1) or AD(2)) the antedependence coefficients (\phi) over time.

Usage

## S3 method for class 'gau_fit'
plot(x, ...)

Arguments

Value

x, invisibly.

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 40, n_time = 6, order = 1, phi = 0.4)
fit <- fit_gau(y, order = 1)
plot(fit)

Plot estimated parameters of an INAD fit

Description

Produces a two-panel base-graphics plot of the estimated thinning parameters (\alpha) and innovation means (\theta) over time.

Usage

## S3 method for class 'inad_fit'
plot(x, ...)

Arguments

Value

x, invisibly.

Examples

set.seed(1)
y <- simulate_inad(n_subjects = 60, n_time = 5, order = 1,
                   thinning = "binom", innovation = "pois",
                   alpha = 0.3, theta = 2)
fit <- fit_inad(y, order = 1, thinning = "binom",
                innovation = "pois", max_iter = 20)
plot(fit)

Heatmap plot for a partial_corr object

Description

Displays the partial correlation matrix as a colour heatmap using base graphics. Upper triangle shows ordinary correlations; lower triangle shows intervenor-adjusted partial correlations; the diagonal shows standardised variances.

Usage

## S3 method for class 'partial_corr'
plot(x, ...)

Arguments

Value

x, invisibly.

Examples

data("bolus_inad")
pc <- partial_corr(bolus_inad$y)
plot(pc)

PRISM Plot (Partial Residual Intervenor Scatterplot Matrix)

Description

Creates a matrix of scatterplots for diagnosing antedependence structure. The upper triangle shows ordinary scatterplots of Y[i] vs Y[j]. The lower triangle shows PRISM plots: residuals from regressing Y[i] and Y[j] on the intervenor variables ⁠Y[i+1], ..., Y[j-1]⁠.

Usage

plot_prism(
  y,
  time_labels = NULL,
  pch = 20,
  cex = 0.6,
  col_upper = "steelblue",
  col_lower = "firebrick",
  main = "PRISM Diagnostic Plot"
)

Arguments

Details

Under an antedependence model of order p, the partial correlation between Y[i] and Y[j] given the intervenors should be zero when |i-j| > p. This means PRISM plots in the lower triangle should show no association for lags greater than p.

Interpretation:

• Upper triangle: Shows marginal associations between time points

• Lower triangle (PRISM): Shows conditional associations after removing effects of intervenor variables

• If AD(1) holds: Only the first sub-diagonal of lower triangle should show association

• If AD(2) holds: First two sub-diagonals should show association

Value

Invisibly returns NULL. Called for side effect (plotting).

References

Zimmerman, D. L. and Nunez-Anton, V. (2009). Antedependence Models for Longitudinal Data. CRC Press. Chapter 2.

See Also

partial_corr for numerical partial correlations

Examples

data("bolus_inad")
plot_prism(bolus_inad$y)

# With custom labels
plot_prism(bolus_inad$y, time_labels = paste0("Hour ", seq(0, 44, by = 4)))

Profile Plot (Spaghetti Plot) for Longitudinal Data

Description

Creates a profile plot showing individual subject trajectories with overlaid mean trajectory and standard deviation bands.

Usage

plot_profile(
  y,
  time_points = NULL,
  blocks = NULL,
  block_labels = NULL,
  title = "Profile Plot",
  xlab = "Time",
  ylab = "Measurement",
  ylim = NULL,
  show_sd = TRUE,
  individual_alpha = 0.3,
  individual_color = "grey50",
  mean_color = "blue",
  sd_color = "red",
  mean_lwd = 2
)

Arguments

Details

This function provides a quick visual summary of longitudinal data showing:

• Individual subject trajectories (light grey lines)

• Mean trajectory across subjects (bold colored line)

• +/- 1 standard deviation bands (error bars)

When blocks is provided, the plot is faceted by block membership, allowing comparison of trajectories across treatment groups or other strata.

Value

A ggplot2 object (invisibly). Called primarily for side effect (plotting).

Examples

data("bolus_inad")

# Basic profile plot
plot_profile(bolus_inad$y)

# With block stratification
plot_profile(bolus_inad$y, blocks = bolus_inad$blocks,
             block_labels = c("2mg", "1mg"))

# Customized
plot_profile(bolus_inad$y,
             time_points = seq(0, 44, by = 4),
             title = "Bolus Counts Over Time",
             xlab = "Hours", ylab = "Count")

Predict from an INAD fit

Description

Generates plug-in predictive means or Monte Carlo sample paths from a fitted integer-valued antedependence model. Version 1 supports order-1 INAD fits and uses fitted maximum-likelihood parameters without integrating over parameter uncertainty.

Usage

## S3 method for class 'inad_fit'
predict(
  object,
  newdata = NULL,
  blocks = NULL,
  h = 1L,
  type = c("mean", "sample"),
  n_sims = 1000L,
  seed = NULL,
  ...
)

Arguments

Details

Predictions are conditional on the fitted parameter estimates, following the usual plug-in convention used by many predict() methods. They do not account for parameter-estimation uncertainty. When newdata is supplied, the forecast starts at ncol(newdata) + 1; therefore ncol(newdata) + h must not exceed the fitted time grid.

Value

For type = "mean", an n_test x h numeric matrix. For type = "sample", an n_test x h x n_sims integer array.

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 30,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
fit <- fit_inad(y, order = 1, thinning = "binom", innovation = "pois")
predict(fit, h = 1, type = "mean", n_sims = 20, seed = 1)

Print Method for Categorical BIC Order Selection

Description

Print Method for Categorical BIC Order Selection

Usage

## S3 method for class 'bic_order_cat'
print(x, digits = 4, ...)

Arguments

Value

The input object, invisibly.

Print method for cat_ci objects

Description

Print method for cat_ci objects

Usage

## S3 method for class 'cat_ci'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for cat_fit objects

Description

Print method for cat_fit objects

Usage

## S3 method for class 'cat_fit'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for cat_lrt objects

Description

Print method for cat_lrt objects

Usage

## S3 method for class 'cat_lrt'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print Method for Categorical Order Test Collections

Description

Print Method for Categorical Order Test Collections

Usage

## S3 method for class 'cat_order_tests'
print(x, digits = 4, ...)

Arguments

Value

The input object, invisibly.

Print method for BIC order selection

Description

Print method for BIC order selection

Usage

## S3 method for class 'gau_bic_order'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for AD confidence intervals

Description

Print method for AD confidence intervals

Usage

## S3 method for class 'gau_ci'
print(x, ...)

Arguments

Value

The input object, invisibly.

Print method for AD contrast test

Description

Print method for AD contrast test

Usage

## S3 method for class 'gau_contrast_test'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print Method for gau_fit Objects

Description

Print Method for gau_fit Objects

Usage

## S3 method for class 'gau_fit'
print(x, ...)

Arguments

Value

The input object, invisibly.

Print method for AD homogeneity test

Description

Print method for AD homogeneity test

Usage

## S3 method for class 'gau_homogeneity_test'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for AD mean test

Description

Print method for AD mean test

Usage

## S3 method for class 'gau_mean_test'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for AD order test

Description

Print method for AD order test

Usage

## S3 method for class 'gau_order_test'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for homogeneity_tests_inad

Description

Print method for homogeneity_tests_inad

Usage

## S3 method for class 'homogeneity_tests_inad'
print(x, digits = 4, ...)

Arguments

Value

Invisibly returns x.

Print method for INAD confidence intervals

Description

Print method for INAD confidence intervals

Usage

## S3 method for class 'inad_ci'
print(x, ...)

Arguments

Value

The input object, invisibly.

Print Method for INAD Model Fits

Description

Print Method for INAD Model Fits

Usage

## S3 method for class 'inad_fit'
print(x, digits = 4, ...)

Arguments

Value

The input object, invisibly.

Print Method for Categorical Stationarity Test Collections

Description

Print Method for Categorical Stationarity Test Collections

Usage

## S3 method for class 'stationarity_tests_cat'
print(x, digits = 4, ...)

Arguments

Value

The input object, invisibly.

Print method for test_homogeneity_inad

Description

Print method for test_homogeneity_inad

Usage

## S3 method for class 'test_homogeneity_inad'
print(x, digits = 4, ...)

Arguments

Value

Invisibly returns x.

100km Race Split-Time Data

Description

Split times (in minutes) from a 100km race example with 10 consecutive sections. This is continuous-response longitudinal data suitable for Gaussian AD modeling.

Usage

race_100km

Format

A list with five components:

y

numeric matrix of dimension N by 10 containing section split times

age

numeric vector of subject ages (may include missing values)

subject

integer subject identifiers

time

integer vector of section indices (1:10)

section_labels

character vector of split-time column names

Source

Combined table extracted from textbook race example data, stored in data-raw/external/100km_race_combined_extracted.csv.

Run All Homogeneity Tests for INAD

Description

Convenience function to run all three homogeneity tests at once and return a summary.

Usage

run_homogeneity_tests_inad(
  y,
  blocks,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  ...
)

Arguments

Value

A list with class "homogeneity_tests_inad" containing results for all three tests and a summary table.

Examples

data("bolus_inad")
# First four occasions, to keep the example short
tests <- run_homogeneity_tests_inad(bolus_inad$y[, 1:4], bolus_inad$blocks,
                                     order = 1)
print(tests)

Run All Pairwise Order Tests for Categorical AD

Description

Performs sequential likelihood ratio tests for AD orders 0 vs 1, 1 vs 2, etc.

Usage

run_order_tests_cat(
  y,
  max_order = 2,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  test = c("lrt", "score", "mlrt", "wald")
)

Arguments

Details

This function performs forward selection: starting from order 0, it tests whether increasing the order provides significant improvement. The selected order is the highest order where the test was significant (at alpha = 0.05).

Value

A list containing:

Examples

y <- simulate_cat(200, 6, order = 1, n_categories = 2)
result <- run_order_tests_cat(y, max_order = 2)
print(result)
cat("Selected order:", result$selected_order, "\n")

Run All Stationarity-Related Tests for Categorical AD

Description

Performs tests for time-invariance and stationarity constraints. For order = 1, the stationarity test corresponds to strict stationarity; for order > 1, it tests marginal-constancy plus time-invariant transitions. Currently supports complete data only.

Usage

run_stationarity_tests_cat(
  y,
  order = 1,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  test = c("lrt", "score", "mlrt")
)

Arguments

Value

A list containing:

time_invariance

Result of test_timeinvariance_cat

stationarity

Result of test_stationarity_cat

table

Summary data frame

Examples

y <- simulate_cat(200, 6, order = 1, n_categories = 2)
result <- run_stationarity_tests_cat(y, order = 1)
print(result)

Run All Stationarity-Related Tests for Gaussian AD

Description

Runs a standard set of stationarity constraints for Gaussian AD models.

Usage

run_stationarity_tests_gau(
  y,
  order = 1L,
  blocks = NULL,
  verbose = FALSE,
  max_iter = 2000L,
  rel_tol = 1e-08,
  ...
)

Arguments

Value

A list with class "stationarity_tests_gau" containing:

fit_unconstrained

Unconstrained Gaussian AD fit

tests

Named list of test_stationarity_gau results

summary

Summary table of all constraints

See Also

test_stationarity_gau

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 80, n_time = 6, order = 1, phi = 0.4, sigma = 1)
out <- run_stationarity_tests_gau(y, order = 1, verbose = FALSE)
out$summary

Run All Stationarity-Related Tests for INAD

Description

Runs the standard battery of time-invariance tests against a common unconstrained fit. For order 1 the constraints are "alpha", "theta", and "both"; with negative binomial innovations they are "alpha", "innovation" (theta and kappa jointly), and "all". For order 2 they are "alpha1", "alpha2", "alpha", "theta" (or "innovation"), and "all".

Usage

run_stationarity_tests_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  blocks = NULL,
  verbose = FALSE,
  ...
)

Arguments

Value

A list with class "stationarity_tests_inad".

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 15,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
out <- run_stationarity_tests_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  verbose = FALSE,
  max_iter = 8
)
out$summary

Simulate Categorical Antedependence Series

Description

Generate simulated longitudinal categorical data from an AD(p) model with specified transition probabilities.

Usage

simulate_cat(
  n_subjects,
  n_time,
  order = 1,
  n_categories = 2,
  marginal = NULL,
  transition = NULL,
  blocks = NULL,
  homogeneous = TRUE,
  seed = NULL
)

Arguments

Details

Data are simulated sequentially:

1. For k = 1: Draw Y(1) from marginal distribution

2. For k = 2 to p: Draw Y(k) conditional on Y(1), ..., Y(k-1)

3. For k = p+1 to n: Draw Y(k) conditional on Y(k-p), ..., Y(k-1)

Parameter structure for marginal:

• Order 0: List with elements t1, t2, ..., tn, each a vector of length c summing to 1

• Order 1: List with element t1 (vector of length c)

• Order 2: List with t1 (vector), t2_given_1to1 (c x c matrix where rows represent conditioning values and columns represent outcomes)

Parameter structure for transition:

• Order 0: Not used (NULL)

• Order 1: List with elements t2, t3, ..., tn, each c x c matrix where rows are previous values and columns are current values (rows sum to 1)

• Order 2: List with elements t3, t4, ..., tn, each c x c x c array where first two indices are conditioning values and third is outcome

Value

An integer matrix of class c("cat_sim", "matrix") with n_subjects rows and n_time columns, where each entry is a category code from 1 to n_categories. Use as.ts(y) (or ts(t(y))) to convert to an mts time-series object where rows correspond to time points.

References

Xie, Y. and Zimmerman, D. L. (2013). Antedependence models for nonstationary categorical longitudinal data with ignorable missingness: likelihood-based inference. Statistics in Medicine, 32, 3274-3289.

Examples

y <- simulate_cat(n_subjects = 30, n_time = 5, order = 1, n_categories = 3, seed = 1)
dim(y)

Simulate Gaussian Antedependence Series

Description

Generate longitudinal continuous data from a Gaussian antedependence (AD) model of order 0, 1, or 2 using a conditional regression on predecessors.

Usage

simulate_gau(
  n_subjects,
  n_time,
  order = 1L,
  mu = NULL,
  phi = NULL,
  sigma = NULL,
  blocks = NULL,
  tau = 0,
  seed = NULL
)

Arguments

Details

For order = 0, each time point is generated independently as Y[, t] = mu[t] + tau[block] + eps, with eps ~ N(0, sigma[t]^2).

For order = 1, for t >= 2: ⁠Y[, t] = m_t + phi[t] * (Y[, t - 1] - m_{t-1}) + eps_t⁠, where m_t = mu[t] + tau[block] and eps_t ~ N(0, sigma[t]^2).

For order = 2, for t >= 3: ⁠Y[, t] = m_t + phi[1, t] * (Y[, t - 1] - m_{t-1}) + phi[2, t] * (Y[, t - 2] - m_{t-2}) + eps_t⁠.

If blocks is provided, each subject s belongs to a block and receives a mean shift tau[blocks[s]]. tau[1] is forced to 0.

Value

A numeric matrix of class c("gau_sim", "matrix") with n_subjects rows and n_time columns. Each row is one subject's simulated trajectory. Use as.ts(y) (or ts(t(y))) to convert to an mts time-series object where rows correspond to time points.

Examples

y <- simulate_gau(
  n_subjects = 20,
  n_time = 6,
  order = 1,
  phi = 0.4,
  seed = 42
)
dim(y)

Simulate INAD Antedependence Series

Description

Generate longitudinal count data from an INAD model using a thinning operator and an innovation distribution.

Usage

simulate_inad(
  n_subjects,
  n_time,
  order = 1L,
  thinning = c("binom", "pois", "nbinom"),
  innovation = c("pois", "bell", "nbinom"),
  alpha = NULL,
  theta = NULL,
  nb_inno_size = NULL,
  blocks = NULL,
  tau = 0,
  seed = NULL
)

Arguments

Details

Time 1 observations are generated from the innovation distribution alone. For times 2 to n_time, counts are generated as thinning of previous counts plus independent innovations. When order = 0, all time points are generated from the innovation distribution and the thinning operator and alpha are ignored.

If blocks is provided, innovations include a block effect. For Poisson and negative binomial innovations, the innovation mean is theta[t] + tau[blocks[i]]. For Bell innovations, the innovation mean is theta[t] * exp(theta[t]) + tau[blocks[i]].

Value

An integer matrix of class c("inad_sim", "matrix") with n_subjects rows and n_time columns. Each row is one subject's simulated count trajectory. Use as.ts(y) (or ts(t(y))) to convert to an mts time-series object where rows correspond to time points.

Examples

y <- simulate_inad(
  n_subjects = 20,
  n_time = 6,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2,
  seed = 42
)
dim(y)

Summary method for cat_ci objects

Description

Summary method for cat_ci objects

Usage

## S3 method for class 'cat_ci'
summary(object, ...)

Arguments

Value

A data frame summarizing all CIs

Summary method for cat_fit objects

Description

Summary method for cat_fit objects

Usage

## S3 method for class 'cat_fit'
summary(object, ...)

Arguments

Value

object, invisibly.

Summary method for gau_ci objects

Description

Summary method for gau_ci objects

Usage

## S3 method for class 'gau_ci'
summary(object, ...)

Arguments

Value

A data frame stacking available confidence intervals with columns param, component, est, se, lower, upper, and level. Returns NULL if no intervals are available.

Summary method for gau_fit objects

Description

Prints a structured table of model-fit statistics and estimated parameters.

Usage

## S3 method for class 'gau_fit'
summary(object, ...)

Arguments

Value

object, invisibly.

Summary method for inad_ci objects

Description

Summary method for inad_ci objects

Usage

## S3 method for class 'inad_ci'
summary(object, ...)

Arguments

Value

A data frame stacking available confidence intervals with columns param, component, est, se, lower, upper, and level. Returns NULL if no intervals are available.

Summary Method for INAD Model Fits

Description

Summary method for inad_fit objects

Usage

## S3 method for class 'inad_fit'
summary(object, ...)

Arguments

Value

object, invisibly.

Summary method for partial_corr objects

Description

Prints a lag-by-lag summary of intervenor-adjusted partial correlations, reporting the mean and range at each lag together with the number of pairs that exceed the rough significance threshold 2/\sqrt{n}.

Usage

## S3 method for class 'partial_corr'
summary(object, ...)

Arguments

Value

Invisibly returns a data frame with one row per lag containing columns lag, n_pairs, mean_abs, min, and max.

See Also

partial_corr, plot.partial_corr

Test Linear Hypotheses on the Mean Under Antedependence

Description

Tests the null hypothesis C * mu = c for a specified contrast matrix C and vector c, under an AD(p) covariance structure. This implements Theorem 7.2 of Zimmerman & Núñez-Antón (2009).

Usage

test_contrast_gau(y, C, c = NULL, p = 1L)

Arguments

Details

The Wald test statistic (Theorem 7.2) is:

(C\bar{Y} - c)^T (C \hat{\Sigma} C^T)^{-1} (C\bar{Y} - c)

where \hat{\Sigma} is the REML estimator of the covariance matrix under the AD(p) model.

Common examples include:

• Testing if mean is constant: C is the first-difference matrix

• Testing for linear trend: C tests deviations from linearity

Value

A list with class gau_contrast_test containing:

method

Inference method used ("wald").

C

Contrast matrix

c

Right-hand side vector

mu_hat

Estimated mean vector

contrast_est

Estimated value of C * mu

statistic

Wald test statistic

df

Degrees of freedom (number of contrasts)

p_value

P-value from chi-square distribution

References

Zimmerman, D.L. and Núñez-Antón, V. (2009). Antedependence Models for Longitudinal Data. Chapman & Hall/CRC. Chapter 7.

Examples

y <- simulate_gau(n_subjects = 50, n_time = 5, order = 1)

# Test if mean is constant (all differences = 0)
# C is 4x5 matrix of first differences
C <- matrix(0, nrow = 4, ncol = 5)
for (i in 1:4) {
  C[i, i] <- 1
  C[i, i+1] <- -1
}
test <- test_contrast_gau(y, C = C, p = 1)
print(test)

Likelihood Ratio Test for Homogeneity Across Groups (Categorical AD Data)

Description

Tests whether multiple groups share the same transition probability parameters in a categorical antedependence model.

Usage

test_homogeneity_cat(
  y = NULL,
  blocks = NULL,
  order = 1,
  n_categories = NULL,
  fit_null = NULL,
  fit_alt = NULL,
  test = c("lrt", "score", "mlrt")
)

Arguments

Details

The null hypothesis is that all G groups share the same transition probability parameters:

H_0: \pi^{(1)} = \pi^{(2)} = \ldots = \pi^{(G)}

The alternative hypothesis allows each group to have its own parameters.

The degrees of freedom are:

df = (G-1) \times k

where G is the number of groups and k is the number of free parameters per population.

Value

A list of class "cat_lrt" containing:

method

Inference method used: one of "lrt", "score", "mlrt", or "wald".

lrt_stat

Likelihood ratio test statistic

df

Degrees of freedom

p_value

P-value from chi-square distribution

fit_null

Fitted homogeneous model (H0)

fit_alt

Fitted heterogeneous model (H1)

n_groups

Number of groups

table

Summary data frame

References

Xie, Y. and Zimmerman, D. L. (2013). Antedependence models for nonstationary categorical longitudinal data with ignorable missingness: likelihood-based inference. Statistics in Medicine, 32, 3274-3289.

See Also

fit_cat, test_order_cat

Examples

# Simulate data with different transition probabilities for two groups
set.seed(123)
marg1 <- list(t1 = c(0.7, 0.3))
marg2 <- list(t1 = c(0.4, 0.6))
trans1 <- list(t2 = matrix(c(0.9, 0.1, 0.2, 0.8), 2, byrow = TRUE),
               t3 = matrix(c(0.9, 0.1, 0.2, 0.8), 2, byrow = TRUE))
trans2 <- list(t2 = matrix(c(0.5, 0.5, 0.5, 0.5), 2, byrow = TRUE),
               t3 = matrix(c(0.5, 0.5, 0.5, 0.5), 2, byrow = TRUE))

y1 <- simulate_cat(100, 3, order = 1, n_categories = 2,
                   marginal = marg1, transition = trans1)
y2 <- simulate_cat(100, 3, order = 1, n_categories = 2,
                   marginal = marg2, transition = trans2)
y <- rbind(y1, y2)
blocks <- c(rep(1, 100), rep(2, 100))

# Test homogeneity
test <- test_homogeneity_cat(y, blocks, order = 1)
print(test)

Likelihood Ratio Test for Homogeneity Across Groups (Gaussian AD Data)

Description

Tests the null hypothesis that G groups have the same AD(p) covariance structure against the alternative that they have different AD(p) structures. This implements Theorem 6.6 of Zimmerman & Núñez-Antón (2009).

Usage

test_homogeneity_gau(y, blocks, p = 1L, use_modified = TRUE)

Arguments

Details

The test compares:

• H0: All G groups have the same AD(p) covariance matrix Sigma(theta)

• H1: Groups have different AD(p) covariance matrices Sigma(theta_g)

The likelihood ratio test statistic (Theorem 6.6) involves comparing pooled and within-group RSS values. The degrees of freedom are (G-1)(2n - p)(p + 1)/2.

This test is useful for determining whether a common covariance structure can be assumed across treatment groups before performing mean comparisons.

Value

A list with class gau_homogeneity_test containing:

method

Inference method used ("lrt").

statistic

Test statistic value

statistic_modified

Modified test statistic (if use_modified = TRUE)

df

Degrees of freedom

p_value

P-value from chi-square distribution

p_value_modified

P-value from modified test

G

Number of groups

group_sizes

Sample sizes for each group

order

Antedependence order

References

Zimmerman, D.L. and Núñez-Antón, V. (2009). Antedependence Models for Longitudinal Data. Chapman & Hall/CRC. Section 6.4.

Kenward, M.G. (1987). A method for comparing profiles of repeated measurements. Applied Statistics, 36, 296-308.

See Also

test_order_gau, test_two_sample_gau

Examples

# Simulate data from two groups with same covariance
n1 <- 30
n2 <- 35
y1 <- simulate_gau(n1, n_time = 6, order = 1, phi = 0.5, sigma = 1)
y2 <- simulate_gau(n2, n_time = 6, order = 1, phi = 0.5, sigma = 1)
y <- rbind(y1, y2)
blocks <- c(rep(1, n1), rep(2, n2))

# Test homogeneity
test <- test_homogeneity_gau(y, blocks, p = 1)
print(test)

Likelihood Ratio Test for Homogeneity Across Groups (INAD Data)

Description

Tests hypotheses about parameter equality across treatment or grouping factors in integer-valued antedependence models. Implements the homogeneity testing framework from Section 3.7 of Li & Zimmerman (2026).

Usage

test_homogeneity_inad(
  y,
  blocks,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  test = c("all", "mean", "dependence"),
  fit_pooled = NULL,
  fit_inadfe = NULL,
  fit_hetero = NULL,
  ...
)

Arguments

Details

The function supports three nested model comparisons as described in the paper:

M1 (Pooled): All parameters are common across groups. This corresponds to fitting fit_inad(y, blocks = NULL).

M2 (INADFE): The thinning parameters \alpha are shared across groups, but innovation means differ via block effects \tau. This is the standard INADFE model fitted via fit_inad(y, blocks = blocks).

M3 (Fully Heterogeneous): Both \alpha and \theta parameters can differ across groups. This is fitted by running separate fit_inad calls for each group.

The three test types correspond to:

• "all": H0: M1 vs H1: M3 (complete homogeneity vs complete heterogeneity)

• "mean": H0: M1 vs H1: M2 (test for group differences in means only)

• "dependence": H0: M2 vs H1: M3 (test for group differences in dependence)

Degrees of freedom are computed as the difference in free parameters between the null and alternative models.

Value

A list with class "test_homogeneity_inad" containing:

method

Inference method used ("lrt").

statistic

Test statistic value

lrt_stat

Likelihood ratio test statistic

df

Degrees of freedom

p_value

P-value from chi-square distribution

test

Type of test performed

fit_null

Fitted model under H0

fit_alt

Fitted model under H1

bic_null

BIC under H0

bic_alt

BIC under H1

bic_selected

Which model BIC prefers

table

Summary data frame

References

Li, C. and Zimmerman, D.L. (2026). Integer-valued antedependence models for longitudinal count data. Manuscript under review.

See Also

fit_inad, test_order_inad, test_stationarity_inad

Examples

data("bolus_inad")
y <- bolus_inad$y[, 1:4]  # first four occasions, to keep the example short
blocks <- bolus_inad$blocks

# Test for any group differences (M1 vs M3)
test_all <- test_homogeneity_inad(y, blocks, order = 1, test = "all")
print(test_all)

# Test only for mean differences (M1 vs M2), reusing the pooled fit
test_mean <- test_homogeneity_inad(y, blocks, order = 1, test = "mean",
                                   fit_pooled = test_all$fit_null)
print(test_mean)

# Test for dependence differences given different means (M2 vs M3),
# reusing the INADFE and heterogeneous fits
test_dep <- test_homogeneity_inad(y, blocks, order = 1, test = "dependence",
                                  fit_inadfe = test_mean$fit_alt,
                                  fit_hetero = test_all$fit_alt)
print(test_dep)

One-Sample Test for Mean Structure Under Antedependence

Description

Tests the null hypothesis that the mean vector equals a specified value mu = mu_0 against the alternative mu != mu_0, under an AD(p) covariance structure. This implements Theorem 7.1 of Zimmerman & Núñez-Antón (2009).

Usage

test_one_sample_gau(y, mu0, p = 1L, order = NULL, use_modified = TRUE)

Arguments

Details

The test exploits the AD structure to gain power over tests that don't assume any covariance structure. The likelihood ratio test statistic (Theorem 7.1) is:

N \sum_{i=1}^{n} [\log RSS_i(\mu_0) - \log RSS_i(\hat{\mu})]

where RSS_i(mu) is the residual sum of squares from the regression of Y_i - mu_i on its p predecessors Y_(i-1) - mu_(i-1), ..., Y_(i-p) - mu_(i-p).

The test has n degrees of freedom (one for each component of mu).

Value

A list with class gau_mean_test containing:

method

Inference method used ("lrt").

test_type

"one-sample"

mu0

Hypothesized mean under null

mu_hat

MLE of mean (sample mean)

statistic

Test statistic value

statistic_modified

Modified test statistic (if use_modified = TRUE)

df

Degrees of freedom (n_time)

p_value

P-value from chi-square distribution

p_value_modified

P-value from modified test

order

Antedependence order used

References

Zimmerman, D.L. and Núñez-Antón, V. (2009). Antedependence Models for Longitudinal Data. Chapman & Hall/CRC. Chapter 7.

See Also

test_two_sample_gau, test_order_gau

Examples

# Simulate data with known mean
mu_true <- c(10, 11, 12, 13, 14, 15)
y <- simulate_gau(n_subjects = 50, n_time = 6, order = 1, mu = mu_true)

# Test if mean is zero
test <- test_one_sample_gau(y, mu0 = rep(0, 6), p = 1)
print(test)

# Test if mean equals true value (should not reject)
test2 <- test_one_sample_gau(y, mu0 = mu_true, p = 1)
print(test2)

Likelihood Ratio Test for Antedependence Order (Categorical AD Data)

Description

Tests whether a higher-order AD model provides significantly better fit than a lower-order model for categorical longitudinal data.

Usage

test_order_cat(
  y = NULL,
  order_null = 0,
  order_alt = 1,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  fit_null = NULL,
  fit_alt = NULL,
  test = c("lrt", "score", "mlrt", "wald")
)

Arguments

Details

The likelihood ratio test statistic is:

\lambda = -2[\ell_0 - \ell_1]

where \ell_0 and \ell_1 are the maximized log-likelihoods under the null and alternative hypotheses.

Under H0, \lambda follows a chi-square distribution with degrees of freedom equal to the difference in the number of free parameters.

For testing AD(p) vs AD(p+1), the degrees of freedom are:

df = (c-1)^2 \times c^p \times (n - p - 1)

where c is the number of categories and n is the number of time points.

If y contains missing values and models are fit internally, this function defaults to na_action = "marginalize" for fitting. Score- and Wald-based variants currently require complete data.

Value

A list of class "cat_lrt" containing:

method

Inference method used: one of "lrt", "score", "mlrt", or "wald".

lrt_stat

Likelihood ratio test statistic

df

Degrees of freedom

p_value

P-value from chi-square distribution

fit_null

Fitted model under H0

fit_alt

Fitted model under H1

order_null

Order under null

order_alt

Order under alternative

table

Summary data frame

References

Xie, Y. and Zimmerman, D. L. (2013). Antedependence models for nonstationary categorical longitudinal data with ignorable missingness: likelihood-based inference. Statistics in Medicine, 32, 3274-3289.

See Also

fit_cat, bic_order_cat

Examples

# Simulate AD(1) data
set.seed(123)
y <- simulate_cat(200, 6, order = 1, n_categories = 2)

# Test AD(0) vs AD(1)
test_01 <- test_order_cat(y, order_null = 0, order_alt = 1)
print(test_01)

# Test AD(1) vs AD(2)
test_12 <- test_order_cat(y, order_null = 1, order_alt = 2)
print(test_12)

# Using pre-fitted models
fit0 <- fit_cat(y, order = 0)
fit1 <- fit_cat(y, order = 1)
test_prefitted <- test_order_cat(fit_null = fit0, fit_alt = fit1)

Likelihood Ratio Test for Antedependence Order (Gaussian AD Data)

Description

Tests the null hypothesis that the data follow an AD(p) model against the alternative that they follow an AD(p+q) model, using the likelihood ratio test described in Theorem 6.4 and 6.5 of Zimmerman & Núñez-Antón (2009).

Usage

test_order_gau(
  y,
  p = 0L,
  q = 1L,
  mu = NULL,
  use_modified = TRUE,
  order_null = NULL,
  order_alt = NULL
)

Arguments

Details

The test is based on the intervenor-adjusted sample partial correlations. Under the null hypothesis AD(p), the partial correlations r_(i,i-k|(i-k+1:i-1)) should be zero for k > p.

The likelihood ratio test statistic (Theorem 6.4) is:

-N \sum_{j=1}^{q} \sum_{i=p+j+1}^{n} \log(1 - r^2_{i,i-p-j\cdot(i-p-j+1:i-1)})

which is asymptotically chi-square with (2n - 2p - q - 1)(q/2) degrees of freedom.

The modified test (formula 6.9) adjusts for small-sample bias using Kenward's (1987) correction.

Value

A list with class gau_order_test containing:

method

Inference method used ("lrt").

p

Order under null hypothesis

q

Order increment

statistic

Test statistic value

statistic_modified

Modified test statistic (if use_modified = TRUE)

df

Degrees of freedom

p_value

P-value from chi-square distribution

p_value_modified

P-value from modified test (if use_modified = TRUE)

n_subjects

Number of subjects

n_time

Number of time points

References

Zimmerman, D.L. and Núñez-Antón, V. (2009). Antedependence Models for Longitudinal Data. Chapman & Hall/CRC. Chapter 6.

Kenward, M.G. (1987). A method for comparing profiles of repeated measurements. Applied Statistics, 36, 296-308.

See Also

test_one_sample_gau, test_homogeneity_gau

Examples

# Simulate AD(1) data
y <- simulate_gau(n_subjects = 50, n_time = 6, order = 1, phi = 0.5)

# Test AD(0) vs AD(1)
test01 <- test_order_gau(y, p = 0, q = 1)
print(test01)

# Test AD(1) vs AD(2)
test12 <- test_order_gau(y, p = 1, q = 1)
print(test12)

Likelihood Ratio Test for Antedependence Order (INAD Data)

Description

Compares INAD models of consecutive orders by their likelihood ratio statistic and BIC.

Usage

test_order_inad(
  y,
  order_null = 0,
  order_alt = 1,
  thinning = "binom",
  innovation = "pois",
  blocks = NULL,
  reference = c("none", "chibar"),
  weights = NULL,
  fit_null = NULL,
  fit_alt = NULL,
  use_chibar = NULL,
  ...
)

Arguments

Details

The statistic is \lambda = 2(\ell_{alt} - \ell_{null}), where \ell_{null} and \ell_{alt} are the maximized log-likelihoods. Under the null hypothesis the df added thinning coefficients are zero, on the boundary of their parameter space, so the ordinary chi-square reference is not valid. The limiting distribution is a chi-bar-square mixture whose weights depend on the efficient information of the added parameters; the binomial weights \binom{r}{j}2^{-r} apply only when these are orthogonal with block-diagonal efficient information. For models with block effects (INADFE) the shared effects couple the time points and these weights have not been established; requesting reference = "chibar" then gives a warning. Calibrated p-values in general require further work, for example a parametric bootstrap, which is not implemented here. BIC comparisons are always reported.

The fits are computed with fit_inad unless supplied. For order 2, fits from em_inad with newton_polish = TRUE are usually more precise and can be passed through fit_null and fit_alt.

Missing-data inputs are supported through the same na_action options available in fit_inad. If y has missing values and na_action is not supplied via ..., this function defaults to na_action = "marginalize".

Value

A list with class c("test_order_inad", "htest") containing:

statistic, lrt_stat

Likelihood ratio statistic

df

Number of added thinning parameters

p.value, p_value

Working chi-bar-square p-value when reference = "chibar"; NA otherwise

p_value_chibar

Same as p_value (kept for compatibility)

reference

The reference used; reference_note describes its assumptions

fit_null, fit_alt

The two fits

bic_null, bic_alt, bic_selected

BIC values and the order BIC prefers

table

Two-row model comparison table

settings

Input and derived settings for the test

References

Self, S.G. and Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. Journal of the American Statistical Association, 82, 605-610.

See Also

fit_inad, em_inad, bic_order_inad, test_stationarity_inad

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 15,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
out <- test_order_inad(
  y,
  order_null = 0,
  order_alt = 1,
  thinning = "binom",
  innovation = "pois",
  max_iter = 8
)
out$statistic
out$bic_selected
# Working chi-bar-square p-value for a baseline model:
test_order_inad(y, order_null = 0, order_alt = 1, reference = "chibar",
                fit_null = out$fit_null, fit_alt = out$fit_alt)$p_value

Likelihood Ratio Test for Stationarity (Categorical AD Data)

Description

Tests whether a categorical antedependence process satisfies stationarity constraints in the AD parameterization.

Usage

test_stationarity_cat(
  y,
  order = 1,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  test = c("lrt", "score", "mlrt")
)

Arguments

Details

The tested constraints are:

1. The marginal distribution P(Yk) is constant for all k

2. The transition probabilities P(Yk | Y(k-p), ..., Y(k-1)) are constant for all k > p

For AD order 1, these two constraints correspond to strict stationarity. For AD order greater than 1, this function should be interpreted as testing marginal-constancy plus time-invariant transitions; these constraints are not, in general, sufficient for full strict stationarity.

This is stronger than time-invariance alone, which only requires condition 2.

This function currently supports complete data only.

The null hypothesis is tested against the general (non-stationary) AD(p) model. The degrees of freedom are computed from the fitted parameter counts:

df = n_{params}(H_1) - n_{params}(H_0)

where H_1 is the unconstrained non-stationary model and H_0 is the stationary model.

Value

A list of class "cat_lrt" containing:

method

Inference method used: one of "lrt", "score", or "mlrt".

lrt_stat

Likelihood ratio test statistic

df

Degrees of freedom

p_value

P-value from chi-square distribution

fit_null

Fitted stationary model (H0)

fit_alt

Fitted non-stationary model (H1)

table

Summary data frame

References

Xie, Y. and Zimmerman, D. L. (2013). Antedependence models for nonstationary categorical longitudinal data with ignorable missingness: likelihood-based inference. Statistics in Medicine, 32, 3274-3289.

See Also

test_timeinvariance_cat, test_order_cat

Examples

# Simulate stationary AD(1) data
set.seed(123)
y <- simulate_cat(200, 6, order = 1, n_categories = 2)

# Test stationarity
test <- test_stationarity_cat(y, order = 1)
print(test)

Likelihood Ratio Test for Stationarity (Gaussian AD Data)

Description

Tests whether time-varying Gaussian AD covariance parameters can be constrained to be constant over time.

Usage

test_stationarity_gau(
  y,
  order = 1L,
  blocks = NULL,
  constrain = "both",
  fit_unconstrained = NULL,
  verbose = FALSE,
  max_iter = 2000L,
  rel_tol = 1e-08,
  ...
)

Arguments

Details

The mean structure is kept unrestricted in both models (time-specific means plus optional block shifts), and the test constrains covariance parameters: innovation standard deviations and/or antedependence coefficients.

The likelihood-ratio statistic is:

\lambda = 2(\ell_{alt} - \ell_{null})

where \ell_{null} and \ell_{alt} are maximized log-likelihoods under the constrained and unconstrained models, respectively.

Degrees of freedom are computed from the number of constraints implied by constrain.

Value

A list with class "test_stationarity_gau" containing:

method

Inference method used ("lrt").

fit_unconstrained

Unconstrained Gaussian AD fit

fit_constrained

Constrained Gaussian AD fit

constraint

Human-readable null constraint description

lrt_stat

Likelihood-ratio statistic

df

Degrees of freedom

p_value

Chi-square p-value

bic_unconstrained

BIC of unconstrained model

bic_constrained

BIC of constrained model

bic_selected

Model selected by BIC

table

Two-row model summary table

References

Zimmerman, D.L. and Nunez-Anton, V. (2009). Antedependence Models for Longitudinal Data. Chapman & Hall/CRC. Chapter 6.

See Also

run_stationarity_tests_gau, test_order_gau, fit_gau

Examples

set.seed(1)
y <- simulate_gau(n_subjects = 80, n_time = 6, order = 1, phi = 0.4, sigma = 1)

# Test jointly constant phi and sigma (order 1)
out <- test_stationarity_gau(y, order = 1, constrain = "both")
out$p_value

Likelihood Ratio Test for Time-Invariance (INAD Data)

Description

Tests whether time-varying INAD parameters can be constrained to be constant over time.

Usage

test_stationarity_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  blocks = NULL,
  constrain = "both",
  fit_unconstrained = NULL,
  verbose = FALSE,
  ...
)

Arguments

Details

The constrained model is refitted by maximizing the observed-data log-likelihood jointly over every parameter that remains free under the null hypothesis, including the negative binomial innovation sizes when they are not themselves constrained. Fitting alternates blockwise BOBYQA updates with a joint BOBYQA refinement until the log-likelihood stabilizes.

For order 2, "alpha1" constrains the lag-1 coefficients \alpha_{i,i-1}, i = 2, \ldots, n, to a common value, and "alpha2" constrains the lag-2 coefficients \alpha_{i,i-2}, i = 3, \ldots, n.

Degrees of freedom are the number of equality constraints imposed under the null model relative to the unconstrained model. Because the constrained parameters are interior under the null, the usual chi-square reference applies.

Missing-data inputs are supported through the same na_action options available in fit_inad. If y has missing values and na_action is not supplied via ..., this function defaults to na_action = "marginalize".

Value

A list with class "test_stationarity_inad" containing:

method

Inference method used ("lrt").

fit_unconstrained

Unconstrained INAD fit

fit_constrained

Constrained INAD fit, an inad_fit object

constraint

Human-readable null constraint description

statistic

Test statistic value

lrt_stat

Likelihood ratio test statistic

df

Degrees of freedom

p_value

Chi-square p-value

bic_unconstrained

BIC of unconstrained model

bic_constrained

BIC of constrained model

bic_selected

Model selected by BIC

table

Two-row model summary table

References

Li, C. and Zimmerman, D.L. (2026). Integer-valued antedependence models for longitudinal count data. Manuscript under review.

See Also

run_stationarity_tests_inad, test_order_inad, fit_inad

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 15,
  n_time = 5,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.3,
  theta = 2
)
out <- test_stationarity_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  constrain = "both",
  max_iter = 8
)
out$p_value

Likelihood Ratio Test for Time-Invariance (Categorical Data)

Description

Tests whether transition probabilities are constant over time in a categorical antedependence model.

Usage

test_timeinvariance_cat(
  y,
  order = 1,
  blocks = NULL,
  homogeneous = TRUE,
  n_categories = NULL,
  test = c("lrt", "score", "mlrt")
)

Arguments