Solution (source code)

= Solution

Add a sex-by-treatment <interaction term>:
$$
\operatorname{logit}(p_i)=\beta_0+\beta_S M_i+\beta_TT_i+\beta_{ST}M_iT_i.
$$
In R the formula becomes `prop ~ sex * treatment`. There are four parameters for four cell probabilities, so this model is saturated and, provided the cell proportions lie strictly between zero and one, fits each observed proportion exactly. It therefore recovers all the stated order relations.