Let's switch from $r_t$ to the more standard notation $Y_t$, which is the potential outcome corresponding to setting treatment $z$ to $t$. The whole point of propensity score analysis is to estimate a causal quantity, usually $E[Y_1-Y_0]$, using observed data. Not all of the potential outcomes are observed because only the potential outcome corresponding to the treatment actually received is observed for each individual.
The causal estimand $E[Y_1-Y_0]$ is known as the ATE. The estimand described in that passage is $E[Y_1-Y_0|b(x)=k]$, which is the conditional average treatment effect for those with balancing score $b(x)=k$. The relevant balancing score is the propensity score, so, for example, the effect of interest might be the causal effect of treatment for those with a propensity score ($b(x)$) of .8 (i.e., $k$).
The estimator is the difference between the conditional means of $Y$ at each level of $z$ adjusted for the balancing score $b(x)$ at $b(x)=k$; more precisely, it's $E[Y|z=1,b(x)=k]-E[Y|z=0,b(x)=k]$. This quantity only relies on observed values (i.e., $Y$, $z$, and $b(x)$), but it is unbiased for the causal quantity $E[Y_1-Y_0|b(x_i)=k]$, which otherwise involves counterfactual values.
Here's how you might get that estimator:
- Regress $Y$ on $b(x)$ for those with $z=1$. Get the coefficients $\hat{\alpha}_1$ and $\hat{\beta}_1$ from this regression.
- Regress $Y$ on $b(x)$ for those with $z=0$. Get the coefficients $\hat{\alpha}_0$ and $\hat{\beta}_0$ from this regression.
- Choose a level of the balancing score $b(X)$, let's say, $k$.
- Plug in those values into the formula $(\hat{\alpha}_1 - \hat{\alpha}_0) + (\hat{\beta}_1 - \hat{\beta}_0)k$.
The estimated quantity is the conditional ATE at $b(x)=k$.
It's possible to get this in one single step using a regression model, like the one you proposed. The problem with that regression model is it assumes $\beta_1=\beta_0$. If that's true, the coefficient $\theta$ in your model is the average treatment effect. If not, $\theta$ is biased.
The correct model to use would be $E[Y|z,b(x)]=\lambda_0 + \tau z + \lambda_1(b(x)-k) + \lambda_2(b(x)-k)z$. Now $\tau$ represents the conditional ATE at $b(x)=k$.
Following up on some matters in the comments...
Regression adjustment is one way to use the propensity score to estimate a causal effect. It requires that the relationship between the propensity score and the outcome is linear (or correctly modeled) at each level of treatment. If there is no treatment effect modification by the propensity score, a regression of the outcome on the treatment and the propensity score, as you have done, will indeed estimate a valid causal effect. The coefficient on treatment ($\theta$ in your model) corresponds to the ATE.
If there is treatment effect modification, which means the effect of treatment depends on one's propensity score, you need to model the interaction between the treatment and the propensity score on the outcome. If you center the propensity score at some value $k$, the effect estimate is the conditional ATE at the centered value (i.e., $k$) of the propensity score, which is equal to the quantity estimated using $(\hat{\alpha}_1 - \hat{\alpha}_0) + (\hat{\beta}_1 - \hat{\beta}_0)k$. With a linear outcome model, the ATE is the conditional ATE evaluated at the mean of the propensity score in the sample. Within a nonlinear outcome model, estimating the ATE requires using a marginal effects procedure like the margins
command in Stata.
You should always allow the treatment effect to vary based on the propensity score or any covariates in your model, which means including an interaction between treatment and that covariate in the outcome model. Centering the covariate at its mean in a linear model allows the coefficient on the main effect of treatment to be equal the ATE. The propensity score literature may not be super clear about this, but it is the best thing to do, and, as you have found, R&R83 even recommend it.