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

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.

• just commenting that you can avoid a constrained problem by passing variables $\pi_i$ through the softmax function (which, frankly, might be easiest in R). Aug 23, 2022 at 16:45
• How many observations? How many predictors, only 3? Just write the likelihood function in R, reparametrize so as to avoid restrictions, and pass it to optim() Aug 23, 2022 at 16:55
• @JohnMadden do you mean that I should fit a logistic regression model (via glm()) and pass the coefficients through the softmax function? Aug 23, 2022 at 16:55
• @kjetilbhalvorsen There are around 100K+ observations and only 3 predictors. I will look into your suggestion and the optim() function. Thank you! Aug 23, 2022 at 16:57
• You should be using a logistic regression rather than least squares. But the least squares problem has an easy solution, as sketched at stats.stackexchange.com/questions/51194 (or, even easier, by just trying all four versions of the model determined by which constraints are bound). That's because you can write it in the form $$Y - X_3 = \pi_1(X_1-X_3) + \pi_2(X_2-X_3) + \epsilon.$$
– whuber
Aug 23, 2022 at 18:24

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)))

• Thank you! Your explanations and code were all extremely helpful. Also, thanks for clarifying that I can still "use the linear probability model with constrained least squares" even if I have binary outcomes. Aug 24, 2022 at 14:40