Such relationships are conventionally summarized with contingency tables, as in this (random) example:
Col 1 Col 2 Col 3 Col 4
Row 1 3 6 40 34
Row 2 18 6 9 1
Typically we are interested in comparing these data to values suggested by some default model, such as a null model of independent row and column proportions. When comparing the data to those values, the actual counts are important because they are proportional to the variances of the differences.
Consequently, a good visualization would clearly show the counts and their expected values, preferably organized to parallel the table.
Studies by psychologists and statisticians indicate that graphical elements like hue and shade do a relatively poor job in depicting quantities like counts. Although length and position tend to be clearest and most accurate, they are suited only for showing relative counts: that is, their proportions. Not good enough.
I therefore propose to represent any count $k$ by drawing $k$ distinct, non-overlapping identically-sized graphical symbols, so that each symbol clearly represents one thing that is counts. To make this work well, my experiments have found the following:
Clustering the symbols into a compact object seems to work better than positioning them randomly within a drawing area.
Overplotting the symbols on a polygon whose area represents the expectation permits a direct visual comparison of the count to its expectation. Rectangles, concentric with the symbol clusters, suffice for this purpose.
As a bonus, the standard error of each count, which is proportional to its square root, is thereby represented by the perimeter of its reference polygon. Although this is subtle, it's nice to see such a useful quantity appear naturally in the graphic.
People gravitate towards colorful graphics, but because colors might not reproduced (think of page charges in a research journal, for instance), I apply color to distinguish the cells but not to represent anything essential.
Here is an example of this solution for the table above:

It is immediately clear which cells have overly large counts and which have overly small ones. We even get a quick impression of how much they exceed or fall short of their expectations. With a little practice, you can learn to eyeball the chi-squared statistic from such a plot.
I have decorated the figure with the usual accompaniments: row and column labels to the left and top; row and column totals to the right and bottom; and the p-value of a test (in this case, Fisher's Exact test of independence as computed with a million simulated datasets).
For comparison, here is the visualization with randomly dispersed symbols:

Because the symbols are no longer clustered, it's useless to draw the reference rectangles. Instead, I have used the cell shading to represent expected values. (Darker is higher.) Although this method still works, I get more out of the first (clustered) version.
When either or both of the variables are ordered, the same visualization is effective provided the rows and columns follow the ordering.
Finally, this works well for $2\times 2$ tables. Here is one that came up in an analysis of an age discrimination case where it was alleged that older workers were preferentially fired. Indeed, the table looks a little incriminating because no younger people were let go at all:
Old Young
Kept 135 26
Fired 14 0
The visualization, however, indicates a close agreement between the observations and the expected values under the null hypothesis of no relationship with age:

The Fisher Exact test p-value of $0.134$ supports the visual impression.
Because I know people will ask for it, here is the R
code used to produce the figures.
m <- 2
n <- 4
set.seed(17)
shape <- .8
mu <- 180 / (m*n)
x <- matrix(rpois(m*n, rgamma(m*n, shape, shape/mu)), m, n)
if (is.null(colnames(x))) colnames(x) <- paste("Col", 1:n)
if (is.null(rownames(x))) rownames(x) <- paste("Row", 1:m)
breaks.x <- seq(0, n, length.out=n+1)
breaks.y <- rev(seq(0, m, length.out=m+1))
#
# Testing.
#
p.value <- signif(fisher.test(x, simulate.p.value=TRUE, B=1e6)$p.value, 3)
print(x)
#
# Set up plotting parameters.
#
random <- TRUE
h <- sample.int(m*n)
colors <- matrix(hsv(h / length(h), 0.9, 0.8, 1/2), nrow(x), ncol(x))
eps <- (1 - 1/(1.08))/2 # (Makes the plotting area exactly the right size.)
lim <- c(eps, 1-eps)
plot(lim*n, lim*m, type="n", xaxt="n", yaxt="n", bty="n", xlab="", ylab="",
xaxs="r", yaxs="r", asp=m/n,
main=substitute(paste("A ", m %*% n, " Table"), list(m=m, n=n)))
mtext(bquote(italic(p)==.(p.value)), side=1, line=2)
#
# Expectations.
#
gamma <- 6/3 # (Values above 1 reduce the background contrast.)
p.row <- rowSums(x)/sum(x)
p.col <- colSums(x)/sum(x)
if (isTRUE(random)) {
for (i in 1:m) {
polygon(c(range(breaks.x), rev(range(breaks.x))), rep(breaks.y[0:1+i], each=2),
col=hsv(0,0,0, p.row[i]^gamma))
}
for (j in 1:n) {
polygon(breaks.x[c(j,j+1,j+1,j)], rep(range(breaks.y), each=2),
col=hsv(0,0,0, p.col[j]^gamma))
}
} else {
for (i in 1:m) {
for (j in 1:n) {
p <- p.row[i] * p.col[j]
h <- (1 - (breaks.y[i] - breaks.y[i+1]) * sqrt(p))/2
w <- (1 - (breaks.x[j+1] - breaks.x[j]) * sqrt(p))/2
polygon(c(breaks.x[j]+w, breaks.x[j+1]-w, breaks.x[j+1]-w, breaks.x[j]+w),
c(breaks.y[i+1]+w, breaks.y[i+1]+w, breaks.y[i]-w, breaks.y[i]-w),
col=hsv(0,0,1/2))
}
}
}
#
# Borders.
#
gray <- hsv(0,0,5/6)
invisible(sapply(breaks.x, function(x) lines(rep(x,2), range(breaks.y), col=gray)))
invisible(sapply(breaks.y, function(y) lines(range(breaks.x), rep(y,2), col=gray)))
polygon(c(range(breaks.x), rev(range(breaks.x))), rep(range(breaks.y), each=2))
#
# Labels.
#
at <- (breaks.y[-1] + breaks.y[-(m+1)])/2
mtext(rownames(x), at=at, side=2, line=1/4)
mtext(rowSums(x), at=at, side=4, line=1/4)
at <- (breaks.x[-1] + breaks.x[-(n+1)])/2
mtext(colnames(x), at=at, side=3, line=0)
mtext(colSums(x), at=at, side=1, line=1/4)
#
# Samples.
#
runif2 <- function(n, ncol, nrow, lower.x=0, upper.x=1, lower.y=0, upper.y=1, random=TRUE) {
if (n > nrow*ncol) {
warning("Unable to generate enough samples")
n <- nrow*ncol
}
if (isTRUE(random)) {
i <- sample.int(nrow*ncol, n) - 1
} else {
# i <- seq_len(n) - 1
k <- order(outer(nrow*(1:ncol-(ncol+1)/2), ncol*(1:nrow-(nrow+1)/2), function(x,y) x^2+y^2))
i <- k[seq_len(n)] - 1
}
j <- (i %% ncol + 1/2) / ncol * (upper.y - lower.y) + lower.y
i <- (i %/% ncol + 1/2) / nrow * (upper.x - lower.x) + lower.x
cbind(i,j)
}
### Adjust the `400` to make the symbols barely overlap ###
cex <- 1 / sqrt(max(x)/400*max(m,n))
eps.x <- eps.y <- 0.05
u <- sqrt(max(x)/ (m*n))
u <- ceiling(u)
for (i in 1:m) {
for (j in 1:n) {
points(runif2(x[i,j], ceiling(m*u), ceiling(n*u),
breaks.x[j]+eps.x, breaks.x[j+1]-eps.x,
breaks.y[i+1]+eps.y, breaks.y[i]-eps.y,
random=random),
pch=22, cex=cex, col=colors[i,j], bg=colors[i,j])
}
}