diff --git a/.Rbuildignore b/.Rbuildignore
index 63a213b..257333f 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -32,3 +32,5 @@
 ^vignettes/.+\.pdf$
 ^vignettes/.+\.tex$
 ^data-raw$
+^.*\.Rproj$
+^\.Rproj\.user$
diff --git a/.gitignore b/.gitignore
index 35edf0e..b4bd72c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,3 +20,4 @@ ignore/
 ignore/*
 revdep/
 revdep/*
+.Rproj.user
diff --git a/DESCRIPTION b/DESCRIPTION
index bfa2fcc..e98c463 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -18,7 +18,11 @@ Authors@R: c(person("Thomas J.", "Leeper",
              person("Jacob A.", "Long", 
                     role = c("ctb"),
                     email = "long.1377@osu.edu",
-                    comment = c(ORCID = "0000-0002-1582-6214"))
+                    comment = c(ORCID = "0000-0002-1582-6214")),
+             person("Grant R.", "McDermott", 
+                    role = c("ctb"),
+                    email = "grantmcd@uoregon.edu",
+                    comment = c(ORCID = "0000-0001-7883-8573"))
              )
 URL: https://github.com/leeper/margins
 BugReports: https://github.com/leeper/margins/issues
@@ -48,4 +52,4 @@ Enhances:
     survey
 ByteCompile: true
 VignetteBuilder: knitr
-RoxygenNote: 7.1.0
+RoxygenNote: 7.1.1
diff --git a/R/margins.R b/R/margins.R
index 2d7133a..2777871 100644
--- a/R/margins.R
+++ b/R/margins.R
@@ -19,6 +19,7 @@
 #' @param iterations If \code{vce = "bootstrap"}, the number of bootstrap iterations. If \code{vce = "simulation"}, the number of simulated effects to draw. Ignored otherwise.
 #' @param unit_ses If \code{vce = "delta"}, a logical specifying whether to calculate and return unit-specific marginal effect variances. This calculation is time consuming and the information is often not needed, so this is set to \code{FALSE} by default.
 #' @param eps A numeric value specifying the \dQuote{step} to use when calculating numerical derivatives.
+#' @param save Logical indicating whether the original model object itself should be saved as an attribute. (This may be useful for further processing of margins objects, e.g., as part of model summaries.)
 #' @param \dots Arguments passed to methods, and onward to \code{\link{dydx}} methods and possibly further to \code{\link[prediction]{prediction}} methods. This can be useful, for example, for setting \code{type} (predicted value type), \code{eps} (precision), or \code{category} (category for multi-category outcome models), etc.
 #' @details Methods for this generic return a \dQuote{margins} object, which is a data frame consisting of the original data, predicted values and standard errors thereof, estimated marginal effects from the model \code{model} (for all variables used in the model, or the subset specified by \code{variables}), along with attributes describing various features of the marginal effects estimates.
 #' 
diff --git a/R/margins_betareg.R b/R/margins_betareg.R
index 0ba371c..dd82bab 100644
--- a/R/margins_betareg.R
+++ b/R/margins_betareg.R
@@ -11,6 +11,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -63,5 +64,6 @@ function(model,
               vcov = vc,
               jacobian = jac,
               weighted = FALSE,
-              iterations = if (vce == "bootstrap") iterations else NULL)
+              iterations = if (vce == "bootstrap") iterations else NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_clm.R b/R/margins_clm.R
index 41e50f3..994bd45 100644
--- a/R/margins_clm.R
+++ b/R/margins_clm.R
@@ -8,6 +8,7 @@ function(model,
          type = c("response", "link"),
          vce = "none",
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -54,5 +55,6 @@ function(model,
               vcov = NULL,
               jacobian = NULL,
               weighted = FALSE,
-              iterations = NULL)
+              iterations = NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_default.R b/R/margins_default.R
index f922c98..71cee76 100644
--- a/R/margins_default.R
+++ b/R/margins_default.R
@@ -11,6 +11,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -59,6 +60,7 @@ function(model,
               type = type,
               call = if ("call" %in% names(model)) model[["call"]] else NULL,
               model_class = class(model),
+              object = if (save) model else NULL,
               vce = vce,
               vcov = vc,
               jacobian = jac,
diff --git a/R/margins_glm.R b/R/margins_glm.R
index 5cc3257..be6d26d 100644
--- a/R/margins_glm.R
+++ b/R/margins_glm.R
@@ -11,6 +11,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -67,5 +68,6 @@ function(model,
               vcov = vc,
               jacobian = jac,
               weighted = FALSE,
-              iterations = if (vce == "bootstrap") iterations else NULL)
+              iterations = if (vce == "bootstrap") iterations else NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_loess.R b/R/margins_loess.R
index 007310c..e9148e5 100644
--- a/R/margins_loess.R
+++ b/R/margins_loess.R
@@ -7,6 +7,7 @@ function(model,
          at = NULL, 
          vce = "none", 
          eps = 1e-7,
+         save = FALSE,
          ...){
     
     # setup data
@@ -54,5 +55,6 @@ function(model,
               vcov = NULL,
               jacobian = NULL,
               weighted = FALSE,
-              iterations = NULL)
+              iterations = NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_merMod.R b/R/margins_merMod.R
index 3ff60ea..3282b1a 100644
--- a/R/margins_merMod.R
+++ b/R/margins_merMod.R
@@ -11,6 +11,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -61,7 +62,8 @@ function(model,
               vcov = vc,
               jacobian = jac,
               weighted = FALSE,
-              iterations = if (vce == "bootstrap") iterations else NULL)
+              iterations = if (vce == "bootstrap") iterations else NULL,
+              object = if (save) model else NULL)
 }
 
 #' @rdname margins
diff --git a/R/margins_multinom.R b/R/margins_multinom.R
index d5abcf3..9dc1598 100644
--- a/R/margins_multinom.R
+++ b/R/margins_multinom.R
@@ -12,6 +12,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -74,5 +75,6 @@ function(model,
               vcov = vc,
               jacobian = jac,
               weighted = FALSE,
-              iterations = if (vce == "bootstrap") iterations else NULL)
+              iterations = if (vce == "bootstrap") iterations else NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_nnet.R b/R/margins_nnet.R
index b128f12..f3ee235 100644
--- a/R/margins_nnet.R
+++ b/R/margins_nnet.R
@@ -7,6 +7,7 @@ function(model,
          at = NULL, 
          vce = "none",
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -49,5 +50,6 @@ function(model,
               vcov = NULL,
               jacobian = NULL,
               weighted = FALSE,
-              iterations = NULL)
+              iterations = NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_polr.R b/R/margins_polr.R
index f5140c1..f5ff789 100644
--- a/R/margins_polr.R
+++ b/R/margins_polr.R
@@ -12,6 +12,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # match.arg()
@@ -74,5 +75,6 @@ function(model,
               vcov = vc,
               jacobian = jac,
               weighted = FALSE,
-              iterations = if (vce == "bootstrap") iterations else NULL)
+              iterations = if (vce == "bootstrap") iterations else NULL,
+              object = if (save) model else NULL)
 }
diff --git a/R/margins_svyglm.R b/R/margins_svyglm.R
index a310482..2b55016 100644
--- a/R/margins_svyglm.R
+++ b/R/margins_svyglm.R
@@ -12,6 +12,7 @@ function(model,
          iterations = 50L, # if vce == "bootstrap" or "simulation"
          unit_ses = FALSE,
          eps = 1e-7,
+         save = FALSE,
          ...) {
     
     # require survey
@@ -81,5 +82,6 @@ function(model,
               vcov = vc,
               jacobian = jac,
               weighted = if (is.null(wts)) FALSE else TRUE,
-              iterations = if (vce == "bootstrap") iterations else NULL)
+              iterations = if (vce == "bootstrap") iterations else NULL,
+              object = if (save) model else NULL)
 }
diff --git a/man/margins.Rd b/man/margins.Rd
index 9cde557..61cbf86 100644
--- a/man/margins.Rd
+++ b/man/margins.Rd
@@ -35,6 +35,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -46,6 +47,7 @@ margins(model, ...)
   type = c("response", "link"),
   vce = "none",
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -60,6 +62,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -74,6 +77,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -88,6 +92,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -98,6 +103,7 @@ margins(model, ...)
   at = NULL,
   vce = "none",
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -112,6 +118,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -126,6 +133,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -140,6 +148,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -150,6 +159,7 @@ margins(model, ...)
   at = NULL,
   vce = "none",
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -164,6 +174,7 @@ margins(model, ...)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 
@@ -181,6 +192,7 @@ margins_summary(model, ..., level = 0.95, by_factor = TRUE)
   iterations = 50L,
   unit_ses = FALSE,
   eps = 1e-07,
+  save = FALSE,
   ...
 )
 }
@@ -207,6 +219,8 @@ margins_summary(model, ..., level = 0.95, by_factor = TRUE)
 
 \item{eps}{A numeric value specifying the \dQuote{step} to use when calculating numerical derivatives.}
 
+\item{save}{Logical indicating whether the original model object itself should be saved as an attribute. (This may be useful for further processing of margins objects, e.g., as part of model summaries.)}
+
 \item{level}{A numeric value specifying the confidence level for calculating p-values and confidence intervals.}
 
 \item{by_factor}{A logical specifying whether to order the output by factor (the default, \code{TRUE}).}
diff --git a/tests/testthat/tests-core.R b/tests/testthat/tests-core.R
index ca7d300..280ebc8 100644
--- a/tests/testthat/tests-core.R
+++ b/tests/testthat/tests-core.R
@@ -96,6 +96,36 @@ test_that("confint() for 'margins' object", {
     expect_true(inherits(confint(m), "matrix"), label = "confint() for margins")
 })
 
+context("Saving model object")
+test_that("save = TRUE works", {
+    x <- lm(mpg ~ wt * hp, data = mtcars)
+    m <- margins(x, save = TRUE)
+    x_df <- with(summary(x),
+                 data.frame(
+                     r.squared = r.squared,
+                    adj.r.squared = adj.r.squared,
+                    sigma = sigma,
+                    logLik = as.numeric(stats::logLik(x)),
+                    AIC = stats::AIC(x),
+                    BIC = stats::BIC(x),
+                    deviance = stats::deviance(x),
+                    df.residual = df.residual(x),
+                    nobs = stats::nobs(x)
+                    ))
+    m_df <- with(summary(attr(m, 'object')),
+                 data.frame(
+                     r.squared = r.squared,
+                     adj.r.squared = adj.r.squared,
+                     sigma = sigma,
+                     logLik = as.numeric(stats::logLik(attr(m, 'object'))),
+                     AIC = stats::AIC(attr(m, 'object')),
+                     BIC = stats::BIC(attr(m, 'object')),
+                     deviance = stats::deviance(attr(m, 'object')),
+                     df.residual = df.residual(attr(m, 'object')),
+                     nobs = stats::nobs(attr(m, 'object'))
+                     ))
+    expect_equal(x_df, m_df, label = "Original lm object correctly saved as attribute")
+})
 
 context("Variance tests")
 test_that("minimum test of variance calculations", {