How can I partition a distribution into two sub-populations with fixed bias? (simulation) I am trying to simulate a selection model for a variable $Y$ dependent on covariate vector $X$, so that two groups/sub-sets $S=(0,1)$ of observations on $Y$ are created, which have a fixed difference in expectations (means) $\theta$ (i.e., selection bias) $$E(Y|S=1)-E(Y|S=0)= \theta$$ 
I am looking for solutions to approach this problem, where in principle I could also simulate two variables $Y_1$ and $Y_2$ and overlay these as a mixture distribution $Y$. However, note that the selection mechanism $S$ should be conditionally independent of $Y$ given $X$ (ignorable selection) $$P(S|Y,X) = P(S|X),$$ but not be marginally independent $$P(S|Y) \ne P(S).$$
Whereas it is easy to parameterize a model for $P(S|X)$, and also relate $E(Y|X)$, it seems very difficult to me to specify this model to achieve a fixed $\theta$ and keep variances of $Y$ across sub-sets approximately constant (the latter is not absolutely necessary, but desirable). If the mixture approach is chosen, I do not know how to overlay distributions in dependence of $X$. Thanks for suggestions!
 A: Maybe I am oversimplifying this but why don't you draw for the unselected group $Y_0 \sim N(\gamma,\sigma_0)$ and for the selected group $Y_1 \sim N(\gamma+\theta,\sigma_1)$ with
$$Y = \left\{ \begin{array}{ll} Y_0 & \mbox{if $S=0 $};\\ Y_1 & \mbox{if $S=1 $}.\end{array} \right.$$
You have $$E(Y|S=1)-E(Y|S=0)= \theta$$
Now either define $X$ as a function of selection $X=f(S)$ or as a scalar transformation $X=cY$. Either of these will give you ignorability $P(S|Y,X) = P(S|X)$ since knowing $X$ either makes knowing $Y$ uninformative about selection.
That said, I doubt either of these options is what you had in mind since it renders your selection model somewhat silly. Beyond these options I find it difficult to imagine an alternative schema in which $Y$ does not provide information with regards selection conditional upon $X$.
Perhaps there are some other transformations of $Y$ that will get you there. Perhaps something like $X$ being a continuous function defined as $Y$ is above and $Y$ being a function which contains some of the information from $X$ though not all of it such as $Y=Round(X)+c$. Thus knowing $X$ gives you all of the information you could have gained by knowing $Y$.
A: If I understood your question correctly, you want observations $(X_i, Y_i, treat_i)$ where $Y$ is different across treatment groups, but where, conditionally on $X$, $Y$ is independent of treatment. If that is the case, this R script works: 
n.obs <- 1e3
prop.1 <- .5
n1 <- ceiling(prop.1 * n.obs)
sub1 <- sample(1:n.obs, n1)
treatment <- (1:n.obs) %in% sub1

x <- rnorm(n.obs, mean= -1)
x[sub1] <- rnorm(n1, mean= 1)
y <- 1 + 2 * x + rnorm(n.obs)

plot(x,y, col= 1 + treatment)
by(y, treatment, mean)
summary(glm(treatment~x+y, family = binomial(logit)))
summary(glm(treatment~y, family = binomial(logit)))
summary(lm(y~treatment+x))

Here $X$ follows a mixture distribution: a different Gaussian in each treatment group (means $\bar x_1 \neq \bar x_2$), so that $P(S|X)\neq P(S)$. Then $Y$ is defined as an arbitrary (noisy) function $f$ of only $X$. Then we have the three equalities you wanted:
$E(Y|S=1)−E(Y|S=0)= E(f(X_1)) - E(f(X_2))= \theta$
$P(S|Y,X)=P(S|X)$ (because $X$ explains away $Y$)
$P(S|Y)≠P(S)$ (because distribution of $Y$ is different across treatment groups).
