Learn R Programming

DMAR (version 1.0.0)

compare_cov_structures: Likelihood-Ratio Comparison of Covariance Structures

Description

Fits a long-format within-subjects regression under a menu of variance-covariance structures, from independence through the unstructured form, and returns a comparison table of log-likelihood, AIC, BIC, and pairwise likelihood-ratio tests against the most general structure (UN). Wraps gls.

Usage

compare_cov_structures(
  data,
  outcome,
  subject,
  time,
  fixed_effects = NULL,
  structures = c("IND", "CS", "CSH", "AR1", "ARH1", "TOEP", "TOEPH", "UN")
)

Value

A data.frame with one row per structure. Columns: structure, log_lik, AIC, BIC,

n_par, LRT_vs_UN_chisq, LRT_vs_UN_df,

LRT_vs_UN_p.

Arguments

data

Long-format data.frame with one row per subject-by-condition observation.

outcome

Character name of the response column.

subject

Character name of the subject-id column.

time

Character name of the time / within-subjects factor column.

fixed_effects

Right-hand-side formula for the fixed effects (default: ~ time).

structures

Character vector of structures to fit. Any subset of c("IND", "CS", "CSH", "AR1", "ARH1", "TOEP", "TOEPH", "UN") (default: all eight). Matching is case insensitive, so lowercase aliases such as "cs", "ar1", "csh", "arh1", "toep", and "un" are accepted and normalized to their canonical uppercase labels.

Author

Ken Kelley [email protected]

Details

Structures. Every structure below is nested in UN, so the likelihood-ratio test against UN is well defined for each.

  • IND: independent observations within subject (correlation = NULL in gls). Provided as a baseline.

  • CS: compound symmetry, a constant correlation and a single variance across time points: nlme::corCompSymm().

  • CSH: heterogeneous compound symmetry, a constant correlation with a separate variance at each time point: nlme::corCompSymm() with nlme::varIdent().

  • AR1: first-order autoregressive correlation with a single variance: nlme::corAR1().

  • ARH1: heterogeneous first-order autoregressive correlation with a separate variance at each time point: nlme::corAR1() with nlme::varIdent().

  • TOEP: Toeplitz (banded), a separate correlation at each lag with a single variance: nlme::corARMA() with autoregressive order one less than the number of time points and no moving-average term.

  • TOEPH: heterogeneous Toeplitz, the Toeplitz correlation with a separate variance at each time point: nlme::corARMA() with nlme::varIdent().

  • UN: unstructured, every variance and covariance free: nlme::corSymm() with nlme::varIdent().

LRT. Each restricted structure is compared against UN by the likelihood-ratio test. Both fits are re-estimated under ML (not REML) for the LRT, following nlme convention. The chi square statistic is \(-2 (\log L_{\mathrm{restricted}} - \log L_{\mathrm{UN}})\) on degrees of freedom equal to the difference in parameter count.

Caveats. The likelihood-ratio test against UN is valid because each listed structure is a restriction of UN. Two structures that are not nested in each other (for example CS and AR(1)) should be compared by AIC or BIC rather than by an LRT.

References

Littell, R. C., Milliken, G. A., Stroup, W. W., Wolfinger, R. D., & Schabenberger, O. (2006). SAS for mixed models (2nd ed.). SAS Institute.

Maxwell, S. E., Delaney, H. D., & Kelley, K. (2027). Designing experiments and analyzing data: A model comparison perspective (4th ed.). Routledge. (See Chapter 15.)

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

See Also

gls, logLik, anova

Other hypothesis tests: adjusted_means(), ancova(), anova_within(), ci_dunnett(), ci_scheffe(), ci_tukey_kramer(), contrast_test(), correlations_test(), equivalence_r(), equivalence_smd(), factorial_anova(), manova_split_plot(), mauchly_test(), mixed_anova(), obrien_test(), pairwise_within(), randomization_test(), randomization_test_paired(), regions_of_significance(), simple_effects_AB(), summary_t_test(), welch_t()

Examples

Run this code
# Four repeated measures on each of 30 subjects.
set.seed(113)
n <- 30; k <- 4
subj <- factor(rep(1:n, each = k))
tm   <- factor(rep(1:k, times = n))
y    <- as.vector(t(matrix(rnorm(n * k), n, k) +
                     rep(rnorm(n, 0, 1), each = k)))
d <- data.frame(y, subj, tm)

# All eight structures at once. Read the table by comparing AIC and BIC
# across rows, and use the likelihood-ratio test only for the nested
# comparison it reports, each structure against UN.
compare_cov_structures(d, outcome = "y", subject = "subj",
                       time = "tm")

# A subset, requested with lowercase aliases (matching is case
# insensitive).
compare_cov_structures(d, outcome = "y", subject = "subj",
                       time = "tm",
                       structures = c("cs", "csh", "ar1", "arh1"))

Run the code above in your browser using DataLab