---
title: "Book Examples in R"
format:
  html:
    toc: true
    number-sections: true
    css: styles.css
execute:
  warning: false
  message: false
---

<!--
%\VignetteIndexEntry{Book Examples in R}
%\VignetteEngine{quarto::html}
%\VignetteEncoding{UTF-8}
-->

```{r}
#| include: false
library(methods)
if (requireNamespace("lme4", quietly = TRUE)) library(lme4)
if (requireNamespace("lmerTest", quietly = TRUE)) library(lmerTest)
load(file.path("..", "data", "ex125.RData"))
load(file.path("..", "data", "ex127.RData"))
load(file.path("..", "data", "ex31.RData"))
```

## Chapter and Example Map

The package implements the following book examples as help pages:

| Package page | Book page | Dataset | Main analysis |
|---|---:|---|---|
| `Examp1.3.2` | 16 | `ex124` | fixed effect ANOVA with split-plot structure |
| `Examp2.2.1.7` | 42 | `ex121` | dose hypothesis test |
| `Examp2.4.2.2` | 64 | `ex125` | ML and REML variance components |
| `Examp2.4.3.1` | 66 | `ex127` | BLUPs for sire effects |
| `Examp2.5.1.1` | 67 | `ex125` | fixed effect estimates |
| `Examp2.5.2.1` | 68 | `ex125` | linear combinations of fixed effects |
| `Examp2.5.3.1` | 70 | `ex125` | confidence intervals |
| `Examp2.5.4.1` | 74 | `ex125` | fixed versus mixed model intervals |
| `Examp2.6.1` | 76 | `ex125` | fixed effect hypothesis tests |
| `Examp3.1` | 80, 84, 86 | `ex31` | three PCV response models |
| `Examp3.2` | 88 | `ex32` | weaning weight mixed model |
| `Examp3.3` | 88 onward | `ex33` | longitudinal PCV mixed models |

## Example 2.4.3.1

The sire random-effect example estimates a sire variance component and predicts
sire effects.

```{r}
if (requireNamespace("lme4", quietly = TRUE)) {
  sire_fit <- lme4::lmer(Ww ~ 1 + (1 | sire), data = ex127, REML = TRUE)
  if (requireNamespace("report", quietly = TRUE)) {
    report::report(sire_fit)
  }

  lme4::fixef(sire_fit)
  as.data.frame(lme4::VarCorr(sire_fit))
  head(lme4::ranef(sire_fit)$sire)
}
```

The readable book targets are reproduced to rounding: mean near `13.97`, sire
variance near `3.685`, residual variance near `3.542`, and sire 4 BLUP near
`3.16`.

## Example 2.6.1

The contrast for the Samorin versus Berenil comparison can be expressed against
the fitted fixed-effect vector.

```{r}
if (requireNamespace("lmerTest", quietly = TRUE)) {
  contrast_fit <- lmerTest::lmer(
    Pcv ~ dose * Drug + (1 | Region / Drug),
    data = ex125,
    REML = TRUE,
    contrasts = list(dose = "contr.SAS", Drug = "contr.SAS")
  )
  if (requireNamespace("report", quietly = TRUE)) {
    report::report(contrast_fit)
  }

  contrast <- matrix(
    c(0, 0, -1, -0.5),
    ncol = 4,
    dimnames = list("drug_difference", rownames(summary(contrast_fit)$coef))
  )

  lmerTest::contest(contrast_fit, contrast, joint = FALSE)
}
```

The estimate and variance match the readable book values (`-5.55` and about
`0.478`). The p-value printed by current R is based on the same F statistic but
can differ from the book because of software and reporting differences.

## Post Hoc Inference with emmeans

After fitting a mixed model, `report` and `emmeans` provide complementary
interpretations. `report` gives a narrative summary of the model object, while
`emmeans` gives estimated marginal means, contrasts, and post hoc comparisons.
For the split-plot example, the marginal means can be computed for dose within
drug. The package helper `emmeans_mixed_model()` returns these same marginal
means or pairwise contrasts through a guarded package-level API.

```{r}
if (requireNamespace("emmeans", quietly = TRUE) && exists("contrast_fit")) {
  example_emm <- emmeans::emmeans(
    contrast_fit,
    ~ dose | Drug,
    lmer.df = "asymptotic"
  )
  example_pairs <- emmeans::contrast(example_emm, method = "pairwise")

  as.data.frame(example_emm)
} else {
  data.frame(
    workflow = "Optional marginal means",
    requirement = "Install emmeans to compute estimated marginal means"
  )
}
```

```{r}
if (exists("example_pairs")) {
  as.data.frame(example_pairs)
}
```

```{r}
if (exists("example_emm") && requireNamespace("ggplot2", quietly = TRUE)) {
  example_emm_df <- as.data.frame(example_emm)
  lower_ci <- intersect(c("lower.CL", "asymp.LCL"), names(example_emm_df))[1L]
  upper_ci <- intersect(c("upper.CL", "asymp.UCL"), names(example_emm_df))[1L]
  example_emm_df$lower_ci <- example_emm_df[[lower_ci]]
  example_emm_df$upper_ci <- example_emm_df[[upper_ci]]

  ggplot2::ggplot(
    example_emm_df,
    ggplot2::aes(x = dose, y = emmean, color = Drug, group = Drug)
  ) +
    ggplot2::geom_line(linewidth = 0.6) +
    ggplot2::geom_point(size = 2.4) +
    ggplot2::geom_errorbar(
      ggplot2::aes(ymin = lower_ci, ymax = upper_ci),
      width = 0.06,
      linewidth = 0.5
    ) +
    ggplot2::labs(
      x = "Dose",
      y = "Estimated marginal mean PCV",
      color = "Drug",
      title = "Model-based marginal means by drug and dose"
    ) +
    ggplot2::theme_minimal()
}
```

## Modern Preprocessing Pattern

Several book examples require factor versions of design columns or indicator
variables for contrasts. The package examples use a single `collapse::fmutate()`
step so that these derived variables are readable and reproducible.

```{r}
if (requireNamespace("collapse", quietly = TRUE)) {
  ex31_prepared <-
    ex31 |>
    collapse::fmutate(
      herd1 = factor(herd),
      drug1 = factor(drug),
      dose1 = factor(dose),
      ber = as.integer(drug == "BERENIL"),
      ber1 = as.integer(drug == "BERENIL" & dose == 1L),
      pcv_ber1 = PCV1 * as.integer(drug == "BERENIL" & dose == 1L),
      pcv_ber2 = PCV1 * as.integer(drug == "BERENIL" & dose == 2L)
    )
  head(ex31_prepared)
}
```

## Example 3.1

Example 3.1 fits increasingly rich models for packed cell volume one month
after treatment. The current package data contain 38 observations, while the
book page image for Model 3 shows denominator degrees of freedom greater than
300. That indicates the printed SAS output and the packaged data are not at the
same observational granularity for that table. The package therefore verifies
the R code logically and statistically but does not claim exact numerical
agreement for the Model 3 printed table.