---
title: "Using lmmix"
output: rmarkdown::html_vignette
bibliography: references.bib
link-citations: true
vignette: >
  %\VignetteIndexEntry{Using lmmix}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Model scope

`lmmix` fits Gaussian linear mixed models with random effects, structured
residual covariance, or both components in the same model. The package
directly optimizes a profiled ML or REML objective and does not delegate
estimation to another mixed-model package.

This vignette focuses on practical use. The mathematical formulation,
covariance parameterizations, likelihood criteria, and inference derivations
are presented separately in the vignette *Estimation and inference in lmmix*.

```{r}
library(lmmix)
```

## Prepare the data

The included `multicentre` data have three treatments, three measurement
occasions, three centers and subjects nested within center and treatment.

```{r data-summary}
str(multicentre)
colSums(is.na(multicentre))
```

The fixed, random and repeated formulas jointly determine the variables used
by the analysis. Missing rows are handled with `na.omit`, `na.exclude`, or
`na.fail`. The default is `na.omit`; `na.exclude` restores omitted positions in
fitted values, residuals, and augmented data.

## Fit a combined covariance model

This model contains a random center intercept and an AR(1) residual covariance
within subject.

```{r combined-fit}
fit <- lmm(
  data = multicentre,
  formula = Y ~ Drug * Time,
  random = ~ 1 | Center,
  repeated = ~ Time | Center:Drug:Subject,
  structure = "ar1",
  method = "REML",
  ddf = "satterthwaite"
)

fit
```

The three formulas have distinct roles:

| Argument | Role in the fitted model |
|:--|:--|
| `formula` | Defines the response and fixed-effects matrix `X` |
| `random` | Defines the random-effects matrix `Z` and grouping factor for `G` |
| `repeated` | Defines the ordering variable and independent blocks for `R` |

For the residual structures other than `id`, each repeated-measures group must
have at most one observation per ordering value.

## Inspect convergence before inference

The optimizer code and Hessian diagnostic are part of the fitted object.

```{r convergence}
fit$convergence[c(
  "code",
  "message",
  "optimizer",
  "iterations",
  "hessian_positive_definite"
)]
fit$convergence$attempts
```

A nonzero convergence code or a non-positive-definite likelihood Hessian
requires investigation before interpreting standard errors or denominator
degrees of freedom. The fitting function emits a warning in either case.

## Fixed effects and type III tests

`summary()` combines coefficient-level inference, type III tests, covariance
parameters and information criteria.

```{r summary}
summary(fit)
```

Individual components are also available through standard methods.

```{r inference}
fixef(fit)
vcov(fit)
anova(fit)
confint(fit)
confint(fit, parm = "theta_")
```

`anova()` performs type III fixed-effect tests for one model. With multiple
models it performs likelihood-ratio comparisons, automatically refitting REML
models with ML when their fixed effects differ. Satterthwaite, Kenward-Roger,
and residual denominator degrees of freedom are available. Kenward-Roger is
restricted to REML fits.

```{r kenward-roger-and-comparison, eval=FALSE}
fit_kr <- lmm(
  data = multicentre,
  formula = Y ~ Drug * Time,
  random = ~ 1 | Center,
  repeated = ~ Time | Center:Drug:Subject,
  structure = "ar1",
  ddf = "kenward-roger"
)

fit_reduced <- lmm(
  data = multicentre,
  formula = Y ~ Drug + Time,
  random = ~ 1 | Center,
  repeated = ~ Time | Center:Drug:Subject,
  structure = "ar1"
)
anova(fit_reduced, fit)

# Use a simulated null reference for a covariance comparison.
anova(
  smaller_ml_fit,
  larger_ml_fit,
  test = "parametric.bootstrap",
  nsim = 999,
  seed = 2026
)
```

## Covariance parameters and random effects

`VarCorr()` returns the estimated covariance components on their natural
scale. `ranef()` always returns a named list with one tibble per random term.
This stable contract also applies when the model has only one random term; a
marginal model returns an empty list.

```{r covariance}
VarCorr(fit)
ranef(fit)
ranef(fit)[[1L]]
```

For a random-slope model, the random-effect covariance is unstructured. Its
variances and correlations appear as separate rows in `VarCorr()`.

## Estimated marginal means

Marginal means use equal weights over nuisance-factor levels. Numeric
covariates not requested in `specs` are held at their observed means. This
construction targets population marginal means in the sense of @searle1980.

```{r lsmeans}
lsmeans(fit, ~Drug)
lsmeans(fit, ~Time)
```

Pairwise contrast matrices are generated automatically. The `adjust` argument
is passed to `stats::p.adjust()`. When p-values are adjusted,
`conf_adjust = "auto"` reports simultaneous Bonferroni confidence intervals.
Use `conf_adjust = "none"` for pointwise intervals.

```{r pairwise}
lsmeans(fit, pairwise ~ Drug, adjust = "holm")
```

Reference values for numeric covariates can be overridden with `at`:

```{r at-example, eval=FALSE}
lsmeans(fit, ~Drug, at = list(Time = "2"))
```

When `emmeans` is installed, `lmmix` supplies the required fixed-effects basis
and model-specific denominator degrees of freedom.

