diff --git a/R/doeAnalysis.R b/R/doeAnalysis.R
index 06201c9a..bd47f4dd 100644
--- a/R/doeAnalysis.R
+++ b/R/doeAnalysis.R
@@ -15,6 +15,21 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 
+.formatCoefEquation <- function(x, digits = .numDecimals) {
+  vapply(x, function(val) {
+    if (val == 0) return("0")
+    absVal <- abs(val)
+    if (absVal < 1e4) {
+      # avoid e-notation for small-to-moderate numbers
+      mag <- floor(log10(absVal))
+      sigDigits <- max(digits, mag + 1 + digits)
+      trimws(formatC(val, digits = sigDigits, format = "g", drop0trailing = TRUE))
+    } else {
+      trimws(formatC(val, digits = digits, format = "g", drop0trailing = TRUE))
+    }
+  }, character(1))
+}
+
 #' @export
 doeAnalysis <- function(jaspResults, dataset, options, ...) {
 
@@ -604,27 +619,28 @@ doeAnalysis <- function(jaspResults, dataset, options, ...) {
     }
 
     # Uncoded Formula
-
+    fmtCoef <- .formatCoefEquation(abs(coefs[-1, 1]), digits = .numDecimals)
     coefFormula <- paste(ifelse(sign(coefs[-1,1])==1, " +", " \u2013"),
-                         round(abs(coefs[-1,1]), .numDecimals),
+                         fmtCoef,
                          coefNames[-1],
                          collapse = "", sep = " ")
 
     # Now add dependent name and intercept
     filledFormula <- paste0(dep, " = ",
-                            round(coefs[1, 1], .numDecimals),
+                            .formatCoefEquation(coefs[1, 1], digits = .numDecimals),
                             coefFormula)
 
     # Coded Formula
+    fmtCoefCoded <- .formatCoefEquation(abs(coefsCoded[-1, 1]), digits = .numDecimals)
     coefFormulaCoded <- paste(ifelse(sign(coefsCoded[-1,1])==1, " +", " \u2013"),
-                         round(abs(coefsCoded[-1,1]), .numDecimals),
-                         coefNames[-1],
-                         collapse = "", sep = " ")
+                              fmtCoefCoded,
+                              coefNames[-1],
+                              collapse = "", sep = " ")
 
 
     filledFormulaCoded <- paste0(dep, " = ",
-                            round(coefsCoded[1, 1], .numDecimals),
-                            coefFormulaCoded)
+                                 .formatCoefEquation(coefsCoded[1, 1], digits = .numDecimals),
+                                 coefFormulaCoded)
 
     result[["regression"]][["filledFormula"]] <- gsubInteractionSymbol(filledFormula)
     jaspResults[[currentDependent]][["doeResult"]] <- createJaspState(result)
