5
$\begingroup$

I would like to solve a linear regression (in R) with weights $w$ and a constraint.

In other words, I would like to find $x$ that minimizes the sum of squares $$\sum_i w_i(b_i-Ax_i)^2$$

On top of that I have an external vector $d$, which I would like to use in a constraint, such that $d \cdot x \le 5$.

Is this something that would be possible to do in R with solve.QP or perhaps some other R script?

Edit: I am adding a bounty for a solution that doesn't require any other custom software except the cran packages. While rstan works perfectly unfortunately I am unable to install it on my production servers due to old versions of some libraries.

$\endgroup$
7
  • 2
    $\begingroup$ Note that the weights are not a complication at all, because they can be absorbed in the values of $b$ and $A$, leading to an ordinary least squares problem with a single linear constraint. That means you problem is solved as described at stats.stackexchange.com/questions/24193 . $\endgroup$
    – whuber
    Sep 8, 2014 at 15:48
  • 1
    $\begingroup$ Very nearly the same question is asked--and answered in four different ways--at stats.stackexchange.com/questions/61733. It differs only in explicitly addressing the two dimensional case. $\endgroup$
    – whuber
    Sep 8, 2014 at 15:54
  • $\begingroup$ @user777 I am sorry about that: I was thinking that if these threads were close enough, we could merge your answer with the ones there. I liked your explicit demonstration of rstan, which has no parallel in the answers within the other threads. $\endgroup$
    – whuber
    Sep 8, 2014 at 16:02
  • $\begingroup$ @user777 I liked your solution as it was addressing this specific problem. If possible can you add it back? If anything it should be informative to someone else. $\endgroup$
    – Datageek
    Sep 8, 2014 at 16:03
  • 1
    $\begingroup$ There is no essential mathematical difference between bounding the slope and bounding the sum of coefficients: both are bounds on linear combinations of the coefficients. The solutions offered to those questions apply with very little change to your slightly more general formulation, thereby immediately giving you access to a variety of approaches to choose from. However, I have not voted to close your question, because although it does appear to be answered elsewhere, evidently it does take a little mathematical manipulation to see that those answers can apply. $\endgroup$
    – whuber
    Sep 8, 2014 at 16:23

2 Answers 2

3
+100
$\begingroup$

You're looking for the mgcv package. With the toy data we used before, it works just fine. (I'm uncertain why rstan is so confident in its results... I'm still looking into it.)

set.seed(1880)

N       <- 1500
d       <- c(1/2, 2/pi, 2/3)
x       <- c(2, 1, 3)
limit   <- 5
d%*%x <= limit

A       <- cbind(1, rnorm(N), rnorm(N))
b.hat   <- A%*%x
wgt     <- rexp(N)
b       <- rnorm(N, mean=b.hat, sd=wgt)

library(mgcv)

pin <- c(1.5, .75, 2.5)
Ain <- matrix(d, nrow=1)

M   <- list(y=b, w=wgt, X=A, p=pin, Ain=-Ain, bin=-limit, C=matrix(1, ncol=0, nrow=0))
pcls(M)

1.8844996 0.9421333 2.9770852

The inequality in this package is flipped the other direction by default. So we have to multiply both sides by $-1$.

$\endgroup$
5
  • $\begingroup$ This is absolutely brilliant, thank you very much @user777! It works flawlessly on my data. I will award bounty tomorrow, as 24h need to pass. Time for me to read "Generalized Additive Models: An Introduction with R" to keep up with your expertise. Hope you don't mind I selected this as a accepted answer instead of rstan. Rstan code is very useful nevertheless. $\endgroup$
    – Datageek
    Sep 10, 2014 at 18:31
  • 1
    $\begingroup$ Thanks! I'm really glad that this solves your problem. As for expertise, I follow CV daily, and it's always an education in how much I have yet to learn. On the other hand, I'm probably punching above the weight class for BA's in political science... ;) $\endgroup$
    – Sycorax
    Sep 10, 2014 at 18:33
  • $\begingroup$ As for rstan, there's honestly little reason to use it; I just wanted to see if I could (and the answer is somewhat mixed!). $\endgroup$
    – Sycorax
    Sep 10, 2014 at 19:49
  • $\begingroup$ I noticed a funny thing. On my data, if I set the same seed as in the above example I do get a correct solution. Otherwise my solution does not seem to converge and instead I am getting values very close to zero. $\endgroup$
    – Datageek
    Sep 12, 2014 at 10:23
  • $\begingroup$ That's weird. Are you sure you're supplying a good starting point? $\endgroup$
    – Sycorax
    Sep 12, 2014 at 12:30
