How to fit this linear regression with constraints?

Question

I would like to fit the following model:

$$Y=\beta_0+\beta_1(\sum_{i=1}^kw_iX_i)+\beta_2(\sum_{i=1}^kw_iX_i)^2+\epsilon$$ where $\beta_0, \beta_1, \beta_2, w_1,...,w_k$ are the parameters, and $\epsilon$ is some normal noise. It does not look like something I have encountered before, and I know it's different from simply including all the second order and interaction terms, since the coefficients are related/fixed in a specific way through the sharing of the weights $w_i$'s. It seems that this is some sort of linear regression with constraints that relate the coefficients. Could someone point me in the right direction how such model can be fitted?

Answer 1

The model is overparametrised: you don't need $\beta_1$, which can be set to anything convenient, like 1.

One thing I thought of was to fit iteratively. Start out with some guess at $w$ and $\beta_2$. Then compute $Z=(\sum_i \hat{w}_iX_i)^2$ and fit the linear model

Y~ X1+X2+...+X_k + Z

The coefficients of the $X$s are the new $\hat{w}_i$, and the coefficient of $Z$ is $\hat \beta_2$. And then recompute Z, iterate and hope it converges. Sadly, it doesn't.

But if $k$ isn't too large, it's easy to just compute the residual sum of squares as a function of the parameters and run it through a general purpose optimiser. In R I'd use minqa::newuoa, but there are lots of alternatives.

> X<-matrix(rnorm(50*100),ncol=5)
> w<-1:5
> Y<- (X%*%w)+2*(X%*%w)^2+rnorm(100)
> 
> 
> rss<-function(theta){
+   beta2<-theta[1]
+   w<-theta[-1]
+   mu<- (X%*%w)+beta2*(X%*%w)^2
+   sum((Y-mu)^2)
+ }
> 
> minqa::newuoa(par=rep(1,6), rss)
parameter estimates: 1.99478699135839, 1.00032043499982, 2.00140284432351, 3.00312315850919, 4.00284240744153, 5.00537517104468 
objective: 1047.51402563294 
number of function evaluations: 1689 

Then use the bootstrap to get standard error estimates.

With $k=50$ it doesn't work (without tuning -- I'm sure it would work if the optimiser defaults were changed or the starting values were better)

Answer 2

If you write out the expression, you get a polynomial in terms of $X_1,X_2,..,X_k$, including their interactions, where the new "coefficients" are all function of $\beta$s and $w$s and twos. For k=2, you get a polynomial that has 5 coefficients (or 6 including the intercept) with 4 unknowns:

$$ \begin{align*} Y &= \beta_0+(\beta_1w_1)X_1+(\beta_1w_2)X_2+(\beta_2w_1^2)X_1^2 + (\beta_2 w_2^2)X_2^2+(2\beta_2 w_1w_2)X_1X_2 +\varepsilon \\ &= \alpha_0+\alpha_1X_1+\alpha_2X_2+\alpha_3X_1^2 + \alpha_4X_2^2+\alpha_5X_1X_2 +\varepsilon \end{align*} $$

If you fit this regression, you will get the new $\alpha$ coefficients, which gives you a system of non-linear equations:

$$ \begin{align*} \alpha_0 &= \beta_0 \\ \alpha_1 &= \beta_1w_1 \\ \alpha_2 &= \beta_1w_2 \\ \alpha_3 & =\beta_2w_1^2\\ \alpha_4 &= \beta_2 w_2^2 \\ \alpha_5 &= 2\beta_2 w_1w_2 \end{align*} $$

In principle, that system of equations should be solvable numerically, at least sometimes. It should remain solvable with $k>3$ since you don't have the curse of dimensionality since each new variable adds only one parameters but multiple new equations that help pin it down.

Here's a toy $k=2$ simulation example using Stata where I ignore the intercept equation since it is trivial:

. clear

. set obs 1000
number of observations (_N) was 0, now 1,000

. set seed 10011979

. gen b0 = 1 

. gen b1 = 2 

. gen b2 = 3

. gen w1 = 4 

. gen w2 = 5

. gen x1  = rnormal(0,1)

. gen x2  = rnormal(10,2)

. gen eps = rnormal()