@@ -770,6 +786,17 @@ get_levels <- function(var, num_levels, dataset) {
 
   tb2 <- createJaspTable(gettext("Response Optimizer Solution"))
   tb2$addColumnInfo(name = "desi", title = "Composite desirability", type = "number")
+  allPredictors <- c(continuousPredictors, unlist(discretePredictors))
+  allPredictors <- allPredictors[allPredictors != ""]
+  for (pred in allPredictors) {
+    predType <- if (pred %in% continuousPredictors) "number" else "string"
+    tb2$addColumnInfo(name = pred, title = pred, type = predType)
+  }
+  for (dep in roDependent) {
+    colName  <- paste0(dep, "_fit")
+    colTitle <- gettextf("%s fit", dep)
+    tb2$addColumnInfo(name = colName, title = colTitle, type = "number")
+  }
   tb2$dependOn(options = c("responsesResponseOptimizer", "responseOptimizerGoal", "responseOptimizerLowerBound", "responseOptimizerTarget",
                              "responseOptimizerUpperBound", "responseOptimizerWeight", "responseOptimizerImportance",
                              "optimizationPlotCustomParameters", "responseOptimizerManualBounds", "responseOptimizerManualTarget",
@@ -779,17 +806,16 @@ get_levels <- function(var, num_levels, dataset) {
   jaspResults[["tableRoSolution"]] <- tb2
 
   rows2 <- data.frame(desi = desi)
-  if (length(continuousPredictors) >= 1 && length(optimParam) > 1 && !identical(continuousPredictors, "")) {
-    optimParam[continuousPredictors] <- round(optimParam[continuousPredictors], .numDecimals)
-  } else if (length(continuousPredictors) == 1 && length(optimParam) == 1 && !identical(continuousPredictors, "")) {
-    optimParam <- round(optimParam, .numDecimals)
+  for (pred in allPredictors) {
+    if (pred %in% continuousPredictors)
+      rows2[[pred]] <- round(as.numeric(optimParam[[pred]]), .numDecimals)
+    else
+      rows2[[pred]] <- as.character(optimParam[[pred]])
+  }
+  for (dep in roDependent) {
+    colName <- paste0(dep, "_fit")
+    rows2[[colName]] <- round(as.numeric(predValues[[dep]]), .numDecimals)
   }
-  rows2 <- cbind(rows2, optimParam)
-  if (length(optimParam) == 1) # transform it this way to have the correct naming
-    colnames(rows2)[ncol(rows2)] <- c(continuousPredictors, discretePredictors)[c(continuousPredictors, discretePredictors) != ""]
-  predValuesFit <- setNames(predValues, paste0(names(predValues), " fit"))
-  predValuesFit <- round(predValuesFit, .numDecimals)
-  rows2 <- cbind(rows2, as.data.frame(t(predValuesFit)))
   tb2$addRows(rows2)
 }
 
@@ -1017,7 +1043,18 @@ get_levels <- function(var, num_levels, dataset) {
 
   tbTitle <- if (options[["optimizationPlotCustomParameters"]]) gettext("Optimization Plot Summary (Manual Predictor Values)") else gettext("Optimization Plot Summary")
   tb <- createJaspTable(tbTitle)
-  tb$addColumnInfo(name = "desi", title = "Composite desirability", type = "number")
+  tb$addColumnInfo(name = "desi", title = gettext("Composite desirability"), type = "number")
+  for (pred in allPredictors) {
+    predType <- if (pred %in% continuousPredictors) "number" else "string"
+    tb$addColumnInfo(name = pred, title = pred, type = predType)
+  }
+
+  for (dep in roDependent) {
+    colName  <- paste0(dep, " fit")
+    colTitle <- gettextf("%s fit", dep)
+    tb$addColumnInfo(name = colName, title = colTitle, type = "number")
+  }
+
   tb$dependOn(options = c("responsesResponseOptimizer", "responseOptimizerGoal", "responseOptimizerLowerBound", "responseOptimizerTarget",
                     "responseOptimizerUpperBound", "responseOptimizerWeight", "responseOptimizerImportance",
                     "optimizationPlotCustomParameters", "responseOptimizerManualBounds", "responseOptimizerManualTarget",
@@ -1034,17 +1071,14 @@ get_levels <- function(var, num_levels, dataset) {
                               currentDiscreteLevels = currentDiscreteLevels, coefficients = coefficients, dependent = roDependent)
 
   rows <- data.frame(desi = desi)
-  predValuesFit <- setNames(predValues, paste0(names(predValues), " fit"))
-  predValuesFit <- round(predValuesFit, .numDecimals)
-  rows <- cbind(rows, as.data.frame(t(predValuesFit)))
-  if (length(continuousPredictors) >= 1 && length(currentSettings) > 1 && !identical(continuousPredictors, "")) {
-    currentSettings[continuousPredictors] <- round(currentSettings[continuousPredictors], .numDecimals)
-  } else if (length(continuousPredictors) == 1 && length(currentSettings) == 1 && !identical(continuousPredictors, "")) {
-    currentSettings <- round(currentSettings, .numDecimals)
+  for (pred in allPredictors) {
+    if (pred %in% continuousPredictors)
+      rows[[pred]] <- round(as.numeric(currentSettings[[pred]]), .numDecimals)
+    else
+      rows[[pred]] <- as.character(currentSettings[[pred]])
   }
-  rows <- cbind(rows, currentSettings)
-  if (length(currentSettings) == 1) # transform it this way to have the correct naming
-    colnames(rows)[ncol(rows)] <- allPredictors
+  for (dep in roDependent)
+    rows[[paste0(dep, " fit")]] <- round(as.numeric(predValues[[dep]]), .numDecimals)
   tb$addRows(rows)
 }