Logistic regression with known probabilities for some datapoints I have a dataset with several features and a binary outcome (0 or 1), from which I want to estimate the probability of getting 1 based on some of those features.
I thought of applying logistic regression, however there’s a constraint I have on some datapoints where the probability is expected to be around 1/2. How do I enforce this in the model training? (Related question: how do I enforce this with sklearn?)
Intuitively I’m thinking of modifying the usual likelihood function by adding the terms with probability at 1/2, but not sure if there’s a better way or if I should use another approach altogether.
 A: I will assume that the objective here is to fit the model subject to the constraint that the predicted probability for one of the observations should have a value of exactly $p_0$ (not just a value "around" $p_0$).
Letting $\mathbf{x}_0$ denote the covariate vector for that observation, this amounts to a linear constraint
$$
\mathbf{x}_0^T\boldsymbol\beta=\eta_0=\operatorname{logit} p_0 \tag{1}
$$
on the parameter vector $\boldsymbol\beta$.  This means that one of the elements of $\boldsymbol\beta$ in effect is redundant and can be eliminated from the model by substitution.  This in turn leads to a modified design matrix and an offset term as follows.
Partitioning $\boldsymbol\beta$ into subvectors containing the redundant element $\beta_j$ and the remaining elements $\boldsymbol\beta_{-j}$, and partitioning $\mathbf{x}_0$ the same way such that $x_{0,j}\neq 0$, (1) can be written as
$$
\begin{bmatrix} \mathbf{x}_{0,-j}^T & x_{0,j}\end{bmatrix}
\begin{bmatrix} \boldsymbol\beta_{-j}\\ \beta_j \end{bmatrix}
=\eta_0
$$
or
$$
\mathbf{x}_{0,-j}^T \boldsymbol\beta_{-j} + x_{0,j}\beta_j=\eta_0
$$
implying that
$$
\beta_j=x_{0,j}^{-1}(\eta_0 - \mathbf{x}_{0,-j}^T \boldsymbol\beta_{-j}).\tag{2}
$$
For the other observations we have
$$
\boldsymbol\eta=\mathbf{X}\boldsymbol\beta,
$$
or, after partitioning of the design matrix $\mathbf{X}$ the same way into its $j$'th column $\mathbf{x}_j$ and the remaining columns $\mathbf{X}_{-j}$,
$$
\boldsymbol\eta=\begin{bmatrix} \mathbf{X}_{-j} & \mathbf{x}_j\end{bmatrix}
\begin{bmatrix} \boldsymbol\beta_{-j}\\ \beta_j \end{bmatrix}=\mathbf{X}_{-j}\boldsymbol\beta_{-j} + \mathbf{x}_j\beta_j. \tag{3}
$$
Substituting (2) into (3) and rearranging terms the model takes the form
$$
\boldsymbol\eta=\underbrace{(\mathbf{X}_{-j}-x_{0,j}^{-1}\mathbf{x}_j\mathbf{x}_{0,-j}^T)}_{\mathbf{X}^*}\boldsymbol\beta_{-j} + \underbrace{x_{0,j}^{-1}\eta_0\mathbf{x}_j }_{\text{offset}}
$$
where $\mathbf{X}^*$ is the appropriate modified design matrix and the second term the offset term needed to fit the model subject to (1).
As a side note, a potentially useful application of the above method is the computation of profile likelihood based confidence intervals for $\eta_0$ and $p_0$ which may have better performance than the current confidence intervals computed by predict.glm in R that relies on asymptotic normality of $\hat{\boldsymbol\beta}$.  The above method can also be straightforwardly extended to fit GLMs under general linear hypotheses on the form $\mathbf{C}\boldsymbol\beta=\mathbf{d}$ to facilitate likilihood ratio tests as opposed to sometimes less accurate Wald tests of such hypotheses.
