---
title: "Saving and Reloading Fits"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteIndexEntry{Saving and Reloading Fits}
  %\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)
```

A mixed-model fit often outlives the R session that produced it.
You fit a model today, save it with the analysis, hand the script to
a collaborator, and reopen the fit six months later for a contrast,
a revision, or a referee response.

`mixeff` stores the fitted values, the random-effects design, the
convergence record, and the inference labels inside the R object —
so each of those tasks works after `readRDS()` without recomputing
the fit, and without depending on the original Rust handle.

## Fit a model

```{r clinic-data, include = FALSE}
set.seed(42)
n_clinic <- 12L
n_per_clinic <- 6L
clinic <- factor(rep(seq_len(n_clinic), each = n_per_clinic))
week <- rep(0:(n_per_clinic - 1L), n_clinic)
treatment <- factor(
  rep(rep(c("usual", "coached"), each = n_per_clinic), length.out = length(week)),
  levels = c("usual", "coached")
)
clinic_shift <- rnorm(n_clinic, sd = 0.45)
score <- 7.2 - 0.28 * week - 0.55 * (treatment == "coached") +
  clinic_shift[as.integer(clinic)] + rnorm(length(week), sd = 0.35)
clinic_visits <- data.frame(score, week, treatment, clinic)
```

```{r fit-model}
fit <- lmm(
  score ~ week + treatment + (1 | clinic),
  clinic_visits,
  control = mm_control(verbose = -1)
)
```

Before saving, the ordinary extractors work as expected.

```{r before-save}
fixef(fit)
reporting_table(fit, "fixed_effects")
```

## Round trip through RDS

```{r rds-roundtrip}
path <- tempfile(fileext = ".rds")
saveRDS(fit, path)

restored <- readRDS(path)
restored <- revive(restored)
```

The restored object still answers the same fitted-model questions.

```{r after-restore}
fixef(restored)
head(predict(restored))
reporting_table(restored, "fixed_effects")
```

```{r restore-checks, include = FALSE}
stopifnot(identical(fixef(restored), fixef(fit)))
stopifnot(length(predict(restored)) == nrow(clinic_visits))
stopifnot(identical(reporting_table(restored, "fixed_effects"),
                    reporting_table(fit, "fixed_effects")))
```

## Rebuild design matrices when needed

Design extractors can be rebuilt from the stored formula and model frame.

```{r design-extractors}
X <- model.matrix(restored, type = "fixed")
Z <- model.matrix(restored, type = "random")

dim(X)
dim(Z)
class(Z)
```

`getME()` provides a small familiar subset for code that expects lme4-style
names.

```{r getme}
getME(restored, c("theta", "beta", "cnms"))
```

```{r design-checks, include = FALSE}
stopifnot(nrow(X) == nrow(clinic_visits))
stopifnot(nrow(Z) == nrow(clinic_visits))
stopifnot(inherits(Z, "dgCMatrix"))
```

## What stays explicit?

Quantities that the Rust inference contract cannot certify are marked rather
than fabricated. For full-rank fits, `vcov()` returns the model-based
fixed-effect covariance from the upstream `fixed_effect_covariance_matrix`
payload. For rank-deficient or otherwise uncertified fits, the matrix carries
an `mm_unavailable_reason` attribute and the values are `NA`.

```{r covariance-status}
V <- vcov(restored)
attr(V, "mm_unavailable_reason")
V
```

Conditional variances for random effects also survive the round trip. With
`condVar = TRUE`, each grouping table carries a finite `postVar` array.

```{r ranef-condvar}
re <- ranef(restored, condVar = TRUE)
attr(re, "mm_unavailable_reason")
dim(attr(re$clinic, "postVar"))
```

```{r unavailable-checks, include = FALSE}
stopifnot(is.null(attr(V, "mm_unavailable_reason")) ||
          is.character(attr(V, "mm_unavailable_reason")))
stopifnot(is.null(attr(re, "mm_unavailable_reason")))
stopifnot(is.array(attr(re$clinic, "postVar")))
stopifnot(all(is.finite(attr(re$clinic, "postVar"))))
```

## What should you report?

Use `reporting_table()` for durable tables and `summary()` for console output.
Both continue to work after a save/load cycle.

```{r final-reporting}
coef_table <- summary(restored, method = "auto")$coefficients
knitr::kable(coef_table, digits = 4)
reporting_table(restored, "overview")
```