---
title: Inference Where Standard Mixed-Model p-values Break Down
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Inference Where Standard Mixed-Model p-values Break Down}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

```{r load-package}
library(mixeff)
```

Mixed models sometimes converge but still leave us in a difficult
inference situation. The fit exists, the estimates are usable, and
the model may be scientifically reasonable, but some familiar
p-values are no longer well justified.

This vignette focuses on one common case: boundary or singular fits.
In a mixed model, this can happen when a variance component is
estimated very close to zero, or when random-effect terms are
estimated as perfectly or nearly perfectly correlated. Equivalently,
the fitted random-effect covariance matrix has lower effective rank
than the model requested.

These are exactly the situations where ordinary large-sample tests
can become least trustworthy. Wald tests, Satterthwaite
approximations, and Kenward-Roger adjustments all rely on regularity
conditions that are strained or broken when a variance parameter is
on the boundary of its parameter space.

The statistical reason is well known. Self and Liang (1987) showed
that likelihood-ratio statistics do not necessarily follow their
usual chi-square reference distributions when a parameter lies on
the boundary. Stram and Lee (1994) specialized this issue to
variance components in linear mixed models. Kenward and Roger (1997)
derived their fixed-effect adjustment from a Taylor expansion of the
variance estimator; that derivation assumes a regular interior
solution for the variance parameters. At a boundary, those
derivatives may not exist or may not behave regularly.

`mixeff` handles this by making the inference contract explicit.
First, it does not silently present every number as equally
reliable. Asymptotic rows are labelled with a
[closed-enum reliability reason](inference-method-glossary.html),
such as `asymptotic_wald_z_fallback` or
`*_unavailable_at_boundary`. Second, when bootstrap inference is the
better-supported route, `mixeff` exposes it as a labelled method
rather than an informal workaround. Parametric bootstrap methods
have a long history in this setting; for example, Halekoh and
Højsgaard (2014) implemented parametric bootstrap tests in
`pbkrtest` for small-sample and boundary cases where Kenward-Roger
methods may need support.

The goal is not to say that a singular fit is automatically useless.
The goal is to keep the covariance state attached to the inference
result, so that available numbers and unavailable numbers are both
documented.

## What fit are we worried about?

This small repeated-measures example has subject-specific intercepts and slopes.
In this data set, those two random effects are nearly collinear, which is enough
to produce a reduced-rank random-effect covariance estimate.

```{r singular-data, include = FALSE}
set.seed(3)
n_subj <- 18L
days <- 0:9
b0 <- rnorm(n_subj, sd = 30)
b1 <- 0.5 * b0 / 30 * 10 + rnorm(n_subj, sd = 0.5)

sleep_like <- do.call(rbind, lapply(seq_len(n_subj), function(i) {
  data.frame(
    subj = factor(i),
    days = days,
    rt = 250 + b0[i] + (10 + b1[i]) * days +
      rnorm(length(days), sd = 20)
  )
}))
```

```{r fit-singular}
fit <- lmm(
  rt ~ days + (1 + days | subj),
  sleep_like,
  control = mm_control(verbose = -1)
)

fit_status(fit)
is_singular(fit)
```

```{r singular-checks, include = FALSE}
stopifnot(inherits(fit, "mm_lmm"))
stopifnot(isTRUE(is_singular(fit)))
stopifnot(identical(fit_status(fit), "converged_reduced_rank"))
```

The fit converged, but the covariance estimate is singular. That matters for
inference, so `mixeff` records it in the model diagnostics.

```{r covariance-diagnostics}
diagnostics(fit)$table[, c("code", "severity", "stage", "message")]
```

The requested random-effect covariance for `(1 + days | subj)` has rank 2: one
dimension for the subject intercept and one for the subject slope. The fitted
covariance has effective rank 1. In plain language, the model asked for two
random-effect dimensions, but the data support only one effective dimension.

That does not automatically invalidate the model. It does mean that inference
methods relying on a regular full-rank covariance estimate need to be treated
carefully.

## Which inference routes are available?

Before testing a term, ask which inference routes are defined for this fit.
`inference_options()` is an audit map, not a recommender. Each row describes one
method that `mixeff` knows about. The row says whether the method can run, why
it can or cannot run, and what the rough computational cost is.

The printed object is the reader-facing map. It keeps the raw enum
columns in `routes$table` for scripts, but explains the route status
with display columns. The complete contract vocabulary is in the
[inference method glossary](inference-method-glossary.html).

```{r route-map}
routes <- inference_options(fit, "days", nsim = 200)
routes
```

```{r route-checks, include = FALSE}
stopifnot("bootstrap" %in% routes$table$method)
satt <- routes$table[routes$table$method == "satterthwaite", ]
stopifnot(identical(satt$expected_status, "not_assessed"))
```

