Visualisation Recipe for 'modelbased' Objects

Visualisation Recipe for 'modelbased' Objects

Usage

# S3 method for class 'estimate_grouplevel'
visualisation_recipe(
  x,
  hline = NULL,
  pointrange = NULL,
  facet_wrap = NULL,
  labs = NULL,
  ...
)

# S3 method for class 'estimate_means'
visualisation_recipe(
  x,
  show_data = "jitter",
  point = NULL,
  jitter = point,
  boxplot = NULL,
  violin = NULL,
  line = NULL,
  pointrange = NULL,
  labs = NULL,
  ...
)

# S3 method for class 'estimate_predicted'
visualisation_recipe(
  x,
  show_data = "points",
  point = NULL,
  density_2d = NULL,
  line = NULL,
  ribbon = NULL,
  labs = NULL,
  ...
)

# S3 method for class 'estimate_slopes'
visualisation_recipe(
  x,
  hline = NULL,
  line = NULL,
  pointrange = NULL,
  ribbon = NULL,
  labs = NULL,
  facet_wrap = NULL,
  ...
)

Arguments

x

A modelbased object.

...

Other arguments passed to other functions.

show_data

Display the "raw" data as a background to the model-based estimation. Can be set to "none" to remove it. When input is the result of estimate_means, show_data can be "points" (the jittered observation points), "boxplot", "violin" a combination of them (see examples). When input is the result of estimate_expectation or estimate_relation, show_data can be "points" (the points of the original data corresponding to the x and y axes), "density_2d", "density_2d_filled", "density_2d_polygon" or "density_2d_raster".

point, jitter, boxplot, violin, pointrange, density_2d, line, hline, ribbon, labs, facet_wrap

Additional aesthetics and parameters for the geoms (see customization example).

Examples

# ==============================================
# 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: hline
#> data = [18 x 8]
#> yintercept = 0
#> linetype = 'dashed'
#> 
#> Layer 2
#> --------
#> Geom type: pointrange
#> data = [18 x 8]
#> aes_string(
#>   y = 'Coefficient'
#>   x = 'Level'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#> )
#> 
#> Layer 3
#> --------
#> Geom type: coord_flip
#> 
#> Layer 4
#> --------
#> Geom type: labs
#> title = 'Group-level Scores'
#> 
plot(layers)


# \donttest{
# 2 random intercepts
model <- lme4::lmer(Reaction ~ Days + (1 | Subject) + (1 | Newfactor), data = data)
x <- estimate_grouplevel(model)
plot(visualisation_recipe(x))


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

# }
# Simple Model ---------------
x <- estimate_means(lm(Sepal.Width ~ Species, data = iris))
#> We selected `by = c("Species")`.
layers <- visualisation_recipe(x)
layers
#> Layer 1
#> --------
#> Geom type: jitter
#> data = [150 x 2]
#> aes_string(
#>   x = 'Species'
#>   y = 'Sepal.Width'
#> )
#> stroke = 0
#> shape = 16
#> width = 0.1
#> height = 0
#> 
#> Layer 2
#> --------
#> Geom type: line
#> data = [3 x 5]
#> aes_string(
#>   y = 'Mean'
#>   x = 'Species'
#>   group = 1
#> )
#> 
#> Layer 3
#> --------
#> Geom type: pointrange
#> data = [3 x 5]
#> aes_string(
#>   y = 'Mean'
#>   x = 'Species'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#> )
#> 
#> Layer 4
#> --------
#> Geom type: labs
#> x = 'Species'
#> y = 'Sepal.Width'
#> title = 'Estimated Means (Sepal.Width ~ Species)'
#> 
plot(layers)

# \dontrun{
# Customize aesthetics
layers <- visualisation_recipe(x,
  jitter = list(width = 0.03, color = "red"),
  line = list(linetype = "dashed")
)
plot(layers)


# Customize raw data
plot(visualisation_recipe(x, show_data = c("violin", "boxplot", "jitter")))


# Two levels ---------------
data <- mtcars
data$cyl <- as.factor(data$cyl)
data$new_factor <- as.factor(rep(c("A", "B"), length.out = nrow(mtcars)))

# Modulations --------------
x <- estimate_means(model, by = c("new_factor", "wt"))
#> Error: Couldn't find which variables were selected in `by`. Check spelling and
#>   specification.
plot(visualisation_recipe(x))


# x <- estimate_means(model, by =c("new_factor", "cyl", "wt"))
# plot(visualisation_recipe(x))  # TODO: broken

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

# }
# ==============================================
# estimate_relation, estimate_expectation, ...
# ==============================================
# Simple Model ---------------
x <- estimate_relation(lm(mpg ~ wt, data = mtcars))
layers <- visualisation_recipe(x)
layers
#> Layer 1
#> --------
#> Geom type: point
#> data = [32 x 2]
#> aes_string(
#>   x = 'wt'
#>   y = 'mpg'
#> )
#> stroke = 0
#> shape = 16
#> 
#> Layer 2
#> --------
#> Geom type: ribbon
#> data = [10 x 5]
#> aes_string(
#>   y = 'Predicted'
#>   x = 'wt'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#>   group = 1
#> )
#> alpha = 0.3333333
#> 
#> Layer 3
#> --------
#> Geom type: line
#> data = [10 x 5]
#> aes_string(
#>   y = 'Predicted'
#>   x = 'wt'
#>   group = 1
#> )
#> 
#> Layer 4
#> --------
#> Geom type: labs
#> x = 'wt'
#> y = 'mpg'
#> title = 'Predicted response (mpg ~ wt)'
#> 
plot(layers)