```{r emmeans, eval=requireNamespace("emmeans", quietly = TRUE)}
emmeans::emmeans(fit, ~Drug)
```

## Predictions, residuals and tidy output

Conditional fitted values include empirical random effects. Marginal fitted
values use fixed effects only.

```{r predictions}
head(fitted(fit, type = "conditional"))
head(fitted(fit, type = "marginal"))
head(residuals(fit, type = "pearson"))
```

For `predict(fit, newdata)`, known grouping levels receive their fitted random
effect. New groups cause an error unless `allow.new.levels = TRUE` is supplied;
when allowed, their random contribution is zero. Use `re.form = NA` or
`re.form = ~ 0` for fixed-effects-only predictions. A fitted random-effects
formula, or a list of fitted formulas, selects a subset of random terms.

Standard errors and intervals are available directly. Confidence intervals
describe the expected response. Prediction intervals additionally include the
fitted residual variance. Both condition on empirical random effects selected
through `re.form`; BLUP estimation uncertainty is not added.

```{r prediction-intervals, eval=FALSE}
predict(fit, se.fit = TRUE, re.form = NA)
predict(fit, interval = "confidence", re.form = NA)
predict(fit, interval = "prediction", level = 0.90, re.form = NA)
```

## Simulation, updating, weights, and offsets

`simulate()` draws marginal Gaussian response vectors from the fitted mean and
covariance. Supplying a seed makes the result reproducible without changing the
caller's random-number state. `update()` refits a modified formula or model
argument.

```{r simulation-update, eval=FALSE}
simulate(fit, nsim = 50, seed = 2026)
update(fit, . ~ . - Drug:Time)
update(fit, structure = "cs")
```

Known relative precision weights rescale the residual covariance on both sides
by `diag(1 / sqrt(weights))`. Formula offsets, explicit offsets, and custom
fixed-effect contrasts are also supported.

```{r weighted-offset-model, eval=FALSE}
weighted_fit <- lmm(
  data = analysis_data,
  formula = response ~ treatment + offset(baseline),
  random = ~ 1 | subject,
  weights = precision,
  contrasts = list(treatment = "contr.sum")
)
```

## Diagnostic plots

The S3 method `plot.lmm()` returns a `ggplot2` object. Three diagnostics are
available through `which`: standardized residuals against fitted values, a
normal Q-Q plot, and observed against fitted values.

```{r diagnostic-plots, fig.show='hold', out.width='48%', fig.alt='Residual and normal Q-Q diagnostic plots'}
plot(fit, which = "residuals")
plot(fit, which = "qq")
```

Because the returned object is a `ggplot`, layers and themes can be added in
the usual way.

```{r customize-plot, fig.alt='Observed values plotted against fitted values'}
plot(fit, which = "fitted") +
  ggplot2::theme_minimal()
```

For correlated residual structures, the plotted Pearson residuals are scaled
by their marginal residual standard deviations. They are not decorrelated
innovations, so residual dependence should also be assessed within the
repeated-measures groups.

The `generics` methods return tibbles.

Printed tables use readable headings without dots, such as `Std Error`,
`p value`, `Conf Low`, and `Conf High`. The underlying tibbles retain the
standard `broom` column names required for programmatic use. Input column names
are preserved by `augment()` and `ranef()`.

```{r broom}
generics::tidy(fit)
generics::glance(fit)
head(generics::augment(fit))
```

## Alternative model specifications

The same interface covers marginal, random-effects and combined models.

```{r alternatives, eval=FALSE}
# Marginal model with unstructured repeated-measures covariance
fit_marginal <- lmm(
  data,
  response ~ treatment * time,
  repeated = ~ time | subject,
  structure = "un"
)

# Random intercept and slope with independent residuals
fit_random <- lmm(
  data,
  response ~ treatment + time,
  random = ~ 1 + time | subject,
  structure = "id"
)

# Independent crossed random intercepts
fit_crossed <- lmm(
  data,
  response ~ treatment + time,
  random = list(
    center = ~ 1 | center,
    subject = ~ 1 | subject
  )
)

# Combined random intercept and Toeplitz residual covariance
fit_combined <- lmm(
  data,
  response ~ treatment * time,
  random = ~ 1 | center,
  repeated = ~ time | center:subject,
  structure = "toep"
)

# Four-band Toeplitz covariance, including the main diagonal
fit_banded <- lmm(
  data,
  response ~ treatment * time,
  random = ~ 1 | center,
  repeated = ~ time | center:subject,
  structure = "toep(4)"
)
```

## Practical limits

Version `0.4.0` identifies independent connected components induced by the
random and repeated grouping structures. The likelihood is evaluated by
factoring one covariance matrix per component and aggregating their GLS
contributions. This substantially reduces work for many subject-level or
nested designs.

The fitted object still stores the complete dense marginal covariance, and
some inference calculations remain dense. Crossed effects can connect all
observations into one component. Consequently, `lmmix` is not yet a
large-scale sparse mixed-model engine and does not promise a fixed observation
limit. The practical cost depends on the size of the largest connected
component, the number of covariance parameters, and the requested inference.
Generalized responses remain outside the Gaussian model scope.