Read this table as follows. The Wald route can run immediately, but it is
labelled as a low-reliability asymptotic fallback. Satterthwaite and
Kenward-Roger are refused because the variance-parameter derivatives needed by
those approximations are not available at the boundary. The parametric bootstrap
route can run and is the documented fixed-effect testing route for this fit. The
bootstrap LRT route is refused here because this model was fit by REML, while
that route requires an ML fit. The cluster bootstrap route is refused for
fixed-effect p-values because its current target is an estimator distribution,
not a null hypothesis test. Profile confidence intervals are not certified for
this boundary fit.

```{r wald-status}
inference_table(fit)$table[, c("term", "method", "status",
                               "reliability", "reliability_reason")]
```

The important point is that `mixeff` distinguishes three cases:

1. A method can run and is considered available.
2. A method can produce a number, but the number is labelled as low reliability.
3. A method is refused because the method's assumptions or target are not
   defined for this fit.

That distinction is especially useful in scripts and reports, because
unavailable inference is recorded with stable reason codes rather than
disappearing as an error message.

## How do I test the term anyway?

For the fixed-effect term `days`, request the parametric bootstrap route
explicitly. Here the bootstrap target is a fixed-effect null model. `mixeff`
simulates data under the null hypothesis for the tested term, refits the model
to each simulated data set, and compares the observed test statistic with the
bootstrap reference distribution.

The small `nsim` value below keeps the vignette fast. For a real analysis, use
more bootstrap replicates.

```{r bootstrap-term}
term_boot <- test_effect(
  fit,
  "days",
  method = "bootstrap",
  bootstrap = bootstrap_control(nsim = 50, seed = 1)
)

term_boot$table[, c("term", "statistic_name", "p_value",
                    "method", "status", "reliability_reason")]
```

```{r bootstrap-term-checks, include = FALSE}
stopifnot(identical(term_boot$table$status, "available"))
stopifnot(identical(term_boot$table$method, "bootstrap"))
stopifnot(is.finite(term_boot$table$p_value))
```

The p-value is returned as part of the labelled result row. It is not
reconstructed by hand in R.

The result also carries a bootstrap payload in `details`, including the
requested number of replicates, the number of successful refits, the boundary
rate among bootstrap refits, and the Monte Carlo standard error.

```{r bootstrap-payload}
run <- term_boot$table$details[[1]]$bootstrap
data.frame(
  requested_replicates = run$requested_replicates,
  successful_replicates = run$successful_replicates,
  boundary_rate = round(run$boundary_rate, 3),
  mcse = round(run$mcse, 4)
)
```

The boundary rate is not a failure count. It reports how often the bootstrap
refits also ended on a covariance boundary. That is useful diagnostic
information: it tells you that the boundary behavior is not just a one-off
feature of the observed data set.

With `nsim = 50`, the bootstrap p-value is also necessarily coarse. For example,
a p-value of about 0.0196 corresponds to the smallest non-zero value possible
under a common plus-one bootstrap adjustment, `1 / (50 + 1)`. Use more
replicates when the p-value itself will be reported or compared with a decision
threshold.

## What about confidence intervals?

For fixed-effect confidence intervals, use `confint(method = "bootstrap")`.
This bootstrap has a different target from the fixed-effect test above. It
simulates from the full fitted model and summarizes the resulting estimator
distribution. That target is appropriate for percentile-style confidence
intervals and diagnostics, but it is not the same as a null distribution for a
fixed-effect p-value.

```{r bootstrap-ci}
ci <- confint(
  fit,
  parm = "days",
  method = "bootstrap",
  bootstrap = bootstrap_control(nsim = 50, seed = 2)
)

ci
```

```{r bootstrap-ci-checks, include = FALSE}
stopifnot(all(is.finite(ci)))
stopifnot(identical(attr(ci, "method"), "bootstrap_full_model_distribution"))
payload <- attr(ci, "bootstrap")[[1]]
stopifnot(identical(payload$metadata$target$kind, "full_model_distribution"))
stopifnot(length(payload$replicate_statistics) >= payload$metadata$successful_replicates,
          length(payload$replicate_statistics) <= payload$metadata$completed_replicates)
```

The attached bootstrap payload records what was simulated and why the interval
should not be reinterpreted as a fixed-effect hypothesis test.

```{r ci-payload}
payload <- attr(ci, "bootstrap")[[1]]
payload$metadata$target$kind
payload$metadata$notes
```

This distinction is central:

- `test_effect(..., method = "bootstrap")` uses a fixed-effect null target and
  returns a p-value.
- `confint(..., method = "bootstrap")` uses a full-model estimator-distribution
  target and returns an interval.
- The second result should not be used to reverse-engineer the first.

## Where does cluster bootstrap stand?

Cluster resampling is available as an estimator-distribution target. It
resamples grouping levels and summarizes how estimates vary across those
resamples.

In the current schema, that target does not certify fixed-effect p-values.
Therefore a request for a fixed-effect cluster-bootstrap p-value refuses cleanly
instead of inventing a null distribution.

