inst/doc/vetresearchlmm-methodology.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"))


## -----------------------------------------------------------------------------
if (requireNamespace("lme4", quietly = TRUE)) {
  fit_ml <- lme4::lmer(
    Pcv ~ dose * Drug + (1 | Region / Drug),
    data = ex125,
    REML = FALSE
  )
  if (requireNamespace("report", quietly = TRUE)) {
    report::report(fit_ml)
  }

  fit_reml <- lme4::lmer(
    Pcv ~ dose * Drug + (1 | Region / Drug),
    data = ex125,
    REML = TRUE
  )
  if (requireNamespace("report", quietly = TRUE)) {
    report::report(fit_reml)
  }

  data.frame(
    component = as.data.frame(lme4::VarCorr(fit_ml))$grp,
    ml = round(as.data.frame(lme4::VarCorr(fit_ml))$vcov, 3),
    reml = round(as.data.frame(lme4::VarCorr(fit_reml))$vcov, 3)
  )
}


## -----------------------------------------------------------------------------
if (requireNamespace("report", quietly = TRUE) && exists("fit_reml")) {
  report::report(fit_reml)
} else {
  data.frame(
    workflow = "Optional model report",
    requirement = "Install report to generate easystats-style summaries"
  )
}


## -----------------------------------------------------------------------------
if (requireNamespace("emmeans", quietly = TRUE) && exists("fit_reml")) {
  dose_emm <- emmeans::emmeans(
    fit_reml,
    ~ dose | Drug,
    lmer.df = "asymptotic"
  )
  dose_pairs <- emmeans::contrast(dose_emm, method = "pairwise")

  as.data.frame(dose_emm)
} else {
  data.frame(
    workflow = "Optional marginal means",
    requirement = "Install emmeans to compute post hoc model summaries"
  )
}


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