Package 'antedep'

Description

Computes AIC using the fitted log likelihood and a parameter count for categorical antedependence parameters.

Usage

aic_cat(fit)

Arguments

Details

The AIC is computed as:

$AIC = -2 \times \ell + 2k$

where $\ell$ is the log-likelihood and $k$ is the number of free parameters.

Value

A numeric scalar AIC value.

Examples

set.seed(1)
y <- simulate_cat(40, 5, order = 1, n_categories = 2)
fit <- fit_cat(y, order = 1)
aic_cat(fit)

Description

Computes AIC using the fitted log likelihood and a parameter count that respects identifiability constraints for the Gaussian antedependence parameters.

Usage

aic_gau(fit)

Arguments

Details

The AIC is computed as:

$AIC = -2 \times \ell + 2k$

where $\ell$ is the log-likelihood and $k$ is the number of free parameters.

This function applies to Gaussian AD fits from fit_gau.

Value

A numeric scalar 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)
aic_gau(fit)

Description

Computes AIC using the fitted log likelihood and a parameter count that respects structural zeros and identifiability constraints.

Usage

aic_inad(fit)

Arguments

Details

The AIC is computed as:

$AIC = -2 \times \ell + 2k$

where $\ell$ is the log-likelihood and $k$ is the number of free parameters.

Value

A numeric scalar AIC value.

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 40,
  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)
aic_inad(fit)

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)

Description

Computes BIC using the fitted log likelihood and a parameter count for categorical antedependence parameters.

Usage

bic_cat(fit, n_subjects = NULL)

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.

Value

A numeric scalar BIC value.

Examples

set.seed(1)
y <- simulate_cat(40, 5, order = 1, n_categories = 2)

# Fit models of different orders
fit0 <- fit_cat(y, order = 0)
fit1 <- fit_cat(y, order = 1)
fit2 <- fit_cat(y, order = 2)

# Compare BIC
c(BIC_0 = bic_cat(fit0), BIC_1 = bic_cat(fit1), BIC_2 = bic_cat(fit2))

Description

Computes BIC using the fitted log likelihood and a parameter count that respects identifiability constraints for the Gaussian antedependence parameters.

Usage

bic_gau(fit, n_subjects = NULL)

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.

This function applies to Gaussian AD fits from fit_gau. For categorical and INAD models, use bic_cat and bic_inad.

Value

A numeric scalar BIC 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))

Description

Computes BIC using the fitted log likelihood and a parameter count that respects structural zeros and identifiability constraints.

Usage

bic_inad(fit, n_subjects = NULL)

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.

Value

A numeric scalar BIC value.

Examples

set.seed(1)
y <- simulate_inad(
  n_subjects = 40,
  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)
bic_inad(fit, n_subjects = nrow(y))

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$table)
print(result$best_order)

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
ord$table

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

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.

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

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")

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

Description

Computes confidence intervals for selected parameters from a fitted INAD model. For the fixed effect case, Wald intervals for time varying alpha and theta are computed via Louis identity for supported thinning-innovation combinations. For block effects tau, profile likelihood intervals are computed by fixing one component of tau and re maximizing the log likelihood over nuisance parameters. For negative binomial innovations, Wald intervals for the innovation size parameter are computed using a one dimensional observed information approximation per time point, holding other parameters fixed at their fitted values.

Usage

ci_inad(
  y,
  fit,
  blocks = NULL,
  level = 0.95,
  idx_time = NULL,
  ridge = 0,
  profile_maxeval = 2500,
  profile_xtol_rel = 1e-06
)

Arguments

Value

An object of class inad_ci, a list with elements settings, level, alpha, theta, nb_inno_size, and tau. Each non NULL interval element is a data frame with columns param, est, lower, upper, and possibly se and width.

Examples

data("bolus_inad", package = "antedep")
y <- bolus_inad$y
blocks <- bolus_inad$blocks
fit <- fit_inad(y, order = 1, thinning = "binom", innovation = "bell", blocks = blocks)
ci <- ci_inad(y, fit, blocks = blocks)
ci$alpha
ci$theta
ci$tau

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

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

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

Description

Fits INAD models using the Expectation-Maximization algorithm. This is an alternative to direct likelihood optimization.

Usage

em_inad(
  y,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  blocks = NULL,
  max_iter = 200,
  tol = 1e-07,
  alpha_init = NULL,
  theta_init = NULL,
  tau_init = NULL,
  nb_inno_size = NULL,
  safeguard = TRUE,
  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.

Value

A list with class "inad_fit" containing estimated parameters.

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

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

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

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")
)

Arguments

Details

No fixed effect: time separable optimization using logL_inad_i with theta eliminated by moment equations for order 1 and 2.

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

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

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.

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
)

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)

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
)

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
)

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

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)))

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")

Description

Print method for cat_ci objects

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for cat_fit objects

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for cat_lrt objects

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for BIC order selection

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for AD confidence intervals

Usage

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

Arguments

Value

The input object, invisibly.

Description

Print method for AD contrast test

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for gau_fit objects.

Usage

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

Arguments

Value

The input object, invisibly.

Description

Print method for AD homogeneity test

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for AD mean test

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for AD order test

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for homogeneity_tests_inad

Usage

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

Arguments

Value

Invisibly returns x.

Description

Print method for INAD confidence intervals

Usage

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

Arguments

Value

The input object, invisibly.

Description

Print method for INAD model fits

Usage

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

Arguments

Value

The input object, invisibly.

Description

Print method for test_homogeneity_inad

Usage

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

Arguments

Value

Invisibly returns x.

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.

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")
tests <- run_homogeneity_tests_inad(bolus_inad$y, bolus_inad$blocks,
                                     order = 1, thinning = "nbinom",
                                     innovation = "bell")
print(tests)

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$table)
cat("Selected order:", result$selected_order, "\n")

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$table)

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

Description

Run all stationarity-related tests for INAD

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 = 25,
  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 = 15
)
out$summary

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

Integer matrix with n_subjects rows and n_time columns, where each entry is a category code from 1 to c.

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)

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

numeric matrix with dimension n_subjects by n_time

Examples

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

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

integer matrix of counts with dimension n_subjects by n_time

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)

Description

Summary method for cat_ci objects

Usage

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

Arguments

Value

A data frame summarizing all CIs

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.

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.

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)

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)

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)

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. Biostatistics. Section 3.7.

See Also

fit_inad, test_order_inad, test_stationarity_inad

Examples

data("bolus_inad")
y <- bolus_inad$y
blocks <- bolus_inad$blocks

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

# Test only for mean differences (M1 vs M2)
test_mean <- test_homogeneity_inad(y, blocks, order = 1,
                                   thinning = "nbinom", innovation = "bell",
                                   test = "mean")
print(test_mean)

# Test for dependence differences given different means (M2 vs M3)
test_dep <- test_homogeneity_inad(y, blocks, order = 1,
                                  thinning = "nbinom", innovation = "bell",
                                  test = "dependence")
print(test_dep)

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)

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$table)

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

# 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)

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)

Description

Performs a likelihood ratio test comparing INAD models of different orders.

Usage

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

Arguments

Details

The test compares nested INAD models of orders order_null and order_alt = order_null + 1 using:

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

where $\ell_{null}$ and $\ell_{alt}$ are maximized log-likelihoods under the null and alternative models.