---
title: "Validation against independent R implementations and PROC MIXED"
output: rmarkdown::html_vignette
bibliography: references.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{Validation against independent R implementations and PROC MIXED}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)

has_nlme <- requireNamespace("nlme", quietly = TRUE)
has_lme4 <- requireNamespace("lme4", quietly = TRUE)
has_lmer_test <- requireNamespace("lmerTest", quietly = TRUE)
has_mmrm <- requireNamespace("mmrm", quietly = TRUE)

official.sas.targets <- utils::read.csv(
  system.file(
    "validation",
    "sas",
    "sas-stat-14.3-targets.csv",
    package = "lmmix"
  ),
  na.strings = "",
  check.names = FALSE,
  stringsAsFactors = FALSE
)

official_sas_target <- function(model, quantity) {
  official.sas.targets[
    official.sas.targets$model == model &
      official.sas.targets$quantity == quantity,
    ,
    drop = FALSE
  ]
}

comparison_table <- function(term, lmmix, reference) {
  tibble::tibble(
    term = term,
    lmmix = as.numeric(lmmix),
    reference = as.numeric(reference),
    abs.diff = abs(as.numeric(lmmix) - as.numeric(reference))
  )
}

show_table <- function(x, digits = 6) {
  names(x) <- trimws(gsub("\\.+", " ", names(x)))
  knitr::kable(x, digits = digits, row.names = FALSE)
}
```

## Purpose and evidentiary standard

This vignette reports the numerical comparisons used to validate `lmmix`.
Every comparison is restricted to a model class for which the two
implementations use the same fixed-effects specification, covariance model,
and estimation criterion. A small difference is meaningful only after those
definitions have been aligned.

Validation is organized in four layers:

1. comparisons with independent R implementations for overlapping models;
2. structural checks for combined random-effect and residual-correlation
   models that no single reference R function covers;
3. comparisons with published output tables from four official PROC MIXED
   examples;
4. regression comparisons with stored PROC MIXED targets for one fully
   specified combined model.

The numerical tables below are produced by executable R code when the
vignette is built. The official PROC MIXED values are transcribed from the
published SAS/STAT output tables. The repository also stores a runnable SAS
program, but it does not contain a SAS execution transcript. These sections
are therefore documentation-based comparisons, not a new independent SAS
run.

## Comparison map

| lmmix model | reference function | quantities compared | tolerance |
|:--|:--|:--|:--|
| Random-intercept REML | `nlme::lme()` | fixed effects and log-likelihood: `1e-6`; variance components: `1e-4` | absolute |
| Random-slope REML | `nlme::lme()` | fixed effects and log-likelihood: `1e-5`; covariance parameters: `1e-4` | absolute |
| Random-intercept ML | `nlme::lme()` | fixed effects, log-likelihood | `1e-6` |
| Marginal CS and AR(1) REML | `nlme::gls()` | fixed effects, residual variance, correlation, log-likelihood | `1e-5` |
| Marginal AR(1) REML | `mmrm::mmrm()` | fixed effects and log-likelihood: `1e-5`; residual covariance matrix: `1e-4` | absolute |
| Random-intercept REML | `lmerTest::lmer()` | estimates and standard errors: `1e-5`; coefficient df: `5e-3` | absolute |
| Type III Satterthwaite tests | `lmerTest::lmer()` | numerator df: `1e-6`; denominator df and F statistics: `5e-3` | absolute |
| Marginal AR(1) Kenward-Roger | `mmrm::mmrm()` | fixed effects: `1e-5`; standard errors: `3e-3`; denominator df: `5e-3` | absolute |
| Crossed random intercepts | `lme4::lmer()` | log-likelihood: `1e-5`; covariance parameters: `1e-4` | absolute |
| Nested fixed-effect models | `lme4::anova()` | likelihood-ratio statistic and p-value: `1e-5` | absolute |
| Official split-plot example 79.1 | PROC MIXED published output | covariance parameters, REML criterion, type III tests | published precision |
| Official repeated-measures example 79.2 | PROC MIXED published output | UN and CS covariance, ML criterion, fixed effects, type III statistics | published precision |
| Official random-coefficients example 79.5 | PROC MIXED published output | random covariance, REML criterion, fixed effects, BLUPs, type III test | published precision |
| Official line-source example 79.6 | PROC MIXED published output | three random variances, `TOEP(4)` residual covariance, REML criterion | published precision |
| Random center plus AR(1) residual | PROC MIXED stored targets | covariance parameters, fixed effects, type III tests, LS-means, differences | target-specific precision |

`nlme` provides independent likelihood comparisons for random-effects and
structured marginal models [@pinheiro2000; @nlme2025]. `lmerTest` provides an
independent Satterthwaite implementation for models fitted by `lme4`
[@kuznetsova2017]. `mmrm` provides a second independent comparison for the
marginal repeated-measures model [@mmrm2026].

## Shared longitudinal data

The R comparisons use the `Orthodont` data distributed with `nlme`. `Occasion`
is created explicitly so that AR(1) visits have a declared order.

```{r orthodont-data, eval=has_nlme}
library(lmmix)

