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$
  • 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 '14 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 '14 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 '14 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 '14 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 '14 at 16:23
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$.

| cite | improve this answer | |
$\endgroup$
  • $\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 '14 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 '14 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 '14 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 '14 at 10:23
  • $\begingroup$ That's weird. Are you sure you're supplying a good starting point? $\endgroup$ – Sycorax Sep 12 '14 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.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ When you have a jackhammer, everything looks like a sidewalk. $\endgroup$ – shadowtalker Sep 8 '14 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 '14 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 '14 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 '14 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 '14 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.