diff --git a/R/class-plotdata-mosaic.R b/R/class-plotdata-mosaic.R
index 862a83f..e207c41 100644
--- a/R/class-plotdata-mosaic.R
+++ b/R/class-plotdata-mosaic.R
@@ -32,7 +32,6 @@ newMosaicPD <- function(.dt = data.table::data.table(),
   
   if (!isEvil) {
     if (statistic == 'all') {
-      # currently only valid for 2x2
       attr$statsTable <- panelAllStats(.pd, x, y, panel, columnReferenceValue, rowReferenceValue)
       veupathUtils::logWithTime('Calculated all relevant statistics.', verbose)
     } else if (statistic == 'chiSq') {
@@ -153,7 +152,7 @@ mosaic.dt <- function(data, variables,
     }
     #na.rm should be safe, since x and y axes will later have NA removed anyhow in the plot.data parent class
     if ((data.table::uniqueN(data[[x]], na.rm = TRUE) > 2 || data.table::uniqueN(data[[y]], na.rm = TRUE) > 2) && statistic == 'all') {
-      stop('Odds ratio and relative risk can only be calculated for 2x2 contingency tables. Please use statistic `chiSq` instead.')
+      warning('Odds ratio and relative risk can only be calculated for 2x2 contingency tables. Only the `chiSq` statistic will be returned.')
     }
   } else {
     if (data.table::uniqueN(data[[x]], na.rm = TRUE) > 2 || data.table::uniqueN(data[[y]], na.rm = TRUE) > 2) {
diff --git a/R/methods-ContingencyTable.R b/R/methods-ContingencyTable.R
index 973796b..d4b1e6f 100644
--- a/R/methods-ContingencyTable.R
+++ b/R/methods-ContingencyTable.R
@@ -59,11 +59,15 @@ setGeneric("orderByReferenceValues",
 )
 
 #' @export
-setMethod("orderByReferenceValues", signature("TwoByTwoTable"), function(object) {
+setMethod("orderByReferenceValues", signature("ContingencyTable"), function(object) {
   tbl <- object@data
   columnReferenceValue <- object@columnReferenceValue
   rowReferenceValue <- object@rowReferenceValue
 
+  # If there are no ref values, just return. This prevents the need for more
+  # complicated handling of contingency tables that are not 2x2s.
+  if (is.na(columnReferenceValue) && is.na(rowReferenceValue)) return(object)
+
   if (!is.na(columnReferenceValue)) {
     if (attributes(tbl)$dimnames[[2]][1] != columnReferenceValue) {
       attr <- attributes(tbl)
@@ -104,8 +108,10 @@ setGeneric("chiSqResults",
 )
 
 #' @export
-setMethod("chiSqResults", signature("TwoByTwoTable"), function(object) {
+setMethod("chiSqResults", signature("ContingencyTable"), function(object) {
+
   object <- orderByReferenceValues(object)
+
   tbl <- object@data
 
   chisq <- try(suppressWarnings(stats::chisq.test(tbl)))
@@ -460,6 +466,16 @@ setGeneric("allStats",
   signature = "object"
 )
 
+#' @importFrom S4Vectors SimpleList
+#' @export
+setMethod("allStats", signature("ContingencyTable"), function(object) {
+   return(veupathUtils::StatisticList(S4Vectors::SimpleList(
+    chiSqResults(object)
+    # fishersTest(object)
+   ))) 
+})
+
+
 #' @importFrom S4Vectors SimpleList
 #' @export
 setMethod("allStats", signature("TwoByTwoTable"), function(object) {
diff --git a/R/panel.R b/R/panel.R
index 1ae97ae..9254ce3 100644
--- a/R/panel.R
+++ b/R/panel.R
@@ -7,31 +7,59 @@ panelAllStats <- function(data, x, y, panel = NULL, columnReferenceValue = NA_ch
   x <- 'x'
   y <- 'y'
 
-  if (is.null(panel)) {
-    tbl <- tableXY(data)
-    tbl <- TwoByTwoTable('data'=tbl, 'columnReferenceValue'=columnReferenceValue, 'rowReferenceValue'=rowReferenceValue)
-    statistics <- allStats(tbl)
-    dt <- veupathUtils::as.data.table(statistics)
-  } else {
-    buildTwoByTwo <- function(tbl) {
-      TwoByTwoTable('data' = tbl, 'columnReferenceValue' = columnReferenceValue, 'rowReferenceValue' = rowReferenceValue)
+  # If both x and y only have two unique values, then we have a 2x2 table, and we should compute
+  # a lot of stats. Otherwise, it's an RxC
+  if (uniqueN(data[[x]]) <= 2 && uniqueN(data[[y]]) <= 2) {
+    # 2x2 Stats
+    if (is.null(panel)) {
+      tbl <- tableXY(data)
+      tbl <- TwoByTwoTable('data'=tbl, 'columnReferenceValue'=columnReferenceValue, 'rowReferenceValue'=rowReferenceValue)
+      statistics <- allStats(tbl)
+      dt <- veupathUtils::as.data.table(statistics)
+    } else {
+      buildTwoByTwo <- function(tbl) {
+        TwoByTwoTable('data' = tbl, 'columnReferenceValue' = columnReferenceValue, 'rowReferenceValue' = rowReferenceValue)
+      }
+
+      dt.list <- split(data, list(data[[panel]]))
+      dt.list <- lapply(dt.list, tableXY)
+      dt.list <- lapply(dt.list, buildTwoByTwo)
+      dt.list <- lapply(dt.list, allStats)
+      dt.list <- lapply(dt.list, veupathUtils::as.data.table)
+      colNames <- names(dt.list[[1]])
+      dt <- data.table::as.data.table(lapply(as.list(colNames), function(name) { lapply( dt.list, function(x) {x[[name]]} ) } ))
+      data.table::setnames(dt, colNames)
+      #dt <- purrr::reduce(dt.list, rbind)
+      dt[[panel]] <- names(dt.list)
     }
+  } else {
+    # RxC Stats. For now just chi squared.
+    if (is.null(panel)) {
+      tbl <- tableXY(data)
+      tbl <- ContingencyTable('data'=tbl, 'columnReferenceValue'=columnReferenceValue, 'rowReferenceValue'=rowReferenceValue)
+      statistics <- allStats(tbl)
+      dt <- veupathUtils::as.data.table(statistics)
+    } else {
+      buildContingency <- function(tbl) {
+        ContingencyTable('data' = tbl, 'columnReferenceValue' = columnReferenceValue, 'rowReferenceValue' = rowReferenceValue)
+      }
 
-    dt.list <- split(data, list(data[[panel]]))
-    dt.list <- lapply(dt.list, tableXY)
-    dt.list <- lapply(dt.list, buildTwoByTwo)
-    dt.list <- lapply(dt.list, allStats)
-    dt.list <- lapply(dt.list, veupathUtils::as.data.table)
-    colNames <- names(dt.list[[1]])
-    dt <- data.table::as.data.table(lapply(as.list(colNames), function(name) { lapply( dt.list, function(x) {x[[name]]} ) } ))
-    data.table::setnames(dt, colNames)
-    #dt <- purrr::reduce(dt.list, rbind)
-    dt[[panel]] <- names(dt.list)
+      dt.list <- split(data, list(data[[panel]]))
+      dt.list <- lapply(dt.list, tableXY)
+      dt.list <- lapply(dt.list, buildContingency)
+      dt.list <- lapply(dt.list, allStats)
+      dt.list <- lapply(dt.list, veupathUtils::as.data.table)
+      colNames <- names(dt.list[[1]])
+      dt <- data.table::as.data.table(lapply(as.list(colNames), function(name) { lapply( dt.list, function(x) {x[[name]]} ) } ))
+      data.table::setnames(dt, colNames)
+      dt[[panel]] <- names(dt.list)
+    }
   }
 
   return(dt) 
 }
 
+# No longer in use?
 panelBothRatios <- function(data, x, y, panel = NULL) {
   names(data)[names(data) == x] <- 'x'
   names(data)[names(data) == y] <- 'y'
diff --git a/tests/testthat/test-mosaic.R b/tests/testthat/test-mosaic.R
index b6ada5a..f778eea 100644
--- a/tests/testthat/test-mosaic.R
+++ b/tests/testthat/test-mosaic.R
@@ -184,6 +184,21 @@ test_that("mosaic.dt() returns a valid plot.data mosaic object", {
   expect_equal(names(namedAttrList$statsTable), c('chisq','pvalue', 'degreesFreedom','entity.cat4'))
 
 
+  # Testing an RxC with 'all' statistics
+  dt <- mosaic.dt(df, variables, statistic = 'all')
+  expect_is(dt, 'plot.data')
+  expect_is(dt, 'mosaic')
+  namedAttrList <- getPDAttributes(dt)
+  expect_equal(names(namedAttrList),c('variables', 'completeCasesAllVars','completeCasesAxesVars','completeCasesTable','sampleSizeTable','statsTable'))
+  completeCases <- completeCasesTable(dt)
+  expect_equal(names(completeCases), c('variableDetails','completeCases'))
+  expect_equal(nrow(completeCases), 3)
+  sampleSizes <- sampleSizeTable(dt)
+  expect_equal(names(sampleSizes), c('entity.cat4','entity.cat7','size'))
+  expect_equal(nrow(sampleSizes), 4)
+  expect_equal(names(namedAttrList$statsTable), c('chiSq','entity.cat4'))
+
+
   # Ensure sampleSizeTable and completeCasesTable do not get returned if we do not ask for them.
   dt <- mosaic.dt(df, variables, sampleSizes = FALSE, completeCases = FALSE)
   expect_is(dt, 'plot.data')
@@ -312,6 +327,7 @@ test_that("mosaic.dt() returns plot data and config of the appropriate types", {
   expect_equal(class(unlist(namedAttrList$statsTable$degreesFreedom)), c('scalar', 'numeric'))
   expect_equal(class(unlist(namedAttrList$statsTable$pvalue)), c('scalar', 'numeric'))
   expect_equal(class(unlist(namedAttrList$statsTable$entity.cat4)), 'character')
+
 })
 
 test_that("mosaic.dt() returns an appropriately sized data.table", {
@@ -522,6 +538,24 @@ test_that("mosaic.dt() returns an appropriately sized data.table", {
   expect_equal(names(sampleSizeTable),c('entity.cat7','size'))
   expect_equal(class(sampleSizeTable$entity.cat7[[1]]), 'character')
 
+  # With 'all' statistics for an RxC
+  dt <- mosaic.dt(df, variables, statistic = 'all')
+  expect_is(dt, 'data.table')
+  expect_is(dt$value, 'list')
+  expect_is(dt$value[[1]], 'list')
+  expect_equal(nrow(dt),1)
+  expect_equal(names(dt),c('xLabel', 'yLabel', 'value'))
+  expect_equal(dt$xLabel[[1]],paste0("cat7_", letters[1:7]))
+  expect_equal(dt$yLabel[[1]][[1]],paste0("cat3_", letters[1:3]))
+  expect_equal(length(dt$value[[1]]),7)
+  expect_equal(length(dt$value[[1]][[1]]),3)
+  statsTable <- statsTable(dt)
+  expect_equal(names(statsTable), c(c('chiSq')))
+  sampleSizeTable <- sampleSizeTable(dt)
+  expect_equal(names(sampleSizeTable),c('entity.cat7','size'))
+  expect_equal(class(sampleSizeTable$entity.cat7[[1]]), 'character')
+
+
   # With factors
   variables <- new("VariableMetadataList", SimpleList(
     new("VariableMetadata",
@@ -878,6 +912,52 @@ test_that("mosaic() returns appropriately formatted json", {
   expect_equal(class(jsonList$sampleSizeTable$xVariableDetails$value[[1]]), 'character')
   expect_equal(names(jsonList$completeCasesTable), c('variableDetails', 'completeCases'))
   expect_equal(names(jsonList$completeCasesTable$variableDetails), c('variableId','entityId')) 
+
+
+  # Check json structure of stats table when statistics = 'all' for rxc
+  variables <- new("VariableMetadataList", SimpleList(
+    new("VariableMetadata",
+      variableClass = new("VariableClass", value = 'native'),
+      variableSpec = new("VariableSpec", variableId = 'int6', entityId = 'entity'),
+      plotReference = new("PlotReference", value = 'yAxis'),
+      dataType = new("DataType", value = 'NUMBER'),
+      dataShape = new("DataShape", value = 'CONTINUOUS')),
+    new("VariableMetadata",
+      variableClass = new("VariableClass", value = 'native'),
+      variableSpec = new("VariableSpec", variableId = 'int7', entityId = 'entity'),
+      plotReference = new("PlotReference", value = 'xAxis'),
+      dataType = new("DataType", value = 'NUMBER'),
+      dataShape = new("DataShape", value = 'CONTINUOUS')),
+    new("VariableMetadata",
+      variableClass = new("VariableClass", value = 'native'),
+      variableSpec = new("VariableSpec", variableId = 'cat3', entityId = 'entity'),
+      plotReference = new("PlotReference", value = 'facet1'),
+      dataType = new("DataType", value = 'STRING'),
+      dataShape = new("DataShape", value = 'CATEGORICAL'))
+  ))
+
+  dt <- mosaic.dt(df, variables, statistic = 'all')
+  outJson <- getJSON(dt, FALSE)
+  jsonList <- jsonlite::fromJSON(outJson)
+
+  expect_equal(names(jsonList),c('mosaic','sampleSizeTable','statsTable','completeCasesTable'))
+  expect_equal(names(jsonList$mosaic),c('data','config'))
+  expect_equal(names(jsonList$mosaic$data),c('xLabel','yLabel','value','facetVariableDetails'))
+  expect_equal(names(jsonList$mosaic$data$facetVariableDetails[[1]]),c('variableId','entityId','value'))
+  expect_equal(length(jsonList$mosaic$data$facetVariableDetails), 3)
+  expect_equal(jsonList$mosaic$data$facetVariableDetails[[1]]$variableId, 'cat3')
+  expect_equal(names(jsonList$mosaic$config),c('variables','completeCasesAllVars','completeCasesAxesVars'))
+  expect_equal(names(jsonList$mosaic$config$variables$variableSpec),c('variableId','entityId'))
+  expect_equal(jsonList$mosaic$config$variables$variableSpec$variableId, c('int6','int7','cat3'))
+  expect_equal(names(jsonList$sampleSizeTable),c('facetVariableDetails','xVariableDetails','size'))
+  expect_equal(class(jsonList$sampleSizeTable$facetVariableDetails[[1]]$value), 'character')
+  expect_equal(class(jsonList$sampleSizeTable$xVariableDetails$value[[1]]), 'character')
+  expect_equal(jsonList$sampleSizeTable$xVariableDetails$variableId[[1]], 'int7')
+  expect_equal(names(jsonList$statsTable),c('chiSq', 'facetVariableDetails'))
+  expect_equal(names(jsonList$completeCasesTable), c('variableDetails', 'completeCases'))
+  expect_equal(names(jsonList$completeCasesTable$variableDetails), c('variableId','entityId'))
+  expect_equal(jsonList$completeCasesTable$variableDetails$variableId, c('int7', 'int6', 'cat3'))
+
 })
 
 test_that("mosaic.dt() returns correct information about missing data", {