Automated plotting for 'modelbased' objects — visualisation_recipe.estimate

Most 'modelbased' objects can be visualized using the plot() function, which internally calls the visualisation_recipe() function. See the examples below for more information and examples on how to create and customize plots.

Usage

# S3 method for class 'estimate_predicted'
visualisation_recipe(
  x,
  show_data = FALSE,
  show_ci = TRUE,
  point = NULL,
  line = NULL,
  pointrange = NULL,
  ribbon = NULL,
  facet = NULL,
  grid = NULL,
  join_dots = TRUE,
  ...
)

# S3 method for class 'estimate_slopes'
visualisation_recipe(
  x,
  show_ci = TRUE,
  line = NULL,
  pointrange = NULL,
  ribbon = NULL,
  facet = NULL,
  grid = NULL,
  ...
)

# S3 method for class 'estimate_grouplevel'
visualisation_recipe(
  x,
  show_ci = TRUE,
  line = NULL,
  pointrange = NULL,
  ribbon = NULL,
  facet = NULL,
  grid = NULL,
  ...
)

Arguments

x: A modelbased object.
show_data: Logical, if TRUE, display the "raw" data as a background to the model-based estimation.
show_ci: Logical, if TRUE, show error bars or uncertainty bands.
point, line, pointrange, ribbon, facet, grid: Additional aesthetics and parameters for the geoms (see customization example).
join_dots: Logical, if TRUE and for categorical focal terms in by, dots (estimates) are connected by lines, i.e. plots will be a combination of dots with error bars and connecting lines. If FALSE, only dots and error bars are shown.
...: Not used.

Details

The plotting works by mapping any predictors from the by argument to the x-axis, colors, alpha (transparency) and facets. Thus, the appearance of the plot depends on the order of the variables that you specify in the by argument. For instance, the plots corresponding to estimate_relation(model, by=c("Species", "Sepal.Length")) and estimate_relation(model, by=c("Sepal.Length", "Species")) will look different.

The automated plotting is primarily meant for convenient visual checks, but for publication-ready figures, we recommend re-creating the figures using the ggplot2 package directly.

Examples

# ==============================================
# estimate_relation, estimate_expectation, ...
# ==============================================
# Simple Model ---------------
x <- estimate_relation(lm(mpg ~ wt, data = mtcars))
layers <- visualisation_recipe(x)
layers
#> Layer 1
#> --------
#> Geom type: ribbon
#> data = [10 x 6]
#> aes_string(
#>   y = 'Predicted'
#>   x = 'wt'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#>   group = '.group'
#> )
#> alpha = 0.3333333
#> 
#> Layer 2
#> --------
#> Geom type: line
#> data = [10 x 6]
#> aes_string(
#>   y = 'Predicted'
#>   x = 'wt'
#>   group = '.group'
#> )
#> 
#> Layer 3
#> --------
#> Geom type: labs
#> y = 'Predicted value of mpg'
#> 
plot(layers)


# visualization_recipe() is called implicitly when you call plot()
plot(estimate_relation(lm(mpg ~ qsec, data = mtcars)))


# And can be used in a pipe workflow
lm(mpg ~ qsec, data = mtcars) |>
  estimate_relation(ci = c(0.5, 0.8, 0.9)) |>
  plot()


# Customize aesthetics ----------

plot(x,
  point = list(color = "red", alpha = 0.6, size = 3),
  line = list(color = "blue", size = 3),
  ribbon = list(fill = "green", alpha = 0.7)
) +
  theme_minimal() +
  labs(title = "Relationship between MPG and WT")



# Customize raw data -------------

plot(x, point = list(geom = "density_2d_filled"), line = list(color = "white")) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(legend.position = "none")


# Single predictors examples -----------

plot(estimate_relation(lm(Sepal.Length ~ Species, data = iris)))


# 2-ways interaction ------------

# Numeric * numeric
x <- estimate_relation(lm(mpg ~ wt * qsec, data = mtcars))
plot(x)


# Numeric * factor
x <- estimate_relation(lm(Sepal.Width ~ Sepal.Length * Species, data = iris))
plot(x)


