How do we visualize the Decomposition of Variance formula $\text{Var}[y] = \text{Var}_x[\text{E}[y|x]] + \text{E}_x[\text{Var}[y|x]]$? The Decomposition of Variance formula is
$\text{Var}[Y] = \text{Var}_X[\text{E}[Y|X]] + \text{E}_X[\text{Var}[Y|X]]$.
This can be described as follows: the variance of y decomposes into the variance of the conditional mean function plus the expected variance around the conditional mean.
Suppose we have the simple linear regression model
$$
Y = \beta_0 + \beta_1 X + \varepsilon,
$$
with normally distributed error. How do we visualize the Decomposition of Variance in this context, that is, where we have a line on a 2D plot with the axes labelled $x$ and $y$.
I can see that the left hand side $\text{Var}[y]$ will be associated with how spread out the distribution of $Y$ is with respect to the $y$ axis but what about the terms on the right hand side, how do we interpret them?
 A: I love graphical displays.  Here are two that illustrate the right hand side of the law of total variance nicely.  First, some code for a linear but heteroscedastic regression.
set.seed(12345)
nsim = 100
X = runif(nsim, 40,120)
Y = 1 + 0.3*X + rnorm(nsim, 0, 0.15*X)

Cond.Mean = 1 + 0.3*X      # Conditional Mean
Cond.SD = 0.15*X           # Conditional Standard Deviation

plot(X,Y, main = "Illustrating Variance of Conditional Mean")
abline(1,.3)
rug(Cond.Mean, side=2) 

The resulting graph is as follows:

The vertical spread of the of the data ticks (the "rug") on the vertical axis represents the variance of the conditional mean values, or $Var_X[E[Y|X]]$.  Notice that this range is a lot smaller than the overall vertical data range, which represents $Var[Y]$.
To visualize the mean of the conditional variance, add the $\pm \sigma_{Y|X}$ bands to the scatter as follows:
plot(X,Y, main = "Illustrating Mean of Conditional Variance")
abline(1,.3)
abline(1,.15, lty=2)
abline(1,.45, lty=2)
rug(X)

The resulting graph is as follows:

Now, for every $x$ value on the "floor" (the "rug"), there is a different vertical spread of potential $Y$ values, as indicated by the $\pm \sigma_{Y|X}$ bands. Each of these spreads represents (via squaring) a conditional variance $Var[Y|X=x]$. The average of all these conditional variances is equal to the other term on the right side, $E_X[Var[Y|X]]$.
You can attempt to verify the equality using
var(Y)
var(Cond.Mean) + mean(Cond.SD^2)

but there is a lot of finite-sample variability, so the results are not that close for this small simulation.  On the other hand, if you keep the same seed and change nsim to 20000000, the results are very close, 204.05 and 204.01.