orthodont <- as.data.frame(nlme::Orthodont)
orthodont$Occasion <- factor(ave(
  seq_len(nrow(orthodont)),
  orthodont$Subject,
  FUN = seq_along
))

tibble::tibble(
  n.obs = nrow(orthodont),
  n.subjects = nlevels(orthodont$Subject),
  n.occasions = nlevels(orthodont$Occasion)
) |>
  show_table(digits = 0)
```

## Random-intercept REML comparison with `nlme::lme()`

The first comparison uses the same fixed formula, random intercept, and REML
criterion in both implementations.

```{r random-intercept-reml, eval=has_nlme}
fit.ri <- lmm(
  orthodont,
  distance ~ age + Sex,
  random = ~ 1 | Subject,
  method = "REML"
)

ref.ri <- nlme::lme(
  distance ~ age + Sex,
  random = ~ 1 | Subject,
  data = orthodont,
  method = "REML"
)

comparison_table(
  term = names(fixef(fit.ri)),
  lmmix = fixef(fit.ri),
  reference = nlme::fixef(ref.ri)
) |>
  show_table()
```

```{r random-intercept-reml-covariance, eval=has_nlme}
ref.ri.cov <- as.matrix(nlme::getVarCov(
  ref.ri,
  type = "random.effects"
))
ref.ri.variance <- c(
  ref.ri.cov[1L, 1L],
  ref.ri$sigma^2
)

comparison_table(
  term = c("random.intercept.var", "residual.var"),
  lmmix = VarCorr(fit.ri)$estimate,
  reference = ref.ri.variance
) |>
  show_table()

comparison_table(
  term = "logLik",
  lmmix = logLik(fit.ri),
  reference = logLik(ref.ri)
) |>
  show_table()
```

## Random-slope REML comparison with `nlme::lme()`

This comparison adds a correlated random slope for age. It checks two random
variances, their correlation, the residual variance, the fixed effects, and
the REML log-likelihood.

```{r random-slope-reml, eval=has_nlme}
fit.rs <- lmm(
  orthodont,
  distance ~ age + Sex,
  random = ~ 1 + age | Subject,
  method = "REML"
)

ref.rs <- nlme::lme(
  distance ~ age + Sex,
  random = ~ age | Subject,
  data = orthodont,
  method = "REML"
)

comparison_table(
  term = names(fixef(fit.rs)),
  lmmix = fixef(fit.rs),
  reference = nlme::fixef(ref.rs)
) |>
  show_table()
```

```{r random-slope-reml-covariance, eval=has_nlme}
ref.rs.cov <- as.matrix(nlme::getVarCov(
  ref.rs,
  type = "random.effects"
))
ref.rs.parameters <- c(
  ref.rs.cov[1L, 1L],
  ref.rs.cov[2L, 2L],
  ref.rs.cov[1L, 2L] / sqrt(ref.rs.cov[1L, 1L] * ref.rs.cov[2L, 2L]),
  ref.rs$sigma^2
)

comparison_table(
  term = c(
    "random.intercept.var",
    "random.age.var",
    "random.intercept.age.cor",
    "residual.var"
  ),
  lmmix = VarCorr(fit.rs)$estimate,
  reference = ref.rs.parameters
) |>
  show_table()

comparison_table(
  term = "logLik",
  lmmix = logLik(fit.rs),
  reference = logLik(ref.rs)
) |>
  show_table()
```

## Maximum-likelihood comparison with `nlme::lme()`

The ML comparison verifies that agreement is not specific to the REML
criterion.

```{r random-intercept-ml, eval=has_nlme}
fit.ml <- lmm(
  orthodont,
  distance ~ age + Sex,
  random = ~ 1 | Subject,
  method = "ML"
)

ref.ml <- nlme::lme(
  distance ~ age + Sex,
  random = ~ 1 | Subject,
  data = orthodont,
  method = "ML"
)

rbind(
  comparison_table(
    term = names(fixef(fit.ml)),
    lmmix = fixef(fit.ml),
    reference = nlme::fixef(ref.ml)
  ),
  comparison_table(
    term = "logLik",
    lmmix = logLik(fit.ml),
    reference = logLik(ref.ml)
  )
) |>
  show_table()
```

## Marginal CS and AR(1) comparisons with `nlme::gls()`

These models contain no random effect. The repeated-measures covariance is
compound symmetric or AR(1), and both implementations optimize a REML
criterion.

```{r marginal-nlme-fits, eval=has_nlme}
marginal.fits <- list()
marginal.refs <- list()

for (structure in c("cs", "ar1")) {
  marginal.fits[[structure]] <- lmm(
    orthodont,
    distance ~ age + Sex,
    repeated = ~ Occasion | Subject,
    structure = structure,
    method = "REML"
  )

  correlation <- if (structure == "cs") {
    nlme::corCompSymm(form = ~ 1 | Subject)
  } else {
    nlme::corAR1(form = ~ 1 | Subject)
  }

  marginal.refs[[structure]] <- nlme::gls(
    distance ~ age + Sex,
    data = orthodont,
    correlation = correlation,
    method = "REML"
  )
}
```

```{r marginal-nlme-fixed, eval=has_nlme}
marginal.fixed <- do.call(rbind, lapply(c("cs", "ar1"), function(structure) {
  fit <- marginal.fits[[structure]]
  ref <- marginal.refs[[structure]]
  table <- comparison_table(
    term = names(fixef(fit)),
    lmmix = fixef(fit),
    reference = stats::coef(ref)
  )
  tibble::add_column(table, structure = structure, .before = 1)
}))

