18
$\begingroup$

I want to perform a very simple linear regression in R. The formula is as simple as $y = ax + b$. However I would like the slope ($a$) to be inside an interval, let's say, between 1.4 and 1.6.

How can this be done?

$\endgroup$
25
$\begingroup$

I want to perform ... linear regression in R. ... I would like the slope to be inside an interval, let's say, between 1.4 and 1.6. How can this be done?

(i) Simple way:

  • fit the regression. If it's in the bounds, you're done.

  • If it's not in the bounds, set the slope to the nearest bound, and

  • estimate the intercept as the average of $(y - ax)$ over all observations.

(ii) More complex way: do least squares with box constraints on the slope; many optimizaton routines implement box constraints, e.g. nlminb (which comes with R) does.

Edit: actually (as mentioned in the example below), in vanilla R, nls can do box constraints; as shown in the example, that's really very easy to do.

You can use constrained regression more directly; I think the pcls function from the package "mgcv" and the nnls function from the package "nnls" both do.

--

Edit to answer followup question -

I was going to show you how to use it with nlminb since that comes with R, but I realized that nls already uses the same routines (the PORT routines) to implement constrained least squares, so my example below does that case.

NB: in my example below, $a$ is the intercept and $b$ is the slope (the more common convention in stats). I realized after I put it in here that you started the other way around; I'm going to leave the example 'backward' relative to your question, though.

First, set up some data with the 'true' slope inside the range:

 set.seed(seed=439812L)
 x=runif(35,10,30)
 y = 5.8 + 1.53*x + rnorm(35,s=5)  # population slope is in range
 plot(x,y)
 lm(y~x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
     12.681        1.217  

... but LS estimate is well outside it, just caused by random variation. So lets use the constrained regression in nls:

 nls(y~a+b*x,algorithm="port",
   start=c(a=0,b=1.5),lower=c(a=-Inf,b=1.4),upper=c(a=Inf,b=1.6))

Nonlinear regression model
  model: y ~ a + b * x
   data: parent.frame()
    a     b 
9.019 1.400 
 residual sum-of-squares: 706.2

Algorithm "port", convergence message: both X-convergence and relative convergence (5)

As you see, you get a slope right on the boundary. If you pass the fitted model to summary it will even produce standard errors and t-values but I am not sure how meaningful/interpretable these are.

So how does my suggestion (1) compare? (i.e. set the slope to the nearest bound and average the residuals $y-bx$ to estimate the intercept)

 b=1.4
 c(a=mean(y-x*b),b=b)
       a        b 
9.019376 1.400000

It's the same estimate ...

In the plot below, the blue line is least squares and the red line is the constrained least squares:

constrained and LS line

$\endgroup$
  • $\begingroup$ Thank you for this answer but... could you give an example using any of these functions? $\endgroup$ – Iñigo Hernáez Corres Jun 14 '13 at 10:39
  • 1
    $\begingroup$ +1 Finding confidence intervals on the parameter estimates is going to be a challenge in any event. $\endgroup$ – whuber Jun 14 '13 at 12:36
  • $\begingroup$ @IñigoHernáezCorres see the update to my answer, where I illustrate using nls to do it. $\endgroup$ – Glen_b -Reinstate Monica Jun 16 '13 at 2:25
  • $\begingroup$ +1 great answer with connections on two ways of doing it! $\endgroup$ – Haitao Du Apr 19 '17 at 14:06
15
$\begingroup$

Glen_b's second method, using least squares with a box constraint can be more easily implemented via ridge regression. The solution to ridge regression can be viewed as the Lagrangian for a regression with a bound on the magnitude of the norm of the weight vector (and hence its slope). So following whuber's suggestion below, the approach would be to subtract a trend of (1.6+1.4)/2 = 1.5 and then apply ridge regression and gradually increase the ridge parameter until the magnitude of the slope is less than or equal to 0.1.

The benefit of this approach is that no fancy optimisation tools are required, just ridge regresson, which is already available in R (and many other packages).

However Glen_b's simple solution (i) seems sensible to me (+1)

$\endgroup$
  • 5
    $\begingroup$ This is clever, but are you sure it will work as described? It seems to me the appropriate approach would be to remove a trend of (1.6+1.4)/2 = 1.5 and then control the ridge parameter until the absolute value of the slope is less than or equal to 0.1. $\endgroup$ – whuber Jun 14 '13 at 12:39
  • 1
    $\begingroup$ yes, that is indeed a better suggestion. The ridge regression approach is really more appropriate if the restriction is on the magnitude of the slope, it sounds like a pretty odd problem! My answer was originally inspired by Glen_b's comment on box constraints, ridge regression is basically just an easier way to implement the box constraints. $\endgroup$ – Dikran Marsupial Jun 14 '13 at 12:52
  • $\begingroup$ Although I appreciate your acknowledgment of my comments, imho it distracts from the content of your answer. We're all in this together to improve our work wherever we can, so it's enough acknowledgment that you acted on my suggestions. For that you deserve the increment in reputation. If you are moved to make additional edits, please consider streamlining the text by removing that superfluous material. $\endgroup$ – whuber Jun 14 '13 at 14:32
  • $\begingroup$ Superfluous material edited, however I enjoy collaborations and always seek to give collaborators the credit they deserve, and still think morally you deserve half the up-votes. ;o) $\endgroup$ – Dikran Marsupial Jun 14 '13 at 15:29
10
$\begingroup$

Another approach would be to use Bayesian methods to fit the regression and choose a prior distribution on $a$ that only has support in the region you want, e.g. a uniform from 1.4 to 1.6, or a beta distribution shifted and scaled to that domain.

There are many examples on the web and in software of using Bayesian methods for regression, you can just follow one of those examples and change the prior on $a$.

This result will still give credible intervals of the parameters of interest (of course the meaningfulness of these intervals will be based on the reasonableness of your prior information about the slope).

$\endgroup$
  • $\begingroup$ +1, this was my first thought as well. I like the other suggestions, but this seems the best to me. $\endgroup$ – gung - Reinstate Monica Jun 16 '13 at 21:05
0
$\begingroup$

Another approach might be to reformulate your regression as an optimization problem and use an optimizer. I'm not sure if it can be reformulated this way, but I thought of this question when I read this blog posting on R optimizers:

http://zoonek.free.fr/blosxom/R/2012-06-01_Optimization.html

$\endgroup$

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.