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?
 A: 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] <- 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
```

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)
library(quadprog)

# 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.
