diff --git a/NAMESPACE b/NAMESPACE
index 4de9f05ab..7f8a759d8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -100,6 +100,7 @@ S3method(nmObjGet,methodOde)
 S3method(nmObjGet,modelName)
 S3method(nmObjGet,notes)
 S3method(nmObjGet,omegaR)
+S3method(nmObjGet,parHist)
 S3method(nmObjGet,parHistStacked)
 S3method(nmObjGet,phiR)
 S3method(nmObjGet,phiRSE)
diff --git a/NEWS.md b/NEWS.md
index c796d8181..4f1ba20b7 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -12,6 +12,8 @@
 - Values in `$parFixed` for BSV without exponential transformation are now
   correctly shown (#366)
 
+- Add back values for $parHistData (#368)
+
 # nlmixr2est 2.1.6
 
 ## Breaking changes
diff --git a/R/complete.R b/R/complete.R
index 8600540df..324e57636 100644
--- a/R/complete.R
+++ b/R/complete.R
@@ -9,9 +9,10 @@
   objDf="Objective Function DF",
   omega="Omega Matrix",
   origData="Original Data",
+  phiC="covariance matrix of each individual's eta (if present)",
   parFixed="Formatted Parameter Values for Fixed effects",
   parFixedDf="Parameter Values for Fixed Effects (data frame)",
-  parHist="Parameter History",
+  parHistData="Parameter History (including gradients)",
   scaleInfo="Scaling Information",
   shrink="Shrinkage data frame",
   table="Table Control Value",
diff --git a/R/focei.R b/R/focei.R
index d7c21c015..9c3e39b42 100644
--- a/R/focei.R
+++ b/R/focei.R
@@ -1547,6 +1547,19 @@ attr(rxUiGet.foceiOptEnv, "desc") <- "Get focei optimization environment"
          deparse1(ui$dvidLine)),
         collapse="\n")
 }
+#' Calculate the parameter history
+#'
+#' @param .ret return data
+#' @return parameter history data frame
+#' @noRd
+#' @author Matthew L. Fidler
+.parHistCalc <- function(.ret) {
+  .tmp <- .ret$parHistData
+  .tmp <- .tmp[.tmp$type == "Unscaled", names(.tmp) != "type"]
+  .iter <- .tmp$iter
+  .tmp <- .tmp[, names(.tmp) != "iter"]
+  data.frame(iter = .iter, .tmp, check.names=FALSE)
+}
 
 #' Setup the par history information
 #'
@@ -1556,13 +1569,7 @@ attr(rxUiGet.foceiOptEnv, "desc") <- "Get focei optimization environment"
 #' @noRd
 .foceiSetupParHistData <- function(.ret) {
   if (exists("parHistData", envir=.ret)) {
-    .tmp <- .ret$parHistData
-    .tmp <- .tmp[.tmp$type == "Unscaled", names(.tmp) != "type"]
-    .iter <- .tmp$iter
-    .tmp <- .tmp[, names(.tmp) != "iter"]
-    ## .ret$parHistStacked <- data.frame(stack(.tmp), iter = .iter)
-    ## names(.ret$parHistStacked) <- c("val", "par", "iter")
-    .ret$parHist <- data.frame(iter = .iter, .tmp)
+    .ret$parHist <- .parHistCalc(.ret)
   }
 }
 
