A regression decomposition for dummy variables Suppose that I have a regression model
\begin{align}
Y = X\beta_X + Z_1 \beta_{Z_1} + Z_2 \beta_{Z_2} + \varepsilon
\end{align}
where $X$ are some controls, and $Z_1$ and $Z_2$ are two sets of dummy ($\{0,1\}$-valued) variables. I assume that $Z_2$ does not depend upon $Z_1$, but $Z_1$ might depend upon $Z_2$ in the sense that the probability that an element of $Z_1$ might be affected by the value of $Z_2$. I would like to decompose the information (variation) in $Z_2$ into a component that affects $Y$ only through $Z_1$, and a residual that affects $Y$ directly. 
If $Z_1$ were continuous, I could do something like regress (each coordinate of) $Z_1$ onto $Z_2$, then regress $Y$ on the fitted values. The natural analogue of this, say fitting a logistic (or similar) regression model for $Z_1$ gives up linearity, losing interpretability, and takes much longer to compute.
I would appreciate any insight on literature addressing this or similar problems.
 A: In the approach you write for the case where the r.v.'s were continuous, you would regress $Z_1$ on $Z_2$ and use the fitted values $\hat Z_1$ in the $Y$-regression. But $\hat Z_1 = \hat E(Z_1\mid Z_2)$. So what you need is the conditional expectation function of $Z_1$ given $Z_2$.
This paper has a nice exposition for the bivariate Bernoulli distribution (as an introduction to the multivariate case).
Denoting  the joint probability $P(Z_1 = 0, Z_2=0)=p(0,0) $ and analogously for the other three possible combinations, the conditional density of $Z_1\mid Z_2$ is (eq. 2.11 of paper)
$$P(Z_1=z_1\mid Z_2=z_2) = \left(\frac {p(1,z_2)}{p(1,z_2)+p(0,z_2)}\right)^{z_1}\left(\frac {p(0,z_2)}{p(1,z_2)+p(0,z_2)}\right)^{1-z_1} $$
Then the conditional expectation function we are looking for is 
$$E(Z_1\mid Z_2=z_2) = \frac {p(1,z_2)}{p(1,z_2)+p(0,z_2)} \equiv Z_1^* $$
You will then have to approximate the joint probabilities, creating a standard $2 \times 2$ contingency table of relative empirical frequencies -from which actually you can also test your dependence suspicion: theoretically, $Z_1$ and $Z_2$ will be independent if and only if $p(0,0)\cdot p(1,1) = p(0,1)\cdot p(1,0)$ and you can test this by a Fisher exact test.
Assuming that independence is rejected, then the contigency table will give you the estimated probabilities that you need and you will obtain the regressor
$$\hat Z_1^* = \frac {\hat p(1,z_2)}{\hat p(1,z_2)+\hat p(0,z_2)}$$
..which will produce a series of numbers, as for each observation $j$ of the sample we will use the ratio that corresponds to the value $z_{2j}$ of $Z_2$, i.e. if $z_{2j}=1$ we will have 
$$\hat z_{1j}^* = \frac {\hat p(1,1)}{\hat p(1,1)+\hat p(0,1)}$$
etc. $\hat Z_1^*$ is still a binary variable, but not taking the $\{0,1\}$ values any more.
