The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Introduction to antedep

Overview

antedep supports three model families:

Packaged datasets include bolus_inad, cattle_growth, cochlear_implant, labor_force_cat, and race_100km.

data("race_100km")
dim(race_100km$y)
colMeans(race_100km$y)

Production-Readiness Matrix

Model Data type Complete-data fit/logLik Missing-data fit/logLik Notes
AD Continuous Ready Ready (fit_gau, logL_gau) Missing-data fit supports EM / observed-data likelihood modes
INAD Counts Ready Ready (fit_inad, logL_inad) Missing-data fit supports na_action = "marginalize"
CAT Categorical states Ready Ready (fit_cat, logL_cat) Missing-data fit supports orders 0, 1, 2

Visualizing Longitudinal Data

data("bolus_inad")
plot_profile(
  bolus_inad$y,
  blocks = bolus_inad$blocks,
  block_labels = c("Group 1", "Group 2"),
  title = "Bolus Longitudinal Profiles"
)

# More detailed dependence diagnostic
plot_prism(bolus_inad$y)

AD Workflow (Continuous Outcomes)

Complete-data fit and order comparison

y_gau <- simulate_gau(n_subjects = 20, n_time = 4, order = 1, phi = 0.5)
fit_gau1 <- fit_gau(y_gau, order = 1)
fit_gau2 <- fit_gau(y_gau, order = 2)

logLik(fit_gau1)
#> 'log Lik.' -105.672 (df=7)
BIC(fit_gau1, fit_gau2)
#>          df      BIC
#> fit_gau1  7 232.3141
#> fit_gau2  9 235.8026
bic_order_gau(y_gau, max_order = 2)
#> BIC-based Order Selection for Gaussian AD
#> ==========================================
#> 
#>  order     log_l      bic selected
#>      0 -112.6549 249.2756    FALSE
#>      1 -105.6720 244.2971     TRUE
#>      2 -104.4205 247.7855    FALSE
#> 
#> Best order by BIC: 1
# Example AD order test
test_order_gau(y_gau, p = 0, q = 1)$p_value
#> [1] 0.002952214

Missing-data AD fit

y_gau_miss <- y_gau
y_gau_miss[sample(length(y_gau_miss), 8)] <- NA

fit_gau_em <- fit_gau(y_gau_miss, order = 1, na_action = "em", em_max_iter = 10)
fit_gau_cc <- fit_gau(y_gau_miss, order = 1, na_action = "complete")
#> Warning: Only 13/20 subjects have complete data (65%).
#> Consider using na_action = 'em' for less information loss.

c(em_logLik = fit_gau_em$log_l, complete_case_logLik = fit_gau_cc$log_l)
#>            em_logLik complete_case_logLik 
#>            -94.00107            -63.11882
summary(fit_gau_em)
#> Gaussian Antedependence Model -- Summary
#> ========================================
#> 
#>   Order            : 1
#>   Subjects (n)     : 20
#>   Time points (T)  : 4
#>   Parameters       : 7
#>   Missing values   : 8 (10.0%)
#> 
#> Fit statistics:
#>   Log-likelihood   : -94.0011
#>   AIC              : 202.0021
#>   BIC              : 208.9723
#> 
#> Estimated mean (mu):
#> [1]  0.0961  0.1565  0.1244 -0.1773
#> 
#> Antedependence parameters (phi):
#> [1] 0.6772 0.1696 0.4342
#> 
#> Innovation std dev (sigma):
#> [1] 0.9905 0.6902 1.0489 0.8335

AD mean and covariance structure tests (complete data)

These AD inference tests currently require complete data.

# One-sample and two-sample mean-profile tests
mu0 <- colMeans(y_gau)
blocks_gau <- rep(1:2, each = nrow(y_gau) / 2)

test_one_sample_gau(y_gau, mu0 = mu0, p = 1)$p_value
#> [1] 1
test_two_sample_gau(y_gau, blocks = blocks_gau, p = 1)$p_value
#> [1] 0.01854549
# Wald contrast test: H0 that adjacent means are equal
C <- diff(diag(ncol(y_gau)))
test_contrast_gau(y_gau, C = C, p = 1)$p_value
#> [1] 0.8183341
# Covariance homogeneity across groups
test_homogeneity_gau(y_gau, blocks = blocks_gau, p = 1)$p_value
#>        g1 
#> 0.3237498

INAD Workflow (Count Outcomes)

Complete-data INAD fit

y_inad <- simulate_inad(
  n_subjects = 20,
  n_time = 4,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  alpha = 0.4,
  theta = 3
)

fit_inad1 <- fit_inad(y_inad, order = 1, thinning = "binom", innovation = "pois")
logLik(fit_inad1)
#> 'log Lik.' -158.3541 (df=7)
bic_order_inad(y_inad, max_order = 2, thinning = "binom", innovation = "pois")
#> 
#> BIC Model Order Comparison for INAD
#> ====================================
#> 
#>  order logLik n_params   BIC
#>      0 -164.7        4 341.3
#>      1 -158.4        7 337.7
#>      2 -158.4        9 343.7
#> 
#> Best order by BIC: 1

