inst/doc/vetresearchlmm-examples.R

In VetResearchLMM: Linear Mixed Models - An Introduction with Applications in Veterinary Research

## -----------------------------------------------------------------------------
#| 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"))


## -----------------------------------------------------------------------------
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)
}


## -----------------------------------------------------------------------------
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)
}


## -----------------------------------------------------------------------------
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"
  )
}


## -----------------------------------------------------------------------------
if (exists("example_pairs")) {
  as.data.frame(example_pairs)
}


## -----------------------------------------------------------------------------
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()
}


## -----------------------------------------------------------------------------
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)
}