show_table(marginal.fixed)
```

```{r marginal-nlme-covariance, eval=has_nlme}
marginal.covariance <- do.call(
  rbind,
  lapply(c("cs", "ar1"), function(structure) {
    fit <- marginal.fits[[structure]]
    ref <- marginal.refs[[structure]]
    reference.parameters <- c(
      ref$sigma^2,
      as.numeric(stats::coef(
        ref$modelStruct$corStruct,
        unconstrained = FALSE
      ))
    )
    table <- comparison_table(
      term = c("residual.var", "residual.cor"),
      lmmix = VarCorr(fit)$estimate,
      reference = reference.parameters
    )
    tibble::add_column(table, structure = structure, .before = 1)
  })
)

show_table(marginal.covariance)

tibble::tibble(
  structure = c("cs", "ar1"),
  lmmix = vapply(
    marginal.fits,
    function(x) as.numeric(logLik(x)),
    numeric(1)
  ),
  reference = vapply(
    marginal.refs,
    function(x) as.numeric(logLik(x)),
    numeric(1)
  )
) |>
  transform(abs.diff = abs(lmmix - reference)) |>
  show_table()
```

## Marginal AR(1) comparison with `mmrm::mmrm()`

The `mmrm` comparison provides a second independent implementation of the
marginal AR(1) model. It checks the complete residual covariance matrix rather
than only the variance and correlation parameters.

```{r marginal-mmrm, eval=has_nlme && has_mmrm}
ref.mmrm <- mmrm::mmrm(
  distance ~ age + Sex,
  data = orthodont,
  covariance = mmrm::cov_struct(
    type = "ar1",
    visits = "Occasion",
    subject = "Subject"
  ),
  reml = TRUE
)

fit.mmrm <- marginal.fits[["ar1"]]

comparison_table(
  term = names(fixef(fit.mmrm)),
  lmmix = fixef(fit.mmrm),
  reference = stats::coef(ref.mmrm)
) |>
  show_table()

tibble::tibble(
  quantity = c("logLik", "residual.covariance.max.abs.diff"),
  lmmix = c(as.numeric(logLik(fit.mmrm)), NA_real_),
  reference = c(as.numeric(logLik(ref.mmrm)), NA_real_),
  abs.diff = c(
    abs(as.numeric(logLik(fit.mmrm)) - as.numeric(logLik(ref.mmrm))),
    max(abs(
      unname(fit.mmrm$covariance$residual_base) -
        unname(mmrm::component(ref.mmrm, "varcor"))
    ))
  )
) |>
  show_table()
```

## Satterthwaite comparison with `lmerTest::lmer()`

The covariance model in this section is a random intercept with independent
residuals, which is supported by both packages. The coefficient comparison
checks estimates, standard errors, and one-dimensional Satterthwaite degrees
of freedom.

```{r lmer-test-coefficients, eval=has_nlme && has_lmer_test}
ref.lmer <- lmerTest::lmer(
  distance ~ age + Sex + (1 | Subject),
  data = orthodont,
  REML = TRUE
)

lmmix.coef <- generics::tidy(fit.ri)
lmer.coef <- summary(ref.lmer)$coefficients

tibble::tibble(
  term = lmmix.coef$term,
  lmmix.estimate = lmmix.coef$estimate,
  lmerTest.estimate = lmer.coef[, "Estimate"],
  estimate.abs.diff = abs(lmmix.estimate - lmerTest.estimate),
  lmmix.std.error = lmmix.coef$std.error,
  lmerTest.std.error = lmer.coef[, "Std. Error"],
  std.error.abs.diff = abs(lmmix.std.error - lmerTest.std.error),
  lmmix.df = lmmix.coef$df,
  lmerTest.df = lmer.coef[, "df"],
  df.abs.diff = abs(lmmix.df - lmerTest.df)
) |>
  show_table()
```

The type III comparison uses the same fixed-effect terms and the same
Satterthwaite denominator-df method.

```{r lmer-test-type3, eval=has_nlme && has_lmer_test}
lmmix.type3 <- anova(fit.ri)
lmer.type3 <- anova(ref.lmer, type = 3)

tibble::tibble(
  term = lmmix.type3$term,
  lmmix.num.df = lmmix.type3$num.df,
  lmerTest.num.df = lmer.type3[, "NumDF"],
  num.df.abs.diff = abs(lmmix.num.df - lmerTest.num.df),
  lmmix.den.df = lmmix.type3$den.df,
  lmerTest.den.df = lmer.type3[, "DenDF"],
  den.df.abs.diff = abs(lmmix.den.df - lmerTest.den.df),
  lmmix.statistic = lmmix.type3$statistic,
  lmerTest.statistic = lmer.type3[, "F value"],
  statistic.abs.diff = abs(lmmix.statistic - lmerTest.statistic)
) |>
  show_table()
