# How to fit a constrained regression in R?

My regression: $$Y = α + β_1X_1 + β_2X_2 + β_3X_3 + β_4X_4 + β_5X_5 + β_6X_6 + β_7X_7 + β_8X_8 + β_9X_9 + β_{10}X_{10} + β_{11}X_{11} + β_{12}X_{12} + β_{13}X_{13} + β_{14}X_{14} + ε$$

Constraints that I need to set:

$$β_3 + β_4 + β_5 = 1 \tag{1}$$ $$β_6 + β_7 + β_8 = 0 \tag{2}$$ $$β_9 + β_{12} + β_{14} = 0 \tag{3}$$ $$β_{10} + β_{12} + β_{13} = 0 \tag{4}$$ $$β_{11} + β_{13} + β_{14} = 0 \tag{5}$$

I thought of adjusting the regression as a solution for constraining the coefficients:

(1) Take $$β_5$$ as $$(1 - β_3 - β_4)$$

(2) Take $$β_8$$ as $$(-β_6 - β_7)$$;

(3+4+5) Problem: If I use a system of equations and adapt the regression in order to estimate 3 of the implied 6 coefficients,then I would come up with multiple solutions depending on my choice of the coefficients to estimate directly

I also tried to use the ConsReg function from the foptim package but I it does not seem to work for equality constraints. Further, I have also been trying to use the constreg function from the coneproj package but unsucessfully since I am having difficulties in setting the constraint matrix. Does anyone have an idea or tip about how I can solve this constrained regression using R?

• I do not understand the sense in which you mean "multiple solutions," because (assuming your model is identifiable) there will be just one solution: all that can vary is how that solution is expressed.
– whuber
May 12, 2021 at 16:18
• I meant for example that if I select the bundles (Beta9, Beta11, Beta13) and (Beta10, Beta 11, Beta13), respectively, to perform the regression and from those coefficients derive (Beta10, Beta12, Beta14) and (Beta9, Beta12, Beta14) respectively, I get different values for some of the coefficients between estimations May 12, 2021 at 16:47

Roughly speaking: if in a linear regression the $$p$$ parameters are subject to $$m$$ (independent) linear constraints, the model is simply a reparameterisation of a linear model with $$p- m$$ unconstrained parameters.

Consider the linear regression in standard matrix form $$\tag{1} \mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}.$$ where $$\mathbf{X}$$ is $$n \times p$$. To stick to the OP notations, take $$\beta_0 = \alpha$$ and assume that the index $$i$$ of $$\beta_i$$ runs from $$0$$ to $$p-1$$. The system of linear constraints writes in matrix form as
$$\tag{2} \mathbf{A}^\top \boldsymbol{\beta} = \mathbf{a}$$

where $$\mathbf{A}$$ is a matrix with size $$p \times m$$ and $$\mathbf{a}$$ is a vector of length $$m$$. We can assume without loss of generality that $$m < p$$ and that $$\mathbf{A}$$ has rank $$m$$. We can complete the $$m$$ columns of $$\mathbf{A}$$ with the $$p-m$$ columns of a $$p\times(p-m)$$ matrix $$\mathbf{B}$$ so that the $$p \times p$$ matrix

$$\mathbf{C} = \begin{bmatrix} \mathbf{A} & \mathbf{B} \end{bmatrix}$$

is invertible. Then we can use $$\boldsymbol{\gamma} := \mathbf{C}^\top \boldsymbol{\beta}$$ as a new parameter vector for the regression model $$\mathbf{Y} = \mathbf{Z}\boldsymbol{\gamma} + \boldsymbol{\varepsilon}$$ where $$\mathbf{Z} := \mathbf{X}\mathbf{C}^{-\top}$$, which is equivalent to (1). Using blocks we can write $$\mathbf{Z} =[\mathbf{Z}_1, \, \mathbf{Z}_2]$$ where the two blocks have $$m$$ and $$p - m$$ columns, and consistently $$\boldsymbol{\gamma}$$ stacks $$\boldsymbol{\gamma}_1$$ and $$\boldsymbol{\gamma}_2$$. The constraint (2) becomes $$\boldsymbol{\gamma}_1 = \mathbf{a}$$ in the new parameterisation. So we just have to omit the $$m$$ first columns of $$\mathbf{Z}$$ and fit the unconstrained linear regression model with parameter $$\boldsymbol{\gamma}_2$$ and an offset

$$\mathbf{Y} = \underset{\text{offset}}{\underbrace{\mathbf{Z}_1 \mathbf{a}}} + \mathbf{Z}_2\boldsymbol{\gamma}_2 + \boldsymbol{\varepsilon},$$