# \dontrun{
# Customize aesthetics ----------

layers <- visualisation_recipe(x,
  point = list(color = "red", alpha = 0.6, size = 3),
  line = list(color = "blue", size = 3),
  ribbon = list(fill = "green", alpha = 0.7),
  labs = list(subtitle = "Oh yeah!")
)
layers
#> Layer 1
#> --------
#> Geom type: point
#> data = [32 x 2]
#> aes_string(
#>   x = 'wt'
#>   y = 'mpg'
#> )
#> stroke = 0
#> shape = 16
#> color = 'red'
#> alpha = 0.6
#> size = 3
#> 
#> Layer 2
#> --------
#> Geom type: ribbon
#> data = [10 x 5]
#> aes_string(
#>   y = 'Predicted'
#>   x = 'wt'
#>   ymin = 'CI_low'
#>   ymax = 'CI_high'
#>   group = 1
#> )
#> alpha = 0.7
#> fill = 'green'
#> 
#> Layer 3
#> --------
#> Geom type: line
#> data = [10 x 5]
#> aes_string(
#>   y = 'Predicted'
#>   x = 'wt'
#>   group = 1
#> )
#> color = 'blue'
#> size = 3
#> 
#> Layer 4
#> --------
#> Geom type: labs
#> x = 'wt'
#> y = 'mpg'
#> title = 'Predicted response (mpg ~ wt)'
#> subtitle = 'Oh yeah!'
#> 
plot(layers)


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

plot(visualisation_recipe(x, show_data = "none"))

plot(visualisation_recipe(x, show_data = c("density_2d", "points")))

plot(visualisation_recipe(x, show_data = "density_2d_filled"))

plot(visualisation_recipe(x, show_data = "density_2d_polygon"))

plot(visualisation_recipe(x, show_data = "density_2d_raster")) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0))


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

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

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


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

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


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


# Factor * numeric
x <- estimate_relation(lm(Sepal.Width ~ Species * Sepal.Length, data = iris))
layers <- visualisation_recipe(x)
plot(layers)


# 3-ways interaction ------------

data <- mtcars
data$vs <- as.factor(data$vs)
data$cyl <- as.factor(data$cyl)
data$new_factor <- as.factor(rep(c("A", "B"), length.out = nrow(mtcars)))

# Numeric * numeric * factor
x <- estimate_relation(lm(mpg ~ wt * am * vs, data = data))
layers <- visualisation_recipe(x)
plot(layers)


# Numeric * factor * factor
x <- estimate_relation(lm(mpg ~ wt * cyl * new_factor, data = data))
layers <- visualisation_recipe(x)
plot(layers)


# Factor * numeric * numeric
x <- estimate_relation(lm(mpg ~ cyl * qsec * hp, data = data))
layers <- visualisation_recipe(x)
plot(layers) +
  scale_size_continuous(range = c(0.2, 1))


# GLMs ---------------------
x <- estimate_relation(glm(vs ~ mpg, data = mtcars, family = "binomial"))
plot(visualisation_recipe(x))

plot(visualisation_recipe(x, show_data = "jitter", point = list(height = 0.03)))


# Multiple CIs ---------------------
plot(estimate_relation(lm(mpg ~ disp, data = mtcars),
  ci = c(.50, .80, .95)
))

plot(estimate_relation(lm(Sepal.Length ~ Species, data = iris),
  ci = c(0.5, 0.7, 0.95)
))


# Bayesian models ---------------------
if (require("ggplot2") && require("rstanarm")) {
  model <- rstanarm::stan_glm(mpg ~ wt, data = mtcars, refresh = 0)

  # Plot individual draws instead of regular ribbon
  x <- estimate_relation(model, keep_iterations = 100)
  layers <- visualisation_recipe(x, ribbon = list(color = "red"))
  plot(layers)

  model <- rstanarm::stan_glm(Sepal.Width ~ Species * Sepal.Length, data = iris, refresh = 0)
  plot(estimate_relation(model, keep_iterations = 100))
}
#> Loading required package: rstanarm
#> Loading required package: Rcpp
#> This is rstanarm version 2.32.1
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())

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

  layers <- visualisation_recipe(x)
  layers
  plot(layers)

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

  model <- lm(Petal.Length ~ Species * poly(Sepal.Width, 3), data = iris)
  x <- estimate_slopes(model, by = c("Sepal.Width", "Species"))
  plot(visualisation_recipe(x))
}
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Sepal.Width"`.
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Sepal.Width"`.
#> Warning: Using alpha for a discrete variable is not advised.

# \donttest{
# TODO: fails with latest emmeans (1.8.0)
if (require("mgcv")) {
  data <- iris
  data$Petal.Length <- data$Petal.Length^2

  model <- mgcv::gam(Sepal.Width ~ t2(Petal.Width, Petal.Length), data = data)
  x <- estimate_slopes(model, by = c("Petal.Width", "Petal.Length"), length = 20)
  plot(visualisation_recipe(x))

  model <- mgcv::gam(Sepal.Width ~ t2(Petal.Width, Petal.Length, by = Species), data = data)
  x <- estimate_slopes(model, by = c("Petal.Width", "Petal.Length", "Species"), length = 10)
  plot(visualisation_recipe(x))
}
#> Loading required package: mgcv
#> Loading required package: nlme
#> 
#> Attaching package: ‘nlme’
#> The following object is masked from ‘package:lme4’:
#> 
#>     lmList
#> This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Petal.Width"`.
#> No numeric variable was specified for slope estimation. Selecting `trend
#>   = "Petal.Width"`.

# }