```

## Kenward-Roger comparison with `mmrm::mmrm()`

The Kenward-Roger comparison uses a marginal AR(1) model, which both packages
fit by REML. It checks fixed effects, adjusted standard errors, and denominator
degrees of freedom.

```{r kenward-roger-mmrm, eval=has_nlme && has_mmrm}
fit.kr <- lmm(
  orthodont,
  distance ~ age + Sex,
  repeated = ~ Occasion | Subject,
  structure = "ar1",
  ddf = "kenward-roger"
)

mmrm.formula <- stats::as.formula(
  "distance ~ age + Sex + ar1(Occasion | Subject)",
  env = asNamespace("mmrm")
)
ref.kr <- mmrm::mmrm(
  mmrm.formula,
  data = orthodont,
  method = "Kenward-Roger"
)
fit.kr.table <- generics::tidy(fit.kr)
ref.kr.table <- coef(summary(ref.kr))

tibble::tibble(
  term = fit.kr.table$term,
  estimate.abs.diff = abs(fit.kr.table$estimate - coef(ref.kr)),
  std.error.abs.diff = abs(
    fit.kr.table$std.error - ref.kr.table[, "Std. Error"]
  ),
  df.abs.diff = abs(fit.kr.table$df - ref.kr.table[, "df"])
) |>
  show_table()
```

## Crossed random effects and model comparison

The crossed-intercept comparison uses two independent grouping factors. The
likelihood-ratio comparison uses nested fixed-effect models and the same ML
covariance specification in both packages.

```{r crossed-and-lrt, eval=has_nlme && has_lme4}
set.seed(42)
crossed <- expand.grid(
  site = factor(seq_len(8)),
  observer = factor(seq_len(6)),
  replicate = seq_len(3)
)
site.effect <- stats::rnorm(8, sd = 0.8)
observer.effect <- stats::rnorm(6, sd = 0.5)
crossed$response <- 2 +
  site.effect[crossed$site] +
  observer.effect[crossed$observer] +
  stats::rnorm(nrow(crossed), sd = 0.4)

fit.crossed <- lmm(
  crossed,
  response ~ 1,
  random = list(site = ~ 1 | site, observer = ~ 1 | observer)
)
ref.crossed <- lme4::lmer(
  response ~ 1 + (1 | site) + (1 | observer),
  data = crossed,
  REML = TRUE
)

tibble::tibble(
  quantity = c(
    "logLik", "site variance", "observer variance", "residual variance"
  ),
  lmmix = c(as.numeric(logLik(fit.crossed)), VarCorr(fit.crossed)$estimate),
  reference = c(
    as.numeric(logLik(ref.crossed)),
    as.data.frame(lme4::VarCorr(ref.crossed))$vcov
  ),
  abs.diff = abs(lmmix - reference)
) |>
  show_table()

fit.reduced <- lmm(
  orthodont,
  distance ~ age,
  random = ~ 1 | Subject
)
fit.full <- lmm(
  orthodont,
  distance ~ age + Sex,
  random = ~ 1 | Subject
)
lmmix.lrt <- suppressMessages(anova(fit.reduced, fit.full))
ref.reduced <- lme4::lmer(
  distance ~ age + (1 | Subject),
  data = orthodont,
  REML = FALSE
)
ref.full <- lme4::lmer(
  distance ~ age + Sex + (1 | Subject),
  data = orthodont,
  REML = FALSE
)
reference.lrt <- anova(ref.reduced, ref.full, refit = FALSE)

tibble::tibble(
  quantity = c("likelihood-ratio statistic", "p-value"),
  lmmix = c(lmmix.lrt$Chisq[2], lmmix.lrt$p.value[2]),
  reference = c(reference.lrt$Chisq[2], reference.lrt$`Pr(>Chisq)`[2]),
  abs.diff = abs(lmmix - reference)
) |>
  show_table()
```

## Combined covariance structures

The included `multicentre` data reproduce the multilocation repeated-measures
example in Section 28.3 of @milliken2009. They contain 153 rows, including 28
missing responses. The complete-case models therefore use 125 observations.

No single external R function used above fits all of the following components
in one object: a random center effect, a structured within-subject residual
covariance, and the `lmmix` Satterthwaite inference layer. For this combined
case, all five residual structures are checked for optimizer convergence, a
positive-definite likelihood Hessian, and a positive-definite marginal
covariance matrix.

```{r combined-structures}
combined.fits <- lapply(
  c("id", "cs", "ar1", "toep", "un"),
  function(structure) {
    lmm(
      multicentre,
      Y ~ Drug * Time,
      random = ~ 1 | Center,
      repeated = ~ Time | Center:Drug:Subject,
      structure = structure,
      method = "REML",
      ddf = "satterthwaite"
    )
  }
)
names(combined.fits) <- c("id", "cs", "ar1", "toep", "un")

tibble::tibble(
  structure = names(combined.fits),
  convergence = vapply(
    combined.fits,
    function(x) x$convergence$code,
    integer(1)
  ),
  hessian.positive.definite = vapply(
    combined.fits,
    function(x) x$convergence$hessian_positive_definite,
    logical(1)
  ),
  min.marginal.eigenvalue = vapply(
    combined.fits,
    function(x) {
      min(eigen(
        x$covariance$v,
        symmetric = TRUE,
        only.values = TRUE
      )$values)
    },
    numeric(1)
  )
) |>
  show_table()
