How do I fit a constrained logistic regression model via quadratic programming in R?

Question

I trying to find $\pi_{1}, \pi_{2}, \pi_{3}$ for model: $$ Y = \pi_{1}X_{1} + \pi_{2}X_{2} + \pi_{3}X_{3} + \epsilon, $$ with constraints: $\Sigma_{k}\pi_{k}=1$ and $\pi_{k}\geq0$. (All $\pi$ are positive and add up to 1)

I found a very similar question here that tried to minimize squared errors: $\Sigma_{i}(Y_{i}-(\pi_{1}X_{i1} + \pi_{2}X_{i2} + \pi_{3}X_{i3}))^{2}$.

However, my dependent variable $Y$ is binary, not continuous. So it's like fitting a logistic regression model rather than a linear regression model.

Is there a way for me to apply a link function (e.g. logit function) to the solution provided in the above link to find better estimates of $\pi_{1}, \pi_{2}, \pi_{3}$ ?

Or, are there any other well-known solutions to solve problems like these (quadratic programming where dependent variable is binary?)

I copied the example code from the above link with the dependent variable rounded to either 1 or 0 in case it helps.

library("quadprog");
X <- matrix(runif(300), ncol=3)
Y <- round(X %*% c(0.2,0.3,0.5))
Rinv <- solve(chol(t(X) %*% X));
C <- cbind(rep(1,3), diag(3))
b <- c(1,rep(0,3))
d <- t(Y) %*% X  
solve.QP(Dmat = Rinv, factorized = TRUE, dvec = d, Amat = C, bvec = b, meq = 1)
# 0.02080963 0.18475508 0.79443529

I tried searching for similar problems on Google, but with my little knowledge on quadratic programs, I couldn't get close to finding information that matched mine.

Answer 1

The reason you are unable to reproduce the coefficients in your example problem is that the data generating model you coded doesn't correspond to the model you fit. Replacing the line

Y <- round(X %*% c(0.2,0.3,0.5))

with

Y <- rbinom(300, 1, X %*% c(0.2,0.3,0.5))

and running the rest does recover (close to) the true values. The former creates a strict threshold at .5, above which all Ys are 1 and below which all are 0. This is not included in your optimization. The fact that you have a binary outcome doesn't mean anything; you can still use the linear probability model with constrained least squares as you did to compute the coefficients.

That said, if you want to do constrained logistic regression, you can use CVXR, which is a general purpose optimization engine, to specify a logistic regression model with constraints, as follows:

library(CVXR)
beta <- Variable(3)
obj <- -sum((1-Y) * (X %*% beta)) - sum(logistic(-X %*% beta))
prob <- Problem(Maximize(obj), constraints = list(sum(beta) == 1, beta >= 0))
result <- solve(prob)
result$getValue(beta)

This relies on maximizing the log likelihood assuming a binomial family and logistic link. This is not a quadratic optimization problem because the objective function is not quadratic. Including these constraints is not useful for this problem because the resulting predictions are already bound between 0 and 1 so the coefficients can be of any value, and you will not be able to recover the data-generating coefficients you specified unless you change them to have a logistic form, e.g.,

Y <- rbinom(300, 1, plogis(X %*% c(0.2,0.3,0.5)))