15

I would like to manually fix a certain coefficient, say $\beta_1=1.0$, then fit coefficients to all other predictors, while keeping $\beta_1=1.0$ in the model.

How can I achieve this using R? I'd particularly like to work with LASSO (glmnet) if possible.

Alternatively, how can I restrict this coefficient to a specific range, say $0.5\le\beta_1\le1.0$?

Nick Stauner
  • 11,558
  • 5
  • 47
  • 105
raco
  • 339
  • 1
  • 3
  • 11
  • To specify box constraints on the fitted coefficients there's the arguments lower.limits and upper.limits in glmnet, right? – Tom Wenseleers Sep 21 '17 at 09:11

4 Answers4

11

Well, let's think. You have:

$Y = b_0 + b_1x_1 + b_2x_2 + e$

(to keep it simple) You want to force $b_1 = 1$ so, you want

$Y = b_0 + x_1 + b_2x_2 + e$

so you could just subtract $x_1$ from each side leaving:

$Ynew = Y-x_1 = b_0 + b_2x_2 + e$

which can then estimate $b_2$.

Peter Flom
  • 94,055
  • 35
  • 143
  • 276
  • 2
    That's the easy part (and has been addressed in other threads as I recall). What about restricting the coefficient to a range? The especially hard part of this problem is obtaining good confidence limits when the estimate lies on the boundary of the constraint region. – whuber Jan 16 '14 at 00:58
  • 2
    That is definitely harder. I missed the end of the post. But I think I should leave my answer up as it does answer part of the question – Peter Flom Jan 16 '14 at 01:10
  • Does this still generalize if $\beta_1 \neq 1$? Let $\beta_1 = 0.75$ instead of 1, $Ynew = Y - .75x_1 = \beta_0 + (\beta_1\prime - 0.75)x_1 + \beta_2x_2 + \epsilon$, where $\beta_1\prime$ is the coefficient chosen by OLS regression. – raco Jan 16 '14 at 17:25
  • 1
    Yes, if it is fixed at .75 then doing what you say will work. But as @whuber points out, that's the easy part of this problem – Peter Flom Jan 16 '14 at 17:52
  • 2
    @whuber, in a Bayesian framework you could throw in a Metropolis step to toss out any coefficients outside your range or alternately you could sample from a truncated multivariate normal distribution. – John Jan 16 '14 at 21:41
  • If you could flesh out that idea it might make a great answer, @John. – whuber Jan 16 '14 at 21:45
  • Can you do that kind of substractions? The variables can be correlated, and you may be not taking into account that transformation to calculate the distribution of the residuals. Furthermore, from the glm help: 'The null model will include the offset, and an intercept if there is one in the model. Note that this will be incorrect if the link function depends on the data other than through the fitted mean: specify a zero offset to force a correct calculation' – skan May 14 '17 at 20:02
  • To specify box constraints on the fitted coefficients there's the arguments lower.limits and upper.limits in glmnet, right? – Tom Wenseleers Sep 21 '17 at 09:12
6

You need to use the offset argument like this:

library(glmnet)
x=matrix(rnorm(100*20),100,20)
x1=matrix(rnorm(100),100,1)
y=rnorm(100)
fit1=glmnet(x,y,offset=x1)
fit1$offset
print(fit1)

About the range ... I don't think that has been implemented in glmnet. If they use some numerical method, you may want to dig into the R code and try to restrict it over there, but you'll need a good, solid programming background.

Stat
  • 7,078
  • 1
  • 24
  • 49
  • 2
    What is `offset` actually doing? How is the value of `1.1*x1` determined from the question? – whuber Jan 16 '14 at 14:12
  • I read the documentation for 'offset' in glmnet, and I'm still not sure what it does. I could not find any great examples, but most reference Poisson processes. Why is 1.1*x1 used? – raco Jan 16 '14 at 19:20
  • I thought he is fixing the coefficients to be $\beta_1=1.1$. I just edited the answer. The offset is the term in which its coefficient is not estimated by the model but is assumed to have the value 1. – Stat Jan 16 '14 at 21:21
  • I'm happy enough with this answer. I can iterate over different offset "coefficients" & compare models. Thanks! – raco Jan 22 '14 at 23:16
  • 1
    With regards to `offset` in the `glmnet` package, the answer provided by Stat does not make sense me. When I run fit1$beta[,ncol(fit1$beta)] I do not see any $\beta_1=1.0$. Could you clarify how offset is working in your example? For the range of the betas, you can use the `lower.limits` and `upper limits` arguments. – Mario Nuñez Jan 14 '16 at 18:59
  • To specify box constraints on the fitted coefficients there's the arguments lower.limits and upper.limits in glmnet, right? – Tom Wenseleers Sep 21 '17 at 09:11
  • Here is an explanation of `offset`, but no reference provided: https://statisticaloddsandends.wordpress.com/2019/01/09/a-deep-dive-into-glmnet-offset/ – mt1022 Apr 27 '21 at 13:14
3

With respect to constraining coefficients to be within a range, a Bayesian approach to estimation is one means to accomplish this.

In particular, one would rely on a Markov Chain Monte Carlo. First, consider a Gibbs sampling algorithm, which is how you would fit the MCMC in a Bayesian framework absent the restriction. In Gibbs sampling, in each step of the algorithm you sample from the posterior distribution of each parameter (or group of parameters) conditional on the data and all other parameters. Wikipedia provides a good summary of the approach.

One way to constrain the range is to apply a Metropolis-Hastings step. The basic idea is to simply throw out any simulated variable that is outside of your bounds. You could then keep re-sampling until that is within your bounds before moving on to the next iteration. The downside to this is that you might get stuck simulating a lot of times, which slows down the MCMC. An alternative approach, originally developed by John Geweke in a few papers and extended in a paper by Rodriguez-Yam, Davis, Sharpe is to simulate from a constrained multivariate normal distribution. This approach can handle linear and non-linear inequality constraints on parameters and I've had some success with it.

John
  • 2,117
  • 16
  • 24
  • To specify box constraints on the fitted coefficients there's the arguments lower.limits and upper.limits in glmnet, right? – Tom Wenseleers Sep 21 '17 at 09:12
  • @TomWenseleers I was answering more generally. Look to some of the other answers with regards to `glmnet`. – John Sep 21 '17 at 13:56
2

I'm not familiar with LASSO or glmnet, but lavaan (short for "latent variable analysis") facilitates multiple regression models with both equality constraints and single-bounded inequality constraints (see the table on page 7 of this PDF, "lavaan: An R package for structural equation modeling"). I don't know if you could have both upper and lower bounds on the coefficient, but maybe you could add each bound with separate lines, e.g:

Coefficient>.49999999
Coefficient<1.0000001

Of course, if you're standardizing everything before fitting the model, you shouldn't have to worry about imposing an upper bound of 1 on your regression coefficients anyway. I'd say you're better off omitting it in this case, just in case something goes wrong! (lavaan is still in beta after all...I've seen some slightly fishy results in my own limited use of it thus far.)

Nick Stauner
  • 11,558
  • 5
  • 47
  • 105
  • To specify box constraints on the fitted coefficients there's the arguments lower.limits and upper.limits in glmnet, right? – Tom Wenseleers Sep 21 '17 at 09:12