```

## Official PROC MIXED examples

The official SAS/STAT 14.3 documentation supplies complete data, model
statements, and published output tables for several PROC MIXED examples
[@sas2017]. Four examples match the current `lmmix` model space exactly. The
package includes their data, a runnable SAS program, and a machine-readable
CSV of the published numerical targets.

| SAS example | `lmmix` coverage | status |
|:--|:--|:--|
| 79.1 Split-Plot Design | two independent random-intercept terms, REML, type III tests | exact overlap |
| 79.2 Repeated Measures | ML with UN and CS residual covariance | exact likelihood and covariance overlap |
| 79.5 Random Coefficients | correlated random intercept and slope, REML, BLUPs | exact overlap |
| 79.6 Line-Source Sprinkler | three random terms and banded `TOEP(4)` residual covariance | exact covariance and likelihood overlap |
| 79.4 Known G and R | externally fixed multivariate G and R matrices | outside the current API |
| 79.7 and 79.8 Influence | deletion diagnostics | influence diagnostics are not implemented |

The source tables use SAS's displayed precision. Comparisons therefore use
the number of decimals printed in the documentation. The executable files are
available through:

```{r official-sas-files}
system.file(
  "validation",
  "sas",
  package = "lmmix"
)
```

### Example 79.1: split-plot design

The split-plot example has random block and block-by-A intercepts. The
Satterthwaite denominator degrees of freedom coincide with the containment
values published by SAS for this balanced design.

```{r official-sas-split-fit}
official.split.fit <- lmm(
  sas_split_plot,
  Y ~ A * B,
  random = list(
    Block = ~ 1 | Block,
    A_Block = ~ 1 | Block:A
  ),
  method = "REML",
  ddf = "satterthwaite"
)

split.cov.target <- official_sas_target("split_plot", "covariance")
split.cov.fit <- VarCorr(official.split.fit)$estimate

