Dichtomising continious treatment variables In the past, I have repeatedly seen studies that dichtomize continuous treatment variables into a binary treatment dummy. I feel that this cannot be good practice. 


*

*we loose important information

*we might underestimate the variance
But can this also result in a bias of the simple OLS estimate? 
 A: Well, as we discuss in the comments, a model that includes a dichotomized predictor is inherently misspecified if the true effect is not discrete. (Or even if the true effect is discrete, but the cutpoint is different from what we believe.)
Thus, the true non-zero parameter is not even in the model, i.e., the model forces it to be zero. It's biased. Conversely, we have a predictor in the model that should not be there, since it is not in the DGP, namely the dichotomized predictor. Since the un-dichotomized predictor does have an effect, the parameter estimate of the dichotomized one will almost certainly be biased.
Here is a little example. I generate $x\sim N(0,1)$ and 
$$ y = 5+x+\epsilon, \quad \epsilon\sim N(0,1). $$
Thus, the true parameters are $\beta_0=5$ and $\beta_x=1$.
I then discretize $x$ based on its sample mean to obtain $x_{\text{discrete}}$. Finally, I compare two models:
$$ y \sim x_{\text{discrete}} \quad\text{and}\quad y\sim x+x_{\text{discrete}}.$$
The second model encompasses the true DGP by adding an extraneous $x_{\text{discrete}}$.
I run this 10,000 times, using a sample size of $n=100$ in each case. Below are beanplots of the parameter estimates, with the true values indicated by red diamonds. We see clearly that both coefficients of the misspecified model are biased, while the ones of the "too large" model aren't. (Of course, we should really use the true model and get smaller parameter variances.)

R code:
n_sims <- 1e4
nn <- 100

coefs_sims_true <- matrix(NA,nrow=n_sims,ncol=3)
coefs_sims_discrete <- matrix(NA,nrow=n_sims,ncol=2)

pb <- winProgressBar(max=n_sims)
for ( ii in 1:n_sims ) {
    setWinProgressBar(pb,ii,paste(ii,"of",n_sims))
    set.seed(ii)

    dataset <- data.frame(x=rnorm(nn))
    dataset$y <- 5+dataset$x+rnorm(nn)
    dataset$x_discrete <- dataset$x<mean(dataset$x)

    model_true <- lm(y~x+x_discrete,dataset)
    coefs_sims_true[ii,] <- coefficients(model_true)
    model_discrete <- lm(y~x_discrete,dataset)
    coefs_sims_discrete[ii,] <- coefficients(model_discrete)
}
close(pb)

summary(coefs_sims_discrete)
summary(coefs_sims_true)

library(beanplot)

opar <- par(mfrow=c(1,2))
    beanplot(data.frame(coefs_sims_discrete),what=c(0,1,0,0),
      col="lightgrey",xaxt="n",main="Discrete model")
    axis(1,at=1:ncol(coefs_sims_discrete),labels=names(coefficients(model_discrete)))
    points(c(1,2),c(5,0),pch=23,col="red",bg="red",cex=1.5)

    beanplot(data.frame(coefs_sims_true),what=c(0,1,0,0),
      col="lightgrey",xaxt="n",main="Full model")
    axis(1,at=1:ncol(coefs_sims_true),labels=names(coefficients(model_true)))
    points(1:3,c(5,1,0),pch=23,col="red",bg="red",cex=1.5)
par(opar)

