inst/doc/introduction_comparisons_2.R

In modelbased: Estimation of Model-Based Predictions, Contrasts and Means

## ----set-options, echo = FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "png",
  fig.width = 7,
  fig.height = 3.5,
  message = FALSE,
  warning = FALSE,
  package.startup.message = FALSE
)
options(width = 800)
options(modelbased_join_dots = FALSE)
ht6 <- m <- NULL

arrow_color <- "#FF00cc"

pkgs <- c("ggplot2", "see", "marginaleffects", "parameters")

if (!all(insight::check_if_installed(pkgs, quietly = TRUE))) {
  knitr::opts_chunk$set(eval = FALSE)
}

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(modelbased)
library(parameters)
data(iris)
m <- lm(Sepal.Width ~ Sepal.Length + Species, data = iris)
model_parameters(m)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_means(m, "Sepal.Length=c(4,5,6,7)")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_slopes(m, "Sepal.Length")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m <- lm(Sepal.Width ~ Sepal.Length * Species, data = iris)
pred <- estimate_means(m, c("Sepal.Length", "Species"))
plot(pred)

## ----echo=FALSE---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
library(ggplot2)
p <- plot(pred)
dat <- as.data.frame(pred)
dat1 <- data.frame(
  x = dat$Sepal.Length[c(7, 19)],
  y = dat$Mean[c(7, 19)],
  group_col = "setosa",
  stringsAsFactors = FALSE
)
dat2 <- data.frame(
  x = dat$Sepal.Length[c(8, 20)],
  y = dat$Mean[c(8, 23)],
  group_col = "versicolor",
  stringsAsFactors = FALSE
)
dat3 <- data.frame(
  x = dat$Sepal.Length[c(15, 27)],
  y = dat$Mean[c(15, 27)],
  group_col = "virginica",
  stringsAsFactors = FALSE
)
p + annotate(
  "segment",
  x = dat1$x[1], xend = dat1$x[2], y = dat1$y[1], yend = dat1$y[1],
  color = "orange", linewidth = 1, linetype = "dashed"
) + annotate(
  "segment",
  x = dat1$x[2], xend = dat1$x[2], y = dat1$y[1], yend = dat1$y[2],
  color = "orange", linewidth = 1, linetype = "dashed"
) + annotate(
  "segment",
  x = dat1$x[1], xend = dat1$x[2], y = dat1$y[1], yend = dat1$y[2],
  color = "orange", linewidth = 1
) + annotate(
  "segment",
  x = dat2$x[1], xend = dat2$x[2], y = dat2$y[1], yend = dat2$y[1],
  color = "orange", linewidth = 1, linetype = "dashed"
) + annotate(
  "segment",
  x = dat2$x[2], xend = dat2$x[2], y = dat2$y[1], yend = dat2$y[2],
  color = "orange", linewidth = 1, linetype = "dashed"
) + annotate(
  "segment",
  x = dat2$x[1], xend = dat2$x[2], y = dat2$y[1], yend = dat2$y[2],
  color = "orange", linewidth = 1
) + annotate(
  "segment",
  x = dat3$x[1], xend = dat3$x[2], y = dat3$y[1], yend = dat3$y[1],
  color = "orange", linewidth = 1, linetype = "dashed"
) + annotate(
  "segment",
  x = dat3$x[2], xend = dat3$x[2], y = dat3$y[1], yend = dat3$y[2],
  color = "orange", linewidth = 1, linetype = "dashed"
) + annotate(
  "segment",
  x = dat3$x[1], xend = dat3$x[2], y = dat3$y[1], yend = dat3$y[2],
  color = "orange", linewidth = 1
) +
  ggtitle("Linear trend of `Sepal.Length` by `Species`")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_slopes(m, "Sepal.Length", by = "Species")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(m, "Sepal.Length", by = "Species")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_slopes(m, "Sepal.Length", by = "Species")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(
  m,
  "Sepal.Length",
  by = "Species",
  comparison = "(b1 - b2) = (b1 - b3)"
)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
m <- lm(Sepal.Width ~ Petal.Length * Petal.Width, data = iris)
pred <- estimate_means(m, c("Petal.Length", "Petal.Width=[sd]"))
plot(pred)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_slopes(m, "Petal.Length", by = "Petal.Width=[sd]")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_contrasts(m, "Petal.Length", by = "Petal.Width=[sd]", digits = 1)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
states <- as.data.frame(state.x77)
states$HSGrad <- states$`HS Grad`
m_mod <- lm(Income ~ HSGrad + Murder * Illiteracy, data = states)

pr <- estimate_means(m_mod, c("Murder", "Illiteracy"))
plot(pr)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
estimate_slopes(m_mod, "Murder", by = "Illiteracy")

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# we will force to calculate slopes at 200 values for "Illiteracy" using `length`
slopes <- estimate_slopes(m_mod, "Murder", by = "Illiteracy", length = 200)
summary(slopes)

## ----message=FALSE------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
plot(slopes)

## -----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
pr <- estimate_means(m_mod, c("Murder", "Illiteracy=c(0.7,1.5,2.8)"))
plot(pr) + ggplot2::facet_wrap(~Illiteracy)