tibble::tibble(
  quantity = split.cov.target$term,
  lmmix = split.cov.fit,
  SAS.target = split.cov.target$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

tibble::tibble(
  quantity = "-2 restricted log likelihood",
  lmmix = deviance(official.split.fit),
  SAS.target = official_sas_target(
    "split_plot",
    "deviance"
  )$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()
```

```{r official-sas-split-tests}
split.test.target <- official_sas_target("split_plot", "type3")
split.test.fit <- anova(official.split.fit)

tibble::tibble(
  term = split.test.fit$term,
  lmmix.den.df = split.test.fit$den.df,
  SAS.den.df = split.test.target$den_df,
  lmmix.F = split.test.fit$statistic,
  SAS.F = split.test.target$statistic,
  lmmix.p.value = split.test.fit$p.value,
  SAS.p.value = split.test.target$p_value
) |>
  show_table()
```

### Example 79.2: repeated measures

The Potthoff-Roy growth data are fitted by ML with the same fixed effects and
the same subject-level covariance blocks used in the SAS example. First, the
unstructured model reproduces all ten entries of the published covariance
matrix.

```{r official-sas-growth-un}
official.growth.un.fit <- lmm(
  sas_growth,
  y ~ Gender * Age,
  repeated = ~ Age | Person,
  structure = "un",
  method = "ML"
)

growth.un.target <- official_sas_target("growth_un", "covariance")
growth.un.matrix <- official.growth.un.fit$covariance$residual_base
growth.un.fit <- c(
  growth.un.matrix[1, 1],
  growth.un.matrix[2, 1], growth.un.matrix[2, 2],
  growth.un.matrix[3, 1], growth.un.matrix[3, 2],
  growth.un.matrix[3, 3],
  growth.un.matrix[4, 1], growth.un.matrix[4, 2],
  growth.un.matrix[4, 3], growth.un.matrix[4, 4]
)

tibble::tibble(
  covariance.parameter = growth.un.target$term,
  lmmix = growth.un.fit,
  SAS.target = growth.un.target$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

tibble::tibble(
  quantity = "-2 log likelihood",
  lmmix = deviance(official.growth.un.fit),
  SAS.target = official_sas_target(
    "growth_un",
    "deviance"
  )$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()
```

SAS uses a last-level reference parameterization, whereas R uses a first-level
treatment contrast by default. The following estimable transformation puts the
`lmmix` coefficients on the published SAS scale before comparison.

```{r official-sas-growth-fixed}
growth.sas.transformation <- rbind(
  Intercept = c(1, 1, 0, 0),
  "Gender F" = c(0, -1, 0, 0),
  Age = c(0, 0, 1, 1),
  "Age*Gender F" = c(0, 0, 0, -1)
)

growth_sas_solution <- function(fit) {
  variance <- growth.sas.transformation %*%
    vcov(fit) %*%
    t(growth.sas.transformation)
  tibble::tibble(
    term = rownames(growth.sas.transformation),
    estimate = as.vector(growth.sas.transformation %*% fixef(fit)),
    std.error = sqrt(diag(variance))
  )
}

growth.un.solution <- growth_sas_solution(official.growth.un.fit)
growth.un.fixed.target <- official_sas_target("growth_un", "fixed")

tibble::tibble(
  term = growth.un.solution$term,
  lmmix.estimate = growth.un.solution$estimate,
  SAS.estimate = growth.un.fixed.target$estimate,
  estimate.abs.diff = abs(lmmix.estimate - SAS.estimate),
  lmmix.std.error = growth.un.solution$std.error,
  SAS.std.error = growth.un.fixed.target$std_error,
  std.error.abs.diff = abs(lmmix.std.error - SAS.std.error)
) |>
  show_table()
```

The compound-symmetry model reproduces the common covariance, diagonal
enhancement, fixed-effects solution, and likelihood criterion.

```{r official-sas-growth-cs}
official.growth.cs.fit <- lmm(
  sas_growth,
  y ~ Gender * Age,
  repeated = ~ Age | Person,
  structure = "cs",
  method = "ML"
)

growth.cs.matrix <- official.growth.cs.fit$covariance$residual_base
growth.cs.fit <- c(
  growth.cs.matrix[1, 2],
  growth.cs.matrix[1, 1] - growth.cs.matrix[1, 2]
)
growth.cs.target <- official_sas_target("growth_cs", "covariance")

tibble::tibble(
  covariance.parameter = growth.cs.target$term,
  lmmix = growth.cs.fit,
  SAS.target = growth.cs.target$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

growth.cs.solution <- growth_sas_solution(official.growth.cs.fit)
growth.cs.fixed.target <- official_sas_target("growth_cs", "fixed")

tibble::tibble(
  term = growth.cs.solution$term,
  lmmix.estimate = growth.cs.solution$estimate,
  SAS.estimate = growth.cs.fixed.target$estimate,
  estimate.abs.diff = abs(lmmix.estimate - SAS.estimate),
  lmmix.std.error = growth.cs.solution$std.error,
  SAS.std.error = growth.cs.fixed.target$std_error,
  std.error.abs.diff = abs(lmmix.std.error - SAS.std.error)
) |>
  show_table()
```

PROC MIXED uses the between-within denominator-df method in this example.
`lmmix` does not currently expose that method, so the validation compares the
likelihood, covariance, fixed-effect estimates, standard errors, and F
statistics, but does not claim equality of denominator degrees of freedom.

```{r official-sas-growth-tests}
growth.type3 <- rbind(
  cbind(
    model = "UN",
    as.data.frame(anova(official.growth.un.fit))
  ),
  cbind(
    model = "CS",
    as.data.frame(anova(official.growth.cs.fit))
  )
)
growth.type3.target <- rbind(
  official_sas_target("growth_un", "type3"),
  official_sas_target("growth_cs", "type3")
)

tibble::tibble(
  model = growth.type3$model,
  term = growth.type3$term,
  lmmix.F = growth.type3$statistic,
  SAS.F = growth.type3.target$statistic,
  abs.diff = abs(lmmix.F - SAS.F)
) |>
  show_table()
```

### Example 79.5: random coefficients

The pharmaceutical stability example fits a correlated random intercept and
slope for each batch. It validates the random covariance matrix, residual
variance, fixed effects, BLUPs, REML criterion, and the Month test.

```{r official-sas-random-coefficients}
official.random.fit <- lmm(
  sas_random_coefficients,
  Y ~ Month,
  random = ~ 1 + Month | Batch,
  method = "REML",
  ddf = "satterthwaite"
)

random.cov.target <- official_sas_target(
  "random_coefficients",
  "covariance"
)
random.cov.fit <- c(
  official.random.fit$covariance$g$Batch[1, 1],
  official.random.fit$covariance$g$Batch[2, 1],
  official.random.fit$covariance$g$Batch[2, 2],
  official.random.fit$covariance$residual_base[1, 1]
)

tibble::tibble(
  covariance.parameter = random.cov.target$term,
  lmmix = random.cov.fit,
  SAS.target = random.cov.target$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

random.fixed.target <- official_sas_target(
  "random_coefficients",
  "fixed"
)
tibble::tibble(
  term = names(fixef(official.random.fit)),
  lmmix.estimate = fixef(official.random.fit),
  SAS.estimate = random.fixed.target$estimate,
  lmmix.std.error = sqrt(diag(vcov(official.random.fit))),
  SAS.std.error = random.fixed.target$std_error
) |>
  show_table()

random.blup.target <- official_sas_target(
  "random_coefficients",
  "random"
)
random.blup <- ranef(official.random.fit)[[1L]]
tibble::tibble(
  term = random.blup.target$term,
  lmmix = as.vector(t(as.matrix(
    random.blup[c("(Intercept)", "Month")]
  ))),
  SAS.target = random.blup.target$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

tibble::tibble(
  quantity = c("-2 restricted log likelihood", "Month F statistic"),
  lmmix = c(
    deviance(official.random.fit),
    anova(official.random.fit)$statistic
  ),
  SAS.target = c(
    official_sas_target(
      "random_coefficients",
      "deviance"
    )$estimate,
    official_sas_target(
      "random_coefficients",
      "type3"
    )$statistic
  ),
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()
```

### Example 79.6: line-source sprinkler irrigation

The line-source example exercises the combined model directly. It contains
three independent random-intercept terms and a four-band Toeplitz residual
covariance within each block-by-cultivar series. In `structure = "toep(4)"`,
the main diagonal and the first three covariance bands are estimated, while
all longer-lag covariances are fixed to zero.

```{r official-sas-line-source}
official.line.fit <- lmm(
  sas_line_source,
  Y ~ (Cult + Dir + Irrig)^2,
  random = list(
    Block = ~ 1 | Block,
    Block.Dir = ~ 1 | Block:Dir,
    Block.Irrig = ~ 1 | Block:Irrig
  ),
  repeated = ~ Sbplt | Block:Cult,
  structure = "toep(4)",
  method = "REML",
  ddf = "residual"
)

line.target <- official_sas_target("line_source", "covariance")
line.residual <- official.line.fit$covariance$residual_base
line.fit <- c(
  official.line.fit$covariance$g$Block[1, 1],
  official.line.fit$covariance$g$Block.Dir[1, 1],
  official.line.fit$covariance$g$Block.Irrig[1, 1],
  line.residual[1, 2],
  line.residual[1, 3],
  line.residual[1, 4],
  line.residual[1, 1]
)

tibble::tibble(
  covariance.parameter = line.target$term,
  lmmix = line.fit,
  SAS.target = line.target$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

tibble::tibble(
  quantity = "-2 restricted log likelihood",
  lmmix = deviance(official.line.fit),
  SAS.target = official_sas_target(
    "line_source",
    "deviance"
  )$estimate,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()
```

The published fixed-effect tests use the containment denominator-df method.
Because `lmmix` does not implement containment degrees of freedom, this
comparison is deliberately restricted to the covariance parameters and the
restricted likelihood criterion.

## Stored PROC MIXED comparison

The stored targets correspond to the following specification. The RANDOM and
REPEATED statements represent distinct covariance sources, and
`DDFM=SATTERTHWAITE` requests approximate denominator degrees of freedom as
documented by @sas2017.

```sas
proc mixed data=multicentre method=reml;
  class Center Drug Subject Time;
  model Y = Drug Time Drug*Time / ddfm=satterth solution;
  random Center;
  repeated Time / type=ar(1) subject=Subject(Center*Drug);
  lsmeans Drug Time / diff;
run;
```

The exact `lmmix` fit used in the comparison is:

```{r sas-fit}
sas.fit <- combined.fits[["ar1"]]
sas.fit
```

### Covariance parameters and likelihood criterion

```{r sas-covariance}
sas.covariance.target <- c(5.1737, 10.6702, 0.9351)
lmmix.covariance <- VarCorr(sas.fit)

tibble::tibble(
  parameter = paste(
    lmmix.covariance$group,
    lmmix.covariance$term,
    lmmix.covariance$component,
    sep = "."
  ),
  lmmix = lmmix.covariance$estimate,
  SAS.target = sas.covariance.target,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()

tibble::tibble(
  quantity = "deviance",
  lmmix = deviance(sas.fit),
  SAS.target = 490.62570532,
  abs.diff = abs(lmmix - SAS.target)
) |>
  show_table()
```

### Fixed effects

```{r sas-fixed-effects}
sas.fixed.target <- c(
  12.0524590,
  4.7647059,
  3.9411765,
  0.9210892,
  2.3922237,
  -0.1427475,
  -0.4735052,
  0.8233152,
  0.3342713
)

sas.fixed.table <- comparison_table(
  term = names(fixef(sas.fit)),
  lmmix = fixef(sas.fit),
  reference = sas.fixed.target
)
names(sas.fixed.table) <- c("term", "lmmix", "SAS.target", "abs.diff")
show_table(sas.fixed.table)
```

### Type III tests

The stored type III targets contain F statistics rounded to two decimal places
and the interaction p-value. The other two SAS p-values were not available in
the stored target set and are therefore reported as `not.stored` rather than
reconstructed.

```{r sas-type3}
lmmix.sas.type3 <- anova(sas.fit)

tibble::tibble(
  term = lmmix.sas.type3$term,
  lmmix.num.df = lmmix.sas.type3$num.df,
  lmmix.den.df = lmmix.sas.type3$den.df,
  lmmix.statistic = lmmix.sas.type3$statistic,
  SAS.statistic = c(11.43, 59.27, 1.35),
  statistic.abs.diff = abs(lmmix.statistic - SAS.statistic),
  lmmix.p.value = lmmix.sas.type3$p.value,
  SAS.p.value = c("not.stored", "not.stored", "0.2597")
) |>
  show_table()
```

### LS-means

The following table compares every stored Drug and Time marginal mean,
standard error, and denominator df.

```{r sas-lsmeans}
drug.means <- lsmeans(sas.fit, ~Drug)
time.means <- lsmeans(sas.fit, ~Time)

lmmix.means <- rbind(
  tibble::tibble(
    factor = "Drug",
    level = as.character(drug.means$Drug),
    estimate = drug.means$estimate,
    std.error = drug.means$std.error,
    df = drug.means$df
  ),
  tibble::tibble(
    factor = "Time",
    level = as.character(time.means$Time),
    estimate = time.means$estimate,
    std.error = time.means$std.error,
    df = time.means$df
  )
)

sas.mean.estimate <- c(
  13.1568, 18.1484, 17.0516,
  14.9543, 15.6700, 17.7324
)
sas.mean.std.error <- c(
  1.5290, 1.5324, 1.5327,
  1.3961, 1.3985, 1.4032
)
sas.mean.df <- c(2.98, 3.01, 3.01, 2.08, 2.09, 2.12)

tibble::tibble(
  factor = lmmix.means$factor,
  level = lmmix.means$level,
  lmmix.estimate = lmmix.means$estimate,
  SAS.estimate = sas.mean.estimate,
  estimate.abs.diff = abs(lmmix.estimate - SAS.estimate),
  lmmix.std.error = lmmix.means$std.error,
  SAS.std.error = sas.mean.std.error,
  std.error.abs.diff = abs(lmmix.std.error - SAS.std.error),
  lmmix.df = lmmix.means$df,
  SAS.df = sas.mean.df
) |>
  show_table()
```

### Pairwise LS-mean differences

```{r sas-lsmeans-pairs}
drug.pairs <- lsmeans(sas.fit, pairwise ~ Drug)$contrasts
time.pairs <- lsmeans(sas.fit, pairwise ~ Time)$contrasts

lmmix.pairs <- rbind(
  tibble::add_column(drug.pairs, factor = "Drug", .before = 1),
  tibble::add_column(time.pairs, factor = "Time", .before = 1)
)

sas.pair.estimate <- c(
  -4.9916, -3.8948, 1.0968,
  -0.7157, -2.7781, -2.0624
)
sas.pair.std.error <- c(
  1.0990, 1.0987, 1.1043,
  0.1845, 0.2714, 0.2056
)
sas.pair.df <- c(46.1, 46.1, 46.9, 68.1, 73.2, 68.6)
sas.pair.p.value <- c(
  "0.00004", "0.0009", "0.3257",
  "0.0002", "<0.0001", "<0.0001"
)

tibble::tibble(
  factor = lmmix.pairs$factor,
  contrast = lmmix.pairs$contrast,
  lmmix.estimate = lmmix.pairs$estimate,
  SAS.estimate = sas.pair.estimate,
  estimate.abs.diff = abs(lmmix.estimate - SAS.estimate),
  lmmix.std.error = lmmix.pairs$std.error,
  SAS.std.error = sas.pair.std.error,
  std.error.abs.diff = abs(lmmix.std.error - SAS.std.error),
  lmmix.df = lmmix.pairs$df,
  SAS.df = sas.pair.df,
  lmmix.p.value = lmmix.pairs$p.value,
  SAS.p.value = sas.pair.p.value
) |>
  show_table()
```

Estimates, covariance parameters, and standard errors are tested with absolute
tolerance `1e-3`. Degrees of freedom are compared at the precision stored in
the targets. Type III F statistics are compared after rounding to two decimal
places. Inequalities such as `<0.0001` are checked as inequalities, not as
invented exact p-values.

## Automated regression tests

The numerical claims in this vignette are enforced in the package test suite.

| test file | enforced comparison |
|:--|:--|
| `test-likelihood.R` | `nlme::lme()`, `nlme::gls()`, and `mmrm::mmrm()` likelihood and covariance comparisons |
| `test-inference.R` | `lmerTest::lmer()` coefficient df and type III Satterthwaite comparisons |
| `test-validation-sas-official.R` | official SAS examples 79.1, 79.2, 79.5, and 79.6, including covariance, likelihood, fixed effects, BLUPs, and supported type III statistics |
| `test-version-030.R` | fixed-band Toeplitz, covariance confidence intervals, and parametric-bootstrap likelihood-ratio tests |
| `test-validation-sas.R` | every stored PROC MIXED covariance, fixed-effect, LS-mean, and difference target |
| `test-covariance.R` | positive definiteness and structure-specific covariance checks |
| `test-methods.R` | S3 outputs, `ggplot2` diagnostics, and readable printed table headings |

## Reproducibility information

The versions used to build this vignette are reported rather than assumed.

```{r validation-versions}
validation.packages <- c("lmmix", "nlme", "lmerTest", "mmrm")

tibble::tibble(
  package = validation.packages,
  version = vapply(
    validation.packages,
    function(package) {
      if (requireNamespace(package, quietly = TRUE)) {
        as.character(utils::packageVersion(package))
      } else {
        "not.installed"
      }
    },
    character(1)
  )
) |>
  show_table()

tibble::tibble(
  R.version.string = R.version.string,
  platform = R.version$platform
) |>
  show_table()
```

## Boundaries of the evidence

The validation covers Kenward-Roger inference for overlapping random-effects
and marginal models, crossed random intercepts, fixed-effect likelihood-ratio
comparisons, and one official combined random-effect and banded-residual model.
It does not validate generalized responses or large-scale sparse performance
because these are outside the model scope. For the official SAS examples,
values are transcribed from published output rather than generated by a SAS
session in this repository. Combined-model Kenward-Roger inference still has
structural tests only. An archived SAS log and output would be required to
upgrade that evidence.