diff --git a/NAMESPACE b/NAMESPACE
index 77a34ecb..6e0d34fe 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method("[",hstats_matrix)
+S3method("dimnames<-",hstats_matrix)
 S3method(average_loss,Learner)
 S3method(average_loss,default)
 S3method(average_loss,explainer)
diff --git a/NEWS.md b/NEWS.md
index 5240d604..a8968029 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -5,7 +5,7 @@
 - `hstats()`: `n_max` has been increased from 300 to 500 rows. This will make estimates of H statistics more stable at the price of longer run time. Reduce to 300 for the old behaviour.
 - `hstats()`: By default, three-way interactions are not calculated anymore. Set `threeway_m` to 5 for the old behaviour.
 - Revised plots: The colors and color palettes have changed and can (also) be controlled via global options. For instance, to change the fill color of all bars, set `options(hstats.fill = new value)`. Value labels are more clear, and there are more options. Varying color/fill scales now use viridis (inferno). This can be modified on the fly or via `options(hstats.viridis_args = list(...))`.
-- "hstats_matrix" object: All statistics functions, e.g., `h2_pairwise()` or `perm_importance()`, now return a "hstats_matrix". The values are stored in `$M` and can be plotted via `plot()`. Other methods are: `dimnames()`, `rownames()`, `colnames()`, `dim()`, `nrow()`, `ncol()`, `head()`, `tail()`, and subsetting like a normal matrix. This allows, e.g, to select and plot only one column of the results.
+- "hstats_matrix" object: All statistics functions, e.g., `h2_pairwise()` or `perm_importance()`, now return a "hstats_matrix". The values are stored in `$M` and can be plotted via `plot()`. Other methods include: `dimnames()`, `rownames()`, `colnames()`, `dim()`, `nrow()`, `ncol()`, `head()`, `tail()`, and subsetting like a normal matrix. This allows, e.g, to select and plot only one column of the results.
 - `perm_importance()`: The `perms` argument has been changed to `m_rep`.
 - `print()` and `summary()` methods have been revised.
 - The arguments `w` (case weights) and `y` (response) can now also be passed as column *names*.
diff --git a/R/H2_pairwise.R b/R/H2_pairwise.R
index 5353747a..74e4ac9f 100644
--- a/R/H2_pairwise.R
+++ b/R/H2_pairwise.R
@@ -69,7 +69,6 @@
 #' s <- hstats(fit, X = iris[3:5], verbose = FALSE)
 #' x <- h2_pairwise(s)
 #' plot(x)
-#' plot(x[, "Sepal.Length"])
 h2_pairwise <- function(object, ...) {
   UseMethod("h2_pairwise")
 }
diff --git a/R/perm_importance.R b/R/perm_importance.R
index 55a97c34..dd262e94 100644
--- a/R/perm_importance.R
+++ b/R/perm_importance.R
@@ -49,7 +49,6 @@
 #' s
 #' plot(s)
 #' plot(s, swap_dim = TRUE, top_m = 2)
-#' plot(s[, "Sepal.Length"])
 perm_importance <- function(object, ...) {
   UseMethod("perm_importance")
 }
diff --git a/R/utils_statistics.R b/R/utils_statistics.R
index de2e2d05..1779f4a4 100644
--- a/R/utils_statistics.R
+++ b/R/utils_statistics.R
@@ -296,6 +296,31 @@ dimnames.hstats_matrix <- function(x) {
   dimnames(x[["M"]])
 }
 
+#' Dimnames (Replacement Method) of "hstats_matrix"
+#'
+#' This implies `colnames(x) <- ...`.
+#'
+#' @param x An object of class "hstats_matrix".
+#' @param value A list with rownames and column names compliant with `$M` (and `$SE`).
+#' @returns Like `x`, but with replaced dimnames.
+#' @examples
+#' fit <- lm(as.matrix(iris[1:2]) ~ Petal.Length + Petal.Width * Species, data = iris)
+#' s <- hstats(fit, X = iris[3:5], verbose = FALSE)
+#' x <- h2_overall(s)
+#' colnames(x) <- c("Sepal Length", "Sepal Width")
+#' plot(x)
+#' 
+#' rownames(x)[2:3] <- c("Petal Width", "Petal Length")
+#' plot(x)
+#' @export
+`dimnames<-.hstats_matrix` <- function(x, value) {
+  dimnames(x[["M"]]) <- value
+  if (!is.null(x[["SE"]])) {
+    dimnames(x[["SE"]]) <- value
+  }
+  x
+}
+
 #' Subsets "hstats_matrix" Object
 #'
 #' Use standard square bracket subsetting to select rows and/or columns of