. gen y = b0 + b1*(w1*x1 + w2*x2) + b2*(w1*x1 + w2*x2)^2 + eps

. reg y (c.x1 c.x2)##(c.x1 c.x2)

      Source |       SS           df       MS      Number of obs   =     1,000
-------------+----------------------------------   F(5, 994)       >  99999.00
       Model |  1.1237e+10         5  2.2475e+09   Prob > F        =    0.0000
    Residual |  1052.11816       994  1.05846897   R-squared       =    1.0000
-------------+----------------------------------   Adj R-squared   =    1.0000
       Total |  1.1237e+10       999  11248523.6   Root MSE        =    1.0288

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          x1 |   8.082131   .1573906    51.35   0.000     7.773275    8.390987
          x2 |   9.852645    .110114    89.48   0.000     9.636562    10.06873
             |
   c.x1#c.x1 |    47.9813   .0233895  2051.40   0.000      47.9354     48.0272
             |
   c.x1#c.x2 |   119.9907   .0153233  7830.59   0.000     119.9606    120.0208
             |
   c.x2#c.x2 |   75.00664   .0053927  1.4e+04   0.000     74.99605    75.01722
             |
       _cons |    1.77947   .5532575     3.22   0.001      .693783    2.865156
------------------------------------------------------------------------------

. 
. clear mata      

. mata:
------------------------------------------------- mata (type end to exit) -----------------------------------------------------------------------------------------------------------------------------------------------
: void mysolver(todo, p, lnf, S, H)
>          {
>                  b1   = p[1]
>                  b2   = p[2]
>                  w1   = p[3]
>                                  w2   = p[4]                 
>                  lnf = (b1*w1 - 8.082131)^2\   
>                        (b1*w2 - 9.852645)^2\
>                                            (b2*w1^2 - 47.9813)^2\
>                                            (b2*w2^2 - 75.00664)^2\
>                                            (2*b2*w1*w2 - 119.9907)^2
>                 }
note: argument todo unused
note: argument S unused
note: argument H unused

: 
: S = optimize_init()

: optimize_init_evaluator(S, &mysolver())

: optimize_init_evaluatortype(S, "v0")

: optimize_init_params(S, (1,1,1,1))

: optimize_init_which(S,  "min" )

: optimize_init_tracelevel(S,"none")

: optimize_init_conv_ptol(S, 1e-16)

: optimize_init_conv_vtol(S, 1e-16)

: p = optimize(S)

: p 
                 1             2             3             4
    +---------------------------------------------------------+
  1 |    2.1561597   3.521534782   3.691630188   4.614939185  |
    +---------------------------------------------------------+

: end
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

The solution is not very good (unless you squint and round to the nearest integer), since $p = (2,3,4,5)$ in the simulation. I am probably doing something wrong when when I solve the equations numerically. But even the intercept is pretty off with $b_0 = 1.77947 \ne 1$.

---

Code:

cls
clear
set obs 1000
set seed 10011979
gen b0 = 1 
gen b1 = 2 
gen b2 = 3
gen w1 = 4 
gen w2 = 5
gen x1  = rnormal(0,1)
gen x2  = rnormal(10,2)
gen eps = rnormal()
gen y = b0 + b1*(w1*x1 + w2*x2) + b2*(w1*x1 + w2*x2)^2 + eps
reg y (c.x1 c.x2)##(c.x1 c.x2)

clear mata  
mata:
void mysolver(todo, p, lnf, S, H)
         {
                 b1   = p[1]
                 b2   = p[2]
                 w1   = p[3]
                 w2   = p[4]                 
                 lnf = (b1*w1 - 8.082131)^2\   
                       (b1*w2 - 9.852645)^2\
                       (b2*w1^2 - 47.9813)^2\
                       (b2*w2^2 - 75.00664)^2\
                       (2*b2*w1*w2 - 119.9907)^2
        }

S = optimize_init()
optimize_init_evaluator(S, &mysolver())
optimize_init_evaluatortype(S, "v0")
optimize_init_params(S, (1,1,1,1))
optimize_init_which(S,  "min" )
optimize_init_tracelevel(S,"none")
optimize_init_conv_ptol(S, 1e-16)
optimize_init_conv_vtol(S, 1e-16)
p = optimize(S)
p 
end