6
$\begingroup$

Whenever I have a complicated model to fit, I usually just fit it directly in rstan because it's great at fitting highly constrained coefficients, and because it's easy to include penalties and transformations of variables. This is true even when I'm not explicitly fitting a Bayesian model.

This is what I've worked up for your particular problem.

library(rstan)

set.seed(1880)

N       <- 1500
d       <- c(1/2, 2/pi, 2/3)
x       <- c(2, 1, 3)
limit   <- 5
d%*%x <= limit
> TRUE
A       <- cbind(1, rnorm(N), rnorm(N))
b.hat   <- A%*%x
tau     <- 5
wgt     <- rexp(N)
Sigma   <- tau*wgt
b       <- rnorm(N, mean=b.hat, sd=Sigma)

constrained.reg <- "
    data{
        int<lower=1>        N;
        int<lower=1>        K;
        vector<lower=0>[N]  wgt;
        matrix[N,K]         A;
        vector[N]       b;
        vector[K]       d;
        real            limit; // s.t. d*x<=limit
    }
parameters{
    real<upper=limit>   c; // this is the largest possible value of x%*%d.
    simplex[K]      sim_x;
    real<lower=0>       tau;
}
transformed parameters {
    vector[K]   x;
    vector[N]   b_hat;
    vector[N]   Sigma;

    x       <- d .*sim_x /c;
    b_hat   <- A*x;
    Sigma   <- tau*wgt;
}
    model{
        b ~ normal(b_hat, Sigma);
        increment_log_prob(-2*log(tau)); // uniform prior on beta, noninformative prior on tau
    }
    generated quantities{
        vector[N]   resid;
        resid   <- (b_hat-b) ./Sigma;
    }
"
fake.data   <- list(N=N, A=A, K=3, b=b, wgt=wgt, d=d, limit=limit)

fit.test    <- stan(model_code=constrained.reg, data=fake.data, iter=10)

system.time(fit     <- stan(fit=fit.test, iter=1000, data=fake.data))
print(fit, c("x", "tau")); x

I realized that I was being dense and that we can enforce the inequality by sampling a value as large as the maximum permissible dot product result and then transforming appropriately.

     mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
x[1] 1.99       0 0.01 1.98 1.98 1.99 1.99  2.00  1645 1.00
x[2] 0.99       0 0.01 0.97 0.98 0.99 0.99  1.00   624 1.00
x[3] 3.00       0 0.01 2.98 2.99 3.00 3.01  3.02   945 1.00
tau  4.82       0 0.09 4.62 4.76 4.82 4.88  5.00   558 1.01

These results look fine to me.

$\endgroup$
12
  • 1
    $\begingroup$ When you have a jackhammer, everything looks like a sidewalk. $\endgroup$ Sep 8, 2014 at 19:24
  • 1
    $\begingroup$ I'm not exactly sure what @ssdecontrol is referring to :-), but as far as the other answers go, one (as I recall) suggested finding the weighted least squares solution. If it does not satisfy the constraint, then find the equality-constrained solution. Observe that (1) this requires only two least squares calculations, so it is efficient and fast and (2) in either case you obtain a covariance matrix for the estimates. Continuing in the analogical vein, that would make this problem about as easy as breaking a saltine into your soup, so you can put away the jackhammer :-). $\endgroup$
    – whuber
    Sep 8, 2014 at 19:30
  • 1
    $\begingroup$ @whuber I'm not sure I fully understand. What happens when neither solutions meets the constraints? $\endgroup$
    – Sycorax
    Sep 8, 2014 at 19:32
  • 3
    $\begingroup$ The second one, if it exists, is guaranteed to meet the constraints, by construction. Otherwise, there is no solution. (But there always will be when the weights are positive.) Here's an example: fit a line $y=ax+b$ to $(x,y)$ data with the constraint $a+b\le 5$. If the OLS estimates do not satisfy that constraint, now fit the model $y=ax+5-a$ to those data. (You could easily do that by fitting the no-intercept model $y-5=a(x-1)$ to the data $(x-1,y-5)$, for instance.) I hope it's obvious the constraint will be met and that no special software is needed. $\endgroup$
    – whuber
    Sep 8, 2014 at 19:54
  • 2
    $\begingroup$ No, I think it means something is not quite right about the code in the present solution. Perhaps I just am not interpreting its output correctly... $\endgroup$
    – whuber
    Sep 10, 2014 at 14:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.