and then complete $$\widehat{\boldsymbol{\gamma}}_2$$ with $$\widehat{\boldsymbol{\gamma}}_1 = \mathbf{a}$$ to get the full $$\widehat{\boldsymbol{\gamma}}$$. The estimate $$\widehat{\boldsymbol{\beta}}$$ for the original vector $$\boldsymbol{\beta}$$ comes as $$\widehat{\boldsymbol{\beta}} = \mathbf{C}^{-\top} \widehat{\boldsymbol{\gamma}}$$. We can infer on $$\boldsymbol{\gamma}_2$$ as usual and derive predictions as well.

Now I shortly explain how to get $$\mathbf{B}$$ in practice. We can use the QR decomposition of the $$p \times m$$ matrix $$\mathbf{A}$$

$$\mathbf{A} = \underset{p \times p }{\underbrace{\mathbf{Q}}} \overset{p \times m }{\overbrace{\mathbf{R}}} = \begin{bmatrix} \mathbf{Q}_1 \, \mathbf{Q}_2 \end{bmatrix} \begin{bmatrix} \mathbf{R}_1 \\ \mathbf{0} \end{bmatrix} = \mathbf{Q}_1 \mathbf{R}_1$$

where $$\mathbf{Q}$$ is an orthogonal matrix and $$\mathbf{R}$$ is upper triangular with non-zero diagonal elements. The blocks $$\mathbf{Q}_1$$ and $$\mathbf{Q}_2$$ have $$m$$ and $$p-m$$ columns, whith $$\mathbf{R}_1$$ is $$m \times m$$ and invertible. We can take $$\mathbf{B} := \mathbf{Q}_2$$. Indeed, since $$\mathbf{A} = \mathbf{Q}_1\mathbf{R}_1$$, the space $$\text{cspan}(\mathbf{A})$$ generated by the columns of $$\mathbf{A}$$ in $$\mathbb{R}^p$$ is the same as $$\text{cspan}(\mathbf{Q}_1)$$, with $$\text{cspan}(\mathbf{Q}_2)$$ being a complement subspace.

For a R implementation, you need the design matrix $$\mathbf{X}$$ and the matrix of restrictions $$\mathbf{A}$$. The qr.Q function works on an object created by qr. It computes by default only the thin QR decomposition $$\mathbf{A} = \mathbf{Q}_1 \mathbf{R}_1$$, and we have to use complete = TRUE to get the whole $$p \times p$$ matrix $$\mathbf{Q}$$. Of course we could use the available QR decomposition to avoid inverting $$\mathbf{C}$$ or its transpose.

## hypothetical design matrix
set.seed(135)
n <- 30; p <- 15; m <- 5
X <- cbind(1, matrix(runif(n * (p-1)), nrow = n))

## Constraints: mind the shift in the row indices due to the constant 'beta0'
A <- matrix(0.0, nrow = p, ncol = m)
a <- rep(0.0, m)
A[4, 1] <- A[5, 1] <- A[6, 1] <- 1.0
A[7, 2] <- A[8, 2] <- A[9, 2] <- 1.0
A[10, 3] <- A[13, 3] <- A[15, 3] <- 1.0
A[11, 4] <- A[13, 4] <- A[14, 4] <- 1.0
A[12, 5] <- A[14, 5] <- A[15, 5] <- 1.0
a <- 1.0

## linear algebra
ind1 <- 1:m
ind2 <- (m + 1):p
QR <- qr(A)
Q <- qr.Q(QR, complete = TRUE)
B <- Q[ , ind2, drop = FALSE]
C <- cbind(A, B)
Cinv <- solve(C)
Z <- X %*% t(Cinv)
offset <- Z[ , ind1, drop = FALSE] %*% a

## Draw true coef and random observations
gamma2 <- rnorm(p - m)
Y <- Z %*% c(a, gamma2) + rnorm(n, sd = 0.1)
fit <- lm(Y ~ Z[ , ind2, drop = FALSE]  - 1, offset = offset)
gammaHat <- c(a, coef(fit))
betaHat <- solve(t(C), gammaHat)

## check
t(A) %*% betaHat - a
$$$$


Maybe the R package pracma, which uses quadprog under the hood, can help.

They are implementations of quadratic programming. In your case the syntax should be something like:

library(pracma)

# Matrix of equality constraints:
Aeq <- rbind(
c(0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),
c(0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0),
c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0),
c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1)
)

# Vector of constraints
beq <- c(1, 0, 0, 0, 0)


Solve:

fit <- lsqlincon(X, y, Aeq= Aeq, beq= beq)


where X is the design matrix of n rows and 14 columns and y is the vector of n responses. fit` are the n coefficients.

• This is an inferior solution because it doesn't offer the usual statistically important capabilities of providing standard errors of estimate, etc. Why not just use the five equations to rewrite the regression in terms of $14-5=9$ variables? By inspection, you can eliminate each of the first variables appearing in the five constraints.
– whuber
May 12, 2021 at 21:35
• @whuber yes - I'm aware of those limitations. I don't mind to delete this answer but I think it is still an option for cases where one cannot or doesn't know how to rewrite the model. May 13, 2021 at 7:48