diff --git a/man/H2_pairwise.Rd b/man/H2_pairwise.Rd
index b41db39f..1be604b7 100644
--- a/man/H2_pairwise.Rd
+++ b/man/H2_pairwise.Rd
@@ -120,7 +120,6 @@ fit <- lm(as.matrix(iris[1:2]) ~ Petal.Length + Petal.Width * Species, data = ir
 s <- hstats(fit, X = iris[3:5], verbose = FALSE)
 x <- h2_pairwise(s)
 plot(x)
-plot(x[, "Sepal.Length"])
 }
 \references{
 Friedman, Jerome H., and Bogdan E. Popescu. \emph{"Predictive Learning via Rule Ensembles."}
diff --git a/man/dimnames-set-.hstats_matrix.Rd b/man/dimnames-set-.hstats_matrix.Rd
new file mode 100644
index 00000000..7c225ae5
--- /dev/null
+++ b/man/dimnames-set-.hstats_matrix.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils_statistics.R
+\name{dimnames<-.hstats_matrix}
+\alias{dimnames<-.hstats_matrix}
+\title{Dimnames (Replacement Method) of "hstats_matrix"}
+\usage{
+\method{dimnames}{hstats_matrix}(x) <- value
+}
+\arguments{
+\item{x}{An object of class "hstats_matrix".}
+
+\item{value}{A list with rownames and column names compliant with \verb{$M} (and \verb{$SE}).}
+}
+\value{
+Like \code{x}, but with replaced dimnames.
+}
+\description{
+This implies \code{colnames(x) <- ...}.
+}
+\examples{
+fit <- lm(as.matrix(iris[1:2]) ~ Petal.Length + Petal.Width * Species, data = iris)
+s <- hstats(fit, X = iris[3:5], verbose = FALSE)
+x <- h2_overall(s)
+colnames(x) <- c("Sepal Length", "Sepal Width")
+plot(x)
+
+rownames(x)[2:3] <- c("Petal Width", "Petal Length")
+plot(x)
+}
diff --git a/man/perm_importance.Rd b/man/perm_importance.Rd
index 7cc422fe..9281bbd3 100644
--- a/man/perm_importance.Rd
+++ b/man/perm_importance.Rd
@@ -210,7 +210,6 @@ s <- perm_importance(fit, X = iris[3:5], y = iris[1:2], normalize = TRUE)
 s
 plot(s)
 plot(s, swap_dim = TRUE, top_m = 2)
-plot(s[, "Sepal.Length"])
 }
 \references{
 Fisher A., Rudin C., Dominici F. (2018). All Models are Wrong but many are Useful:
diff --git a/tests/testthat/test_statistics.R b/tests/testthat/test_statistics.R
index 4056ad16..a43d4232 100644
--- a/tests/testthat/test_statistics.R
+++ b/tests/testthat/test_statistics.R
@@ -47,17 +47,34 @@ test_that(".zap_small() works for matrix input", {
 fit <- lm(cbind(up = uptake, up2 = 2 * uptake) ~ Type * Treatment * conc, data = CO2)
 H <- hstats(fit, X = CO2[2:4], verbose = FALSE)
 s <- h2_pairwise(H)
+imp <- perm_importance(fit, CO2, v = c("Type", "Treatment", "conc"), y = "uptake")
 
 test_that("print() method does not give error", {
   capture_output(expect_no_error(print(s)))
+  capture_output(expect_no_error(print(s)))
 })
 
 test_that("dim() is correct", {
   expect_equal(dim(s), c(3L, 2L))
+  expect_equal(dim(imp), c(3L, 2L))
 })
 
 test_that("dimnames() is correct", {
   expect_equal(dimnames(s), list(rownames(s$M), colnames(s$M)))
+  expect_equal(dimnames(imp), list(rownames(imp$SE), colnames(imp$SE)))
+})
+
+test_that("dimnames() (replacement) works", {
+  s2 <- s
+  colnames(s2) <- c("y", "x")
+  rownames(s2) <- c("A", "B", "C")
+  expect_equal(colnames(s2), c("y", "x"))
+  expect_equal(rownames(s2), c("A", "B", "C"))
+  
+  imp2 <- imp
+  dimnames(imp2) <- list(c("A", "B", "C"), c("y", "x"))
+  expect_equal(colnames(imp2), c("y", "x"))
+  expect_equal(rownames(imp2), c("A", "B", "C"))
 })
 
 test_that("subsetting works", {