diff --git a/DESCRIPTION b/DESCRIPTION
index 585a51285a..30689a9307 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -44,6 +44,7 @@ Imports:
     vctrs (>= 0.6.0),
     withr (>= 2.5.0)
 Suggests:
+    broom,
     covr,
     dplyr,
     ggplot2movies,
@@ -126,9 +127,8 @@ Collate:
     'facet-grid-.R'
     'facet-null.R'
     'facet-wrap.R'
-    'fortify-lm.R'
     'fortify-map.R'
-    'fortify-multcomp.R'
+    'fortify-models.R'
     'fortify-spatial.R'
     'fortify.R'
     'stat-.R'
diff --git a/NEWS.md b/NEWS.md
index 0c493a8f58..42b1bfaa98 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,5 +1,9 @@
 # ggplot2 (development version)
 
+* The following methods have been deprecated: `fortify.lm()`, `fortify.glht()`,
+  `fortify.confint.glht()`, `fortify.summary.glht()` and `fortify.cld()`. It
+  is recommend to use `broom::augment()` and `broom::tidy()` instead 
+  (@teunbrand, #3816).
 * Custom and raster annotation now respond to scale transformations, and can
   use AsIs variables for relative placement (@teunbrand based on 
   @yutannihilation's prior work, #3120)
diff --git a/R/fortify-lm.R b/R/fortify-lm.R
deleted file mode 100644
index ebfb0334b2..0000000000
--- a/R/fortify-lm.R
+++ /dev/null
@@ -1,87 +0,0 @@
-#' Supplement the data fitted to a linear model with model fit statistics.
-#'
-#' If you have missing values in your model data, you may need to refit
-#' the model with `na.action = na.exclude`.
-#'
-#' @return The original data with extra columns:
-#'   \item{.hat}{Diagonal of the hat matrix}
-#'   \item{.sigma}{Estimate of residual standard deviation when
-#'     corresponding observation is dropped from model}
-#'   \item{.cooksd}{Cooks distance, [cooks.distance()]}
-#'   \item{.fitted}{Fitted values of model}
-#'   \item{.resid}{Residuals}
-#'   \item{.stdresid}{Standardised residuals}
-#' @param model linear model
-#' @param data data set, defaults to data used to fit model
-#' @param ... not used by this method
-#' @keywords internal
-#' @export
-#' @examples
-#' mod <- lm(mpg ~ wt, data = mtcars)
-#' head(fortify(mod))
-#' head(fortify(mod, mtcars))
-#'
-#' plot(mod, which = 1)
-#'
-#' ggplot(mod, aes(.fitted, .resid)) +
-#'   geom_point() +
-#'   geom_hline(yintercept = 0) +
-#'   geom_smooth(se = FALSE)
-#'
-#' ggplot(mod, aes(.fitted, .stdresid)) +
-#'   geom_point() +
-#'   geom_hline(yintercept = 0) +
-#'   geom_smooth(se = FALSE)
-#'
-#' ggplot(fortify(mod, mtcars), aes(.fitted, .stdresid)) +
-#'   geom_point(aes(colour = factor(cyl)))
-#'
-#' ggplot(fortify(mod, mtcars), aes(mpg, .stdresid)) +
-#'   geom_point(aes(colour = factor(cyl)))
-#'
-#' plot(mod, which = 2)
-#' ggplot(mod) +
-#'   stat_qq(aes(sample = .stdresid)) +
-#'   geom_abline()
-#'
-#' plot(mod, which = 3)
-#' ggplot(mod, aes(.fitted, sqrt(abs(.stdresid)))) +
-#'   geom_point() +
-#'   geom_smooth(se = FALSE)
-#'
-#' plot(mod, which = 4)
-#' ggplot(mod, aes(seq_along(.cooksd), .cooksd)) +
-#'   geom_col()
-#'
-#' plot(mod, which = 5)
-#' ggplot(mod, aes(.hat, .stdresid)) +
-#'   geom_vline(linewidth = 2, colour = "white", xintercept = 0) +
-#'   geom_hline(linewidth = 2, colour = "white", yintercept = 0) +
-#'   geom_point() + geom_smooth(se = FALSE)
-#'
-#' ggplot(mod, aes(.hat, .stdresid)) +
-#'   geom_point(aes(size = .cooksd)) +
-#'   geom_smooth(se = FALSE, linewidth = 0.5)
-#'
-#' plot(mod, which = 6)
-#' ggplot(mod, aes(.hat, .cooksd)) +
-#'   geom_vline(xintercept = 0, colour = NA) +
-#'   geom_abline(slope = seq(0, 3, by = 0.5), colour = "white") +
-#'   geom_smooth(se = FALSE) +
-#'   geom_point()
-#'
-#' ggplot(mod, aes(.hat, .cooksd)) +
-#'   geom_point(aes(size = .cooksd / .hat)) +
-#'   scale_size_area()
-fortify.lm <- function(model, data = model$model, ...) {
-  infl <- stats::influence(model, do.coef = FALSE)
-  data$.hat <- infl$hat
-  data$.sigma <- infl$sigma
-  data$.cooksd <- stats::cooks.distance(model, infl)
-
-  data$.fitted <- stats::predict(model)
-  data$.resid <- stats::resid(model)
-  data$.stdresid <- stats::rstandard(model, infl)
-
-  data
-}
diff --git a/R/fortify-models.R b/R/fortify-models.R
new file mode 100644
index 0000000000..5a0e95199a
--- /dev/null
+++ b/R/fortify-models.R
@@ -0,0 +1,169 @@
+#' Supplement the data fitted to a linear model with model fit statistics.
+#'
+#' @description
+#' `r lifecycle::badge("deprecated")`
+#'
+#' This method is deprecated because using `broom::augment()` is a better
+#' solution to supplement data from a linear model.
+#' If you have missing values in your model data, you may need to refit
+#' the model with `na.action = na.exclude`.
+#'
+#' @return The original data with extra columns:
+#'   \item{.hat}{Diagonal of the hat matrix}
+#'   \item{.sigma}{Estimate of residual standard deviation when
+#'     corresponding observation is dropped from model}
+#'   \item{.cooksd}{Cooks distance, [cooks.distance()]}
+#'   \item{.fitted}{Fitted values of model}
+#'   \item{.resid}{Residuals}
+#'   \item{.stdresid}{Standardised residuals}
+#' @param model linear model
+#' @param data data set, defaults to data used to fit model
+#' @param ... not used by this method
+#' @keywords internal
+#' @export
+#' @examplesIf require("broom")
+#' mod <- lm(mpg ~ wt, data = mtcars)
+#'
+#' # Show augmented model
+#' head(augment(mod))
+#' head(fortify(mod))
+#'
+#' # Using augment to convert model to ready-to-plot data
+#' ggplot(augment(mod), aes(.fitted, .resid)) +
+#'   geom_point() +
+#'   geom_hline(yintercept = 0) +
+#'   geom_smooth(se = FALSE)
+#'
+#' # Colouring by original data not included in the model
+#' ggplot(augment(mod, mtcars), aes(.fitted, .std.resid, colour = factor(cyl))) +
+#'   geom_point()
+fortify.lm <- function(model, data = model$model, ...) {
+  lifecycle::deprecate_warn(
+    "3.6.0", I("`fortify(<lm>)`"), I("`broom::augment(<lm>)`")
+  )
+  infl <- stats::influence(model, do.coef = FALSE)
+  data$.hat <- infl$hat
+  data$.sigma <- infl$sigma
+  data$.cooksd <- stats::cooks.distance(model, infl)
+
+  data$.fitted <- stats::predict(model)
+  data$.resid <- stats::resid(model)
+  data$.stdresid <- stats::rstandard(model, infl)
+
+  data
+}
+
+#' Fortify methods for objects produced by \pkg{multcomp}
+#'
+#' @description
+#' `r lifecycle::badge("deprecated")`
+#'
+#' This function is deprecated because using `broom::tidy()` is a better
+#' solution to convert model objects.
+#'
+#' @param model an object of class `glht`, `confint.glht`,
+#'  `summary.glht` or [multcomp::cld()]
+#' @param data,... other arguments to the generic ignored in this method.
+#' @name fortify-multcomp
+#' @keywords internal
+#' @examplesIf require("multcomp") && require("broom")
+#' amod <- aov(breaks ~ wool + tension, data = warpbreaks)
+#' wht <- multcomp::glht(amod, linfct = multcomp::mcp(tension = "Tukey"))
+#'
+#' tidy(wht) # recommended
+#' fortify(wht)
+#'
+#' ggplot(tidy(wht), aes(contrast, estimate)) + geom_point()
+#'
+#' ci <- confint(wht)
+#' tidy(ci) # recommended
+#' fortify(ci)
+#'
+#' ggplot(tidy(confint(wht)),
+#'        aes(contrast, estimate, ymin = conf.low, ymax = conf.high)) +
+#'    geom_pointrange()
+#'
+#' smry <- summary(wht)
+#' tidy(smry) # recommended
+#' fortify(smry)
+#'
+#' ggplot(mapping = aes(contrast, estimate)) +
+#'    geom_linerange(aes(ymin = conf.low, ymax = conf.high), data = tidy(ci)) +
+#'    geom_point(aes(size = adj.p.value), data = tidy(smry)) +
+#'    scale_size(transform = "reverse")
+#'
+#' cld <- multcomp::cld(wht)
+#' tidy(cld) # recommended
+#' fortify(cld)
+NULL
+
+#' @method fortify glht
+#' @rdname fortify-multcomp
+#' @export
+fortify.glht <- function(model, data, ...) {
+  lifecycle::deprecate_warn(
+    "3.6.0", I("`fortify(<glht>)`"), I("`broom::tidy(<glht>)`")
+  )
+  base::data.frame(
+    lhs = rownames(model$linfct),
+    rhs = model$rhs,
+    estimate = stats::coef(model),
+    check.names = FALSE,
+    stringsAsFactors = FALSE
+  )
+}
+
+#' @rdname fortify-multcomp
+#' @method fortify confint.glht
+#' @export
+fortify.confint.glht <- function(model, data, ...) {
+  lifecycle::deprecate_warn(
+    "3.6.0", I("`fortify(<confint.glht>)`"), I("`broom::tidy(<confint.glht>)`")
+  )
+  coef <- model$confint
+  colnames(coef) <- to_lower_ascii(colnames(coef))
+
+  base::data.frame(
+    lhs = rownames(coef),
+    rhs = model$rhs,
+    coef,
+    check.names = FALSE,
+    stringsAsFactors = FALSE
+  )
+}
+
+#' @method fortify summary.glht
+#' @rdname fortify-multcomp
+#' @export
+fortify.summary.glht <- function(model, data, ...) {
+  lifecycle::deprecate_warn(
+    "3.6.0", I("`fortify(<summary.glht>)`"), I("`broom::tidy(<summary.glht>)`")
+  )
+  coef <- as.data.frame(
+    model$test[c("coefficients", "sigma", "tstat", "pvalues")])
+  names(coef) <- c("estimate", "se", "t", "p")
+
+  base::data.frame(
+    lhs = rownames(coef),
+    rhs = model$rhs,
+    coef,
+    check.names = FALSE,
+    stringsAsFactors = FALSE
+  )
+}
+
+
+#' @method fortify cld
+#' @rdname fortify-multcomp
+#' @export
+fortify.cld <- function(model, data, ...) {
+  lifecycle::deprecate_warn(
+    "3.6.0", I("`fortify(<summary.glht>)`"), I("`broom::tidy(<summary.glht>)`")
+  )
+  base::data.frame(
+    lhs = names(model$mcletters$Letters),
+    letters = model$mcletters$Letters,
+    check.names = FALSE,
+    stringsAsFactors = FALSE
+  )
+}
diff --git a/R/fortify-multcomp.R b/R/fortify-multcomp.R
deleted file mode 100644
index 79714b2a68..0000000000
--- a/R/fortify-multcomp.R
+++ /dev/null
@@ -1,89 +0,0 @@
-#' Fortify methods for objects produced by \pkg{multcomp}
-#'
-#' @param model an object of class `glht`, `confint.glht`,
-#'  `summary.glht` or [multcomp::cld()]
-#' @param data,... other arguments to the generic ignored in this method.
-#' @name fortify-multcomp
-#' @keywords internal
-#' @examples
-#' if (require("multcomp")) {
-#' amod <- aov(breaks ~ wool + tension, data = warpbreaks)
-#' wht <- glht(amod, linfct = mcp(tension = "Tukey"))
-#'
-#' fortify(wht)
-#' ggplot(wht, aes(lhs, estimate)) + geom_point()
-#'
-#' CI <- confint(wht)
-#' fortify(CI)
-#' ggplot(CI, aes(lhs, estimate, ymin = lwr, ymax = upr)) +
-#'    geom_pointrange()
-#'
-#' fortify(summary(wht))
-#' ggplot(mapping = aes(lhs, estimate)) +
-#'    geom_linerange(aes(ymin = lwr, ymax = upr), data = CI) +
-#'    geom_point(aes(size = p), data = summary(wht)) +
-#'    scale_size(transform = "reverse")
-#'
-#' cld <- cld(wht)
-#' fortify(cld)
-#' }
-NULL
-
-#' @method fortify glht
-#' @rdname fortify-multcomp
-#' @export
-fortify.glht <- function(model, data, ...) {
-  base::data.frame(
-    lhs = rownames(model$linfct),
-    rhs = model$rhs,
-    estimate = stats::coef(model),
-    check.names = FALSE,
-    stringsAsFactors = FALSE
-  )
-}
-
-#' @rdname fortify-multcomp
-#' @method fortify confint.glht
-#' @export
-fortify.confint.glht <- function(model, data, ...) {
-  coef <- model$confint
-  colnames(coef) <- to_lower_ascii(colnames(coef))
-
-  base::data.frame(
-    lhs = rownames(coef),
-    rhs = model$rhs,
-    coef,
-    check.names = FALSE,
-    stringsAsFactors = FALSE
-  )
-}
-
-#' @method fortify summary.glht
-#' @rdname fortify-multcomp
-#' @export
-fortify.summary.glht <- function(model, data, ...) {
-  coef <- as.data.frame(
-    model$test[c("coefficients", "sigma", "tstat", "pvalues")])
-  names(coef) <- c("estimate", "se", "t", "p")
-
-  base::data.frame(
-    lhs = rownames(coef),
-    rhs = model$rhs,
-    coef,
-    check.names = FALSE,
-    stringsAsFactors = FALSE
-  )
-}
-
-
-#' @method fortify cld
-#' @rdname fortify-multcomp
-#' @export
-fortify.cld <- function(model, data, ...) {
-  base::data.frame(
-    lhs = names(model$mcletters$Letters),
-    letters = model$mcletters$Letters,
-    check.names = FALSE,
-    stringsAsFactors = FALSE
-  )
-}
diff --git a/man/fortify-multcomp.Rd b/man/fortify-multcomp.Rd
index a52dec001c..654e5bbe9a 100644
--- a/man/fortify-multcomp.Rd
+++ b/man/fortify-multcomp.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fortify-multcomp.R
+% Please edit documentation in R/fortify-models.R
 \name{fortify-multcomp}
 \alias{fortify-multcomp}
 \alias{fortify.glht}
@@ -23,29 +23,41 @@
 \item{data, ...}{other arguments to the generic ignored in this method.}
 }
 \description{
-Fortify methods for objects produced by \pkg{multcomp}
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}}
+
+This function is deprecated because using \code{broom::tidy()} is a better
+solution to convert model objects.
 }
 \examples{
-if (require("multcomp")) {
+\dontshow{if (require("multcomp") && require("broom")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
 amod <- aov(breaks ~ wool + tension, data = warpbreaks)
-wht <- glht(amod, linfct = mcp(tension = "Tukey"))
+wht <- multcomp::glht(amod, linfct = multcomp::mcp(tension = "Tukey"))
 
+tidy(wht) # recommended
 fortify(wht)
-ggplot(wht, aes(lhs, estimate)) + geom_point()
 
-CI <- confint(wht)
-fortify(CI)
-ggplot(CI, aes(lhs, estimate, ymin = lwr, ymax = upr)) +
+ggplot(tidy(wht), aes(contrast, estimate)) + geom_point()
+
+ci <- confint(wht)
+tidy(ci) # recommended
+fortify(ci)
+
+ggplot(tidy(confint(wht)),
+       aes(contrast, estimate, ymin = conf.low, ymax = conf.high)) +
    geom_pointrange()
 
-fortify(summary(wht))
-ggplot(mapping = aes(lhs, estimate)) +
-   geom_linerange(aes(ymin = lwr, ymax = upr), data = CI) +
-   geom_point(aes(size = p), data = summary(wht)) +
+smry <- summary(wht)
+tidy(smry) # recommended
+fortify(smry)
+
+ggplot(mapping = aes(contrast, estimate)) +
+   geom_linerange(aes(ymin = conf.low, ymax = conf.high), data = tidy(ci)) +
+   geom_point(aes(size = adj.p.value), data = tidy(smry)) +
    scale_size(transform = "reverse")
 
-cld <- cld(wht)
+cld <- multcomp::cld(wht)
+tidy(cld) # recommended
 fortify(cld)
-}
+\dontshow{\}) # examplesIf}
 }
 \keyword{internal}
diff --git a/man/fortify.lm.Rd b/man/fortify.lm.Rd
index 4a994a6c56..d98b28a07f 100644
--- a/man/fortify.lm.Rd
+++ b/man/fortify.lm.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fortify-lm.R
+% Please edit documentation in R/fortify-models.R
 \name{fortify.lm}
 \alias{fortify.lm}
 \title{Supplement the data fitted to a linear model with model fit statistics.}
@@ -24,65 +24,30 @@ corresponding observation is dropped from model}
 \item{.stdresid}{Standardised residuals}
 }
 \description{
+\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}}
+
+This method is deprecated because using \code{broom::augment()} is a better
+solution to supplement data from a linear model.
 If you have missing values in your model data, you may need to refit
 the model with \code{na.action = na.exclude}.
 }
 \examples{
+\dontshow{if (require("broom")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf}
 mod <- lm(mpg ~ wt, data = mtcars)
-head(fortify(mod))
-head(fortify(mod, mtcars))
-
-plot(mod, which = 1)
 
-ggplot(mod, aes(.fitted, .resid)) +
-  geom_point() +
-  geom_hline(yintercept = 0) +
-  geom_smooth(se = FALSE)
+# Show augmented model
+head(augment(mod))
+head(fortify(mod))
 
-ggplot(mod, aes(.fitted, .stdresid)) +
+# Using augment to convert model to ready-to-plot data
+ggplot(augment(mod), aes(.fitted, .resid)) +
   geom_point() +
   geom_hline(yintercept = 0) +
   geom_smooth(se = FALSE)
 
-ggplot(fortify(mod, mtcars), aes(.fitted, .stdresid)) +
-  geom_point(aes(colour = factor(cyl)))
-
-ggplot(fortify(mod, mtcars), aes(mpg, .stdresid)) +
-  geom_point(aes(colour = factor(cyl)))
-
-plot(mod, which = 2)
-ggplot(mod) +
-  stat_qq(aes(sample = .stdresid)) +
-  geom_abline()
-
-plot(mod, which = 3)
-ggplot(mod, aes(.fitted, sqrt(abs(.stdresid)))) +
-  geom_point() +
-  geom_smooth(se = FALSE)
-
-plot(mod, which = 4)
-ggplot(mod, aes(seq_along(.cooksd), .cooksd)) +
-  geom_col()
-
-plot(mod, which = 5)
-ggplot(mod, aes(.hat, .stdresid)) +
-  geom_vline(linewidth = 2, colour = "white", xintercept = 0) +
-  geom_hline(linewidth = 2, colour = "white", yintercept = 0) +
-  geom_point() + geom_smooth(se = FALSE)
-
-ggplot(mod, aes(.hat, .stdresid)) +
-  geom_point(aes(size = .cooksd)) +
-  geom_smooth(se = FALSE, linewidth = 0.5)
-
-plot(mod, which = 6)
-ggplot(mod, aes(.hat, .cooksd)) +
-  geom_vline(xintercept = 0, colour = NA) +
-  geom_abline(slope = seq(0, 3, by = 0.5), colour = "white") +
-  geom_smooth(se = FALSE) +
+# Colouring by original data not included in the model
+ggplot(augment(mod, mtcars), aes(.fitted, .std.resid, colour = factor(cyl))) +
   geom_point()
-
-ggplot(mod, aes(.hat, .cooksd)) +
-  geom_point(aes(size = .cooksd / .hat)) +
-  scale_size_area()
+\dontshow{\}) # examplesIf}
 }
 \keyword{internal}