# ==============================================
# estimate_means
# ==============================================
# Simple Model ---------------
x <- estimate_means(lm(Sepal.Width ~ Species, data = iris), by = "Species")
layers <- visualisation_recipe(x)
layers
#> Layer 1
#> --------
#> Geom type: line
#> data = [3 x 8]
#> aes_string(
#>   y = 'Mean'
#>   x = 'Species'
#>   group = '.group'
#> )
#> 
#> Layer 2
#> --------
#> Geom type: pointrange
#> data = [3 x 8]
#> aes_string(
#>   y = 'Mean'
#>   x = 'Species'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#>   group = '.group'
#> )
#> 
#> Layer 3
#> --------
#> Geom type: labs
#> y = 'Mean of Sepal.Width'
#> 
plot(layers)


# Customize aesthetics
layers <- visualisation_recipe(x,
  point = list(width = 0.03, color = "red"),
  pointrange = list(size = 2, linewidth = 2),
  line = list(linetype = "dashed", color = "blue")
)
plot(layers)


# Two levels ---------------
data <- mtcars
data$cyl <- as.factor(data$cyl)

model <- lm(mpg ~ cyl * wt, data = data)

x <- estimate_means(model, by = c("cyl", "wt"))
plot(x)



# GLMs ---------------------
data <- data.frame(vs = mtcars$vs, cyl = as.factor(mtcars$cyl))
x <- estimate_means(glm(vs ~ cyl, data = data, family = "binomial"), by = c("cyl"))
plot(x)

# ==============================================
# estimate_slopes
# ==============================================
model <- lm(Sepal.Width ~ Species * Petal.Length, data = iris)
x <- estimate_slopes(model, trend = "Petal.Length", by = "Species")

layers <- visualisation_recipe(x)
layers
#> Layer 1
#> --------
#> Geom type: line
#> data = [3 x 8]
#> aes_string(
#>   y = 'Slope'
#>   x = 'Species'
#>   group = '.group'
#> )
#> 
#> Layer 2
#> --------
#> Geom type: pointrange
#> data = [3 x 8]
#> aes_string(
#>   y = 'Slope'
#>   x = 'Species'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#>   group = '.group'
#> )
#> 
#> Layer 3
#> --------
#> Geom type: labs
#> y = 'Slope of Sepal.Width'
#> 
plot(layers)


# Customize aesthetics and add horizontal line and theme
layers <- visualisation_recipe(x, pointrange = list(size = 2, linewidth = 2))
plot(layers) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(y = "Effect of Petal.Length", title = "Marginal Effects")


model <- lm(Petal.Length ~ poly(Sepal.Width, 4), data = iris)
x <- estimate_slopes(model, trend = "Sepal.Width", by = "Sepal.Width", length = 20)
plot(visualisation_recipe(x))


model <- lm(Petal.Length ~ Species * poly(Sepal.Width, 3), data = iris)
x <- estimate_slopes(model, trend = "Sepal.Width", by = c("Sepal.Width", "Species"))
plot(visualisation_recipe(x))

# ==============================================
# estimate_grouplevel
# ==============================================
data <- lme4::sleepstudy
data <- rbind(data, data)
data$Newfactor <- rep(c("A", "B", "C", "D"))

# 1 random intercept
model <- lme4::lmer(Reaction ~ Days + (1 | Subject), data = data)
x <- estimate_grouplevel(model)
layers <- visualisation_recipe(x)
layers
#> Layer 1
#> --------
#> Geom type: pointrange
#> data = [18 x 9]
#> aes_string(
#>   y = 'Coefficient'
#>   x = 'Level'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#>   group = '.group'
#> )
#> 
#> Layer 2
#> --------
#> Geom type: coord_flip
#> 
plot(layers)


# 2 random intercepts
model <- lme4::lmer(Reaction ~ Days + (1 | Subject) + (1 | Newfactor), data = data)
x <- estimate_grouplevel(model)
plot(x) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  theme_minimal()

# Note: we need to use hline instead of vline because the axes is flipped

model <- lme4::lmer(Reaction ~ Days + (1 + Days | Subject) + (1 | Newfactor), data = data)
x <- estimate_grouplevel(model)
plot(x)