diff --git a/R/minc_anatomy.R b/R/minc_anatomy.R
index b71000f1..1ca3d3d8 100644
--- a/R/minc_anatomy.R
+++ b/R/minc_anatomy.R
@@ -878,7 +878,6 @@ anatLm <- function(formula, data, anat, subset=NULL, weights = NULL) {
       matrixName = formula[[2]]
       matrixFound = TRUE
       data.matrix.left <- t(anat[mf[,"(rowcount)"],])
-      #data.matrix.left <- vertexTable(filenames)
       data.matrix.right <- t(term)
     }  
   }
@@ -909,7 +908,11 @@ anatLm <- function(formula, data, anat, subset=NULL, weights = NULL) {
     rows = append(rows,matrixName)
   }	
   
-  
+  if(matrixFound){
+    stop("A term in your model is a matrix, this use of `anatLm` is deprecated, "
+       , "Please merge the two anatomy matrices using rowbind or similar and use a group "
+       , "add a group indicator to your `data.frame`")
+  }
   
   # Call subroutine based on whether matrix was found
   if(!matrixFound) {    	
diff --git a/tests/testthat/test_anatLm.R b/tests/testthat/test_anatLm.R
index 0291eb57..fd9dc09f 100644
--- a/tests/testthat/test_anatLm.R
+++ b/tests/testthat/test_anatLm.R
@@ -119,3 +119,12 @@ test_that("Weighted anatLm works", {
   expect_equivalent(sapply(lmods, AIC), AIC(alm))
 })
 
+context("Matrix mode")
+
+test_that("Test that passing a matrix on the RHS fails", {
+  d <- data.frame(y = rnorm(10))
+  d$x <- scale(rnorm(10))
+
+  expect_error(anatLm(~ x, data = d, anat = as.matrix(d$y)), "matrix")
+})
+