Missing-data INAD fit

y_inad_miss <- y_inad
y_inad_miss[sample(length(y_inad_miss), 10)] <- NA

fit_inad_miss <- fit_inad(
  y_inad_miss,
  order = 1,
  thinning = "binom",
  innovation = "pois",
  na_action = "marginalize",
  max_iter = 10
)

summary(fit_inad_miss)
#> INAD Antedependence Model -- Summary
#> ====================================
#> 
#>   Order            : 1
#>   Thinning         : binom
#>   Innovation       : pois
#>   Subjects (n)     : 20
#>   Time points (T)  : 4
#>   Parameters       : 7
#>   Missing values   : 10 (12.5%)
#> 
#> Fit statistics:
#>   Log-likelihood   : -142.2786
#>   AIC              : 298.5573
#>   BIC              : 305.5274
#> 
#> Estimated innovation mean (theta):
#> [1] 3.3775 2.0946 2.4217 3.2136
#> 
#> Estimated thinning parameter (alpha):
#> [1] 0.0000 0.7107 0.5782 0.1964
# Missing-data LRTs (observed-data likelihood)
blocks_inad <- rep(1:2, each = nrow(y_inad_miss) / 2)
# Order comparisons are boundary problems: compare BIC (no calibrated p-value)
ord_inad <- test_order_inad(y_inad_miss, order_null = 0, order_alt = 1, blocks = blocks_inad)
c(LR = ord_inad$lrt_stat, BIC_order0 = ord_inad$bic_null, BIC_order1 = ord_inad$bic_alt)
#>         LR BIC_order0 BIC_order1 
#>   6.743751 306.205451 308.448896
test_homogeneity_inad(y_inad_miss, blocks = blocks_inad, order = 1, test = "mean")$p_value
#> [1] 0.7852656

CAT Workflow (Categorical Outcomes)

Complete-data CAT fit and inference

y_cat <- simulate_cat(n_subjects = 20, n_time = 4, order = 1, n_categories = 3)
fit_cat1 <- fit_cat(y_cat, order = 1)
logLik(fit_cat1)
#> 'log Lik.' -79.08119 (df=20)

# Wald CIs (complete data)
head(confint(fit_cat1), 3)
#>                lower     upper
#> pi_t1_cat1 0.0602273 0.4397727
#> pi_t1_cat2 0.1852967 0.6147033
#> pi_t1_cat3 0.1409627 0.5590373
# Order tests for CAT
run_order_tests_cat(y_cat, max_order = 2)
#> Sequential Order Tests for Categorical AD
#> =========================================
#> 
#>      comparison order_null order_alt
#>  AD(0) vs AD(1)          0         1
#>  AD(1) vs AD(2)          1         2
#>                                                                   method
#>  Likelihood Ratio Test for Categorical AD Order (H0: AD(0) vs H1: AD(1))
#>  Likelihood Ratio Test for Categorical AD Order (H0: AD(1) vs H1: AD(2))
#>  lrt_stat df p_value significant
#>   13.0444 12  0.3658       FALSE
#>   37.2866 24  0.0410        TRUE
#> 
#> Selected order: 0

Missing-data CAT fit

y_cat_miss <- y_cat
y_cat_miss[sample(length(y_cat_miss), 12)] <- NA

fit_cat_miss <- fit_cat(y_cat_miss, order = 1, na_action = "marginalize")
logLik(fit_cat_miss)
#> 'log Lik.' -65.91133 (df=20)
fit_cat_miss$settings$na_action_effective
#> [1] "marginalize"
# EM alternative for CAT missing data (orders 0/1)
fit_cat_miss_em <- fit_cat(y_cat_miss, order = 1, na_action = "em", em_max_iter = 10)
#> Warning in em_cat(y = y, order = p, blocks = blocks_input, homogeneous =
#> homogeneous, : EM did not converge after 10 iterations (change = 0.000809).
logLik(fit_cat_miss_em)
#> 'log Lik.' -65.91256 (df=20)
# Missing-data order/homogeneity tests are supported
blocks_cat <- rep(1:2, each = nrow(y_cat_miss) / 2)
test_order_cat(y_cat_miss, order_null = 0, order_alt = 1)$p_value
#> [1] 0.2500927
test_homogeneity_cat(y_cat_miss, blocks = blocks_cat, order = 1)$p_value
#> [1] 0.3780152

Additional Inference Helpers (Complete Data)

# INAD homogeneity/stationarity tools
blocks <- rep(1:2, each = nrow(y_inad) / 2)
run_homogeneity_tests_inad(y_inad, blocks = blocks, order = 1)
run_stationarity_tests_inad(y_inad, order = 1, blocks = blocks)

# CAT stationarity/time-invariance tools
test_timeinvariance_cat(y_cat, order = 1)
test_stationarity_cat(y_cat, order = 1)

Known Limitations

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.