@@ -1637,7 +1644,7 @@ attr(rxUiGet.foceiOptEnv, "desc") <- "Get focei optimization environment"
       .env$saem0 <- .saem
     }
     if (.control$compress) {
-      for (.item in c("origData", "phiM", "parHist", "saem0")) {
+      for (.item in c("origData", "phiM", "parHistData", "saem0")) {
         if (exists(.item, .env)) {
           .obj <- get(.item, envir=.env)
           .size <- utils::object.size(.obj)
@@ -1659,7 +1666,7 @@ attr(rxUiGet.foceiOptEnv, "desc") <- "Get focei optimization environment"
                     "xType", "IDlabel", "ODEmodel",
                     # times
                     "optimTime", "setupTime", "covTime",
-                    "parHistData", "dataSav", "idLvl", "theta",
+                    "parHist", "dataSav", "idLvl", "theta",
                     "missingTable", "missingControl", "missingEst")) {
       if (exists(.item, .env)) {
         rm(list=.item, envir=.env)
@@ -1913,4 +1920,3 @@ nlmixr2CreateOutputFromUi <- function(ui, data=NULL, control=NULL, table=NULL, e
   class(.env) <- c("output", "nlmixr2Est")
   nlmixr2Est(.env)
 }
-
diff --git a/R/nmObjGet.R b/R/nmObjGet.R
index 8e3db96df..d093c3db7 100644
--- a/R/nmObjGet.R
+++ b/R/nmObjGet.R
@@ -488,22 +488,32 @@ nmObjGet.saemTransformedData <- function(x, ...) {
 }
 #attr(nmObjGet.saemTransformedData, "desc") <- "data that saem sees for optimization"
 
+#' @rdname nmObjGet
+#' @export
+nmObjGet.parHist <- function(x, ...) {
+  .obj <- x[[1]]
+  .env <- .obj$env
+  if (exists("parHistData", envir=.env)) {
+    return(.parHistCalc(.env))
+  }
+  NULL
+}
+attr(nmObjGet.parHist, "desc") <- "Parameter History"
 
 #' @rdname nmObjGet
 #' @export
 nmObjGet.parHistStacked <- function(x, ...) {
   .obj <- x[[1]]
   .env <- .obj$env
-  if (exists("parHist", envir=.env)) {
-    .parHist <- .env$parHist
+  if (exists("parHistData", envir=.env)) {
+    .parHist <- .parHistCalc(.env)
     .iter <- .parHist$iter
     .ret <- data.frame(iter=.iter, stack(.parHist[, -1]))
     names(.ret) <- sub("values", "val",
                        sub("ind", "par", names(.ret)))
-    .ret
-  } else {
-    NULL
+    return(.ret)
   }
+  NULL
 }
 attr(nmObjGet.parHistStacked, "desc") <- "stacked parameter history"
 
diff --git a/R/saem.R b/R/saem.R
index 9b277e35b..6b2aa7f92 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -396,12 +396,13 @@
   if (ncol(.m) > length(.allThetaNames)) {
     .m <- .m[, seq_along(.allThetaNames)]
   }
-  .ph <- data.frame(iter = rep(1:nrow(.m)), as.data.frame(.m))
-  names(.ph) <- c("iter", .allThetaNames)
+  .ph <- data.frame(iter = rep(seq_len(nrow(.m))), as.data.frame(.m),
+                    type="Unscaled", check.names=FALSE)
+  names(.ph) <- c("iter", .allThetaNames, "type")
   .cls <- class(.ph)
   attr(.cls, "niter") <- env$saemControl$mcmc$niter[1]
   class(.ph) <- .cls
-  assign("parHist", .ph, envir=env)
+  assign("parHistData", .ph, envir=env)
 }
 #' Calculate the covariance term
 #'
diff --git a/man/nmObjGet.Rd b/man/nmObjGet.Rd
index b58a07ba6..a469e6b0e 100644
--- a/man/nmObjGet.Rd
+++ b/man/nmObjGet.Rd
@@ -13,6 +13,7 @@
 \alias{nmObjGet.idLvl}
 \alias{nmObjGet.covLvl}
 \alias{nmObjGet.saemTransformedData}
+\alias{nmObjGet.parHist}
 \alias{nmObjGet.parHistStacked}
 \alias{nmObjGet.md5}
 \alias{nmObjGet.notes}
@@ -66,6 +67,8 @@ nmObjGet(x, ...)
 
 \method{nmObjGet}{saemTransformedData}(x, ...)
 
+\method{nmObjGet}{parHist}(x, ...)
+
 \method{nmObjGet}{parHistStacked}(x, ...)
 
 \method{nmObjGet}{md5}(x, ...)
diff --git a/tests/testthat/test-addCwres.R b/tests/testthat/test-addCwres.R
index 3ed4c87e8..0c408ffde 100644
--- a/tests/testthat/test-addCwres.R
+++ b/tests/testthat/test-addCwres.R
@@ -1,6 +1,7 @@
 nmTest({
 
   test_that("addCwres", {
+
     one.compartment <- function() {
       ini({
         tka <- log(1.57)
@@ -22,6 +23,7 @@ nmTest({
     suppressMessages(
       fitNoEta <- nlmixr2(one.compartment, theo_sd,  est="focei", control = list(print=0))
     )
+    expect_true(inherits(fitNoEta$parHistData, "data.frame"))
     expect_error(
       addCwres(fitNoEta),
       regexp = "cannot add CWRES to a model without etas"