```{r cluster-refusal}
cluster_row <- test_effect(
  fit,
  "days",
  method = "cluster_bootstrap",
  group = "subj"
)

cluster_row$table[, c("term", "method", "status", "p_value", "reason_code")]
```

```{r cluster-checks, include = FALSE}
stopifnot(identical(cluster_row$table$status, "not_assessed"))
stopifnot(is.na(cluster_row$table$p_value))
stopifnot(identical(cluster_row$table$reason_code,
                    "bootstrap_cluster_resample_p_value_unavailable"))
```

That refusal is part of the same inference contract as the available bootstrap
result. Available numbers are labelled with their method and reliability reason.
Unavailable numbers are labelled with stable reason codes explaining why they
are unavailable.

## Variance-component boundary tests

The boundary likelihood-ratio route is for random-effect variance components,
not ordinary fixed effects.

This is a different testing problem. For a variance component, the null
hypothesis can put the parameter exactly on the boundary of the parameter space,
because variances cannot be negative. In the simplest one-component case, the
reference distribution is the Self-Liang 50:50 mixture:
`0.5 * chi-square(0) + 0.5 * chi-square(1)`.

Fit a simpler random-intercept model by ML:

```{r boundary-lrt}
ri_fit <- lmm(
  rt ~ days + (1 | subj),
  sleep_like,
  REML = FALSE,
  control = mm_control(verbose = -1)
)

re_lrt <- test_random_effect(ri_fit, "subj", method = "boundary_lrt")
re_lrt
```

```{r boundary-lrt-table}
reporting_table(re_lrt)$table[, c("term", "statistic", "p_value",
                                  "reference_distribution", "status")]
```

```{r boundary-lrt-checks, include = FALSE}
stopifnot(inherits(re_lrt, "mm_random_effect_test"))
stopifnot(identical(re_lrt$table$status, "available"))
stopifnot(identical(re_lrt$table$reference_distribution,
                    "0.5 * chi-square(0) + 0.5 * chi-square(1)"))
stopifnot(is.finite(re_lrt$table$p_value))
```

The reporting table keeps the reference distribution attached to the result.

Asking for the same method on a fixed effect returns a typed refusal.

```{r boundary-fixed-refusal}
fixed_boundary <- test_effect(ri_fit, "days", method = "boundary_lrt")
fixed_boundary$table[, c("term", "method", "status", "reason_code")]
```

This is not a software gap. It is a methodological distinction. Boundary
likelihood-ratio theory applies naturally to variance-component tests, where
the null hypothesis places a variance on the edge of the parameter space. It is
not the reference distribution for an ordinary fixed-effect Wald row.

```{r boundary-fixed-checks, include = FALSE}
stopifnot(identical(fixed_boundary$table$status, "unsupported"))
stopifnot(identical(fixed_boundary$table$reason_code,
                    "boundary_lrt_not_applicable_to_fixed_effects"))
```

## Summary

A singular mixed-model fit can still be useful, but it changes the inference
problem.

For this fit:

- the model converges, but the random-effect covariance is reduced rank;
- Wald inference is available but labelled as low reliability;
- Satterthwaite and Kenward-Roger routes are refused because the
  variance-parameter derivatives are not defined at the boundary;
- parametric bootstrap testing is available for the fixed effect;
- bootstrap confidence intervals are available, but they use a full-model
  estimator-distribution target, not a fixed-effect null target;
- cluster bootstrap currently reports estimator distributions, not fixed-effect
  p-values;
- boundary LRTs are available for variance-component tests, not fixed-effect
  tests.

The main design principle is that `mixeff` does not make the analyst guess.
Every route either returns a labelled result or refuses with a stable reason.

## References

- Crainiceanu, C. M., and Ruppert, D. (2004). Likelihood ratio tests
  in linear mixed models with one variance component. *Journal of the
  Royal Statistical Society: Series B*, 66(1), 165–185.
  <https://doi.org/10.1111/j.1467-9868.2004.00438.x>
- Halekoh, U., and Højsgaard, S. (2014). A Kenward-Roger approximation
  and parametric bootstrap methods for tests in linear mixed models —
  the `R` package `pbkrtest`. *Journal of Statistical Software*,
  59(9), 1–32. <https://doi.org/10.18637/jss.v059.i09>
- Kenward, M. G., and Roger, J. H. (1997). Small sample inference for
  fixed effects from restricted maximum likelihood. *Biometrics*,
  53(3), 983–997. <https://doi.org/10.2307/2533558>
- Self, S. G., and Liang, K.-Y. (1987). Asymptotic properties of
  maximum likelihood estimators and likelihood ratio tests under
  nonstandard conditions. *Journal of the American Statistical
  Association*, 82(398), 605–610.
  <https://doi.org/10.1080/01621459.1987.10478472>
- Stram, D. O., and Lee, J. W. (1994). Variance components testing in
  the longitudinal mixed effects model. *Biometrics*, 50(4),
  1171–1177. <https://doi.org/10.2307/2533455>
