I have some domain knowledge I want to use in a regression problem.

Problem statement

The dependent variable $y$ is continuous.
The independent variables are $x_1$ and $x_2$.

  • Variable $x_1$ is continuous and positive.
  • Variable $x_2$ is categorical, ordered, and takes only few different values (i.e, less than 10).

I want to fit a regression function $f$ so that $y = f(x_1,x_2)$, with the constraints

  1. $f(x_1,x_2)$ is monotonically increasing in $x_1$
  2. $f(x_1,x_2)$ is bounded in $[0,1]$
  3. $f(x_1,x_2)$ is "smooth" in $x_2$

These constraints come from domain knowledge of the problem.

Samples are evenly distributed among $x_2$ but not $x_1$.

My question: which techniques do you recommend for such regression problem?

I'm currently ignoring the last two constraints, and using the monreg package. (I run one regression for each possible value of $x_2$)

I can not give a formal definition of "smooth" in this context. I can assume that $f(x_1,x_2)$ does not change much along $x_2$ values that are consecutive.

I have found some SO questions regarding this issue but it looks like they have not raised much attention, or are very focused on specific R packages. Q1 Q2

Problem context: This regression will be used as a function approximation of (a component) the value function in a reinforcement learning algorithm. Because of that the constraints have to be enforced by the regression model, and can not be hand controlled. Moreover, the regression will be run several times with increased number of samples.

  • $\begingroup$ the requirements of boundedness makes things more complicated. Rational function fitting could solve all of your problems, but it's definitely not a robust approach: it will likely succeed on a data set and fail on another one, without warning. The second question you quote is actually quite a good suggestion: mgcv is a great package, and unless constrained polynomial regression is good enough for your application, you won't get out of this without using some specialized packages (or reinventing A LOT of wheels). Why don't you show us your data, so that we can help you better? $\endgroup$ – DeltaIV Apr 10 '17 at 14:54
  • 1
    $\begingroup$ If $x_2$ is categorical, what does it mean that $f$ is "smooth" in $x_2$? Suppose $x_2\in[1,2,3]$. Do you mean that $f(x_1,x_2)$ must be smooth for values of $x_2$ between 1 and 2, and 2 and 3, even if such values are never seen in your data? Also, are you opposed to R packages because you are not going to code the model in R? Logistic regression with monotonicity constraints is easily implemented in R with the package glmnet. You may implement it by hand in any languages, but you need at least access to some library/tools for optimization. $\endgroup$ – DeltaIV Apr 12 '17 at 13:35
  • 1
    $\begingroup$ @DeltaV you are right the definition of "smooth" is not clear. I just wanted to point out that $f(x_1,x_2=1)$ is more similar to $f(x_1,x_2=2)$ than to $f(x_1,x_2=10)$. Regarding the logistic regression, I understand that logistic regression is expecting binomial dependent variable. In my problem $y$ is continuous and bounded between $[0,1]$. $\endgroup$ – PolBM Apr 18 '17 at 10:13
  • 1
    $\begingroup$ Yes: I would accept an answer based on an R package. $x_2$ measures time (in particular hours, but only few hours each day are measured). The questions mentions that $x_2$ is categorical and ordered. Feel free to add your answer! $\endgroup$ – PolBM Apr 18 '17 at 11:39
  • 1
    $\begingroup$ just for the record. I have decided to implement a nonparametric regression based on the mgcv R package. At the end the it boils down to properly state the problem as an pcls problem. $\endgroup$ – PolBM Apr 21 '17 at 7:12

Logistic regression with box constraints served me well in the past for problems similar to yours. You are interested in prediction, not in inference, so, as long a suitable estimate of generalization error (for example, 5-fold cross-validation error) is low enough, you should be ok.

Let's consider a logit link function and a linear model:

$$\beta_0+\beta_1 x_1 +\beta_2 x_2 =\log{\frac{\mu}{1-\mu}}$$

where $\mu=\mathbb{E}[y|x_1]$


$$\frac{\partial \mu}{\partial x_1}=\beta_1 \frac{\exp{(-\beta_0-\beta_1 x_1 -\beta_2 x_2)}}{(1+\exp{(-\beta_0-\beta_1 x_1 -\beta_2 x_2)})^2}>0 \iff \beta_1>0 $$

Thus constraint 1 and 2 are satisfied if you just use logistic regression with the constraint that $\beta_1>0$. In general, monotonicity constraints with respect to one or more variables are relatively easy to enforce with GLMs (Generalized Linear Models) such as logistic regression, because the monotonicity of the link function and the fact that it's expressed it as a linear function of the predictors imply that $\mu$ is always monotonic with respect to the continuous predictors.

An R package which supports logistic regression with box constraints (constraints of the type $a_i\leq\beta_i\leq b_i$) is glmnet. Its usage is a bit different from other regression functions in R, so have a look at ?glmnet. Constraint 3 wouldn't need specific attention in most cases, because most R regression functions will automatically encode categorical variables using dummy variables. Unfortunately, glmnet is one of the few functions which doesn't do that. You need to use model.matrix to solve this: if my_data holds your observations $X=\{x_{1i},x_{2i}\}_{i=1}^N$, then

M <- model.matrix(~ x1 + x2, my_data)

will build a design matrix suitable for use with glmnet.

The only limitation of this approach lies in the fact that we have modeled the logit function as a linear function of the predictors. This may prove not flexible enough for your problem: in other words, you could get a large cross-validation error. If this is so, you should look into nonparametric logistic regression - here, however, you need to fit GAMs (Generalized Additive Models), not GLMs, and imposing monotonicity becomes more complicated. The package mgcv and the function mono.con are your friends here - you'll need to read quite a lot of documentation. Gavin Simpson's answer to question

How to smooth data and force monotonicity

which you linked in your question, has a good example.

Finally, I reiterate that this approach (as well as all other approaches which rely on logistic regression, whether Bayesian or frequentist) only makes sense because you need a quick tool to approximate in an automated way multiple unknown functions inside your reinforcement learning workflow. $y|\mathbf{x}$ doesn't really have the binomial distribution, so you cannot expect to get realistic estimates of standard errors, confidence intervals, etc. If you need a real statistical model, which would give you not only point estimates but also realistic prediction intervals, then you need to take into account the real conditional distribution of your output. This question might help:

Judging the quality of a statistical model for a percentage

  • $\begingroup$ @PolinLnd thanks, but what's the constraint I'm not answering? Constraint 3 is taken care of by model.matrix. If you add sample data to your question, I can show you a practical example. $\endgroup$ – DeltaIV Apr 21 '17 at 7:31
  • $\begingroup$ you are right i ahv deleted the comment. $\endgroup$ – PolBM Apr 21 '17 at 8:30

For constraint 3: dummy code $x_2$ so that you can model all categories within one regression function. EDIT: I can think of two ways to enforce similarity between coefficients of those dummy variables but they are rather cumbersome and I do not know available implementations. If you want to go down that road, you could of course always program it yourself:

a) You could penalize the difference between the coefficients in the objective: $\min \sum_{i}(y_{i} - \hat{y}_{i})^2 + \lambda |\beta_{2}-\beta_{3}| + \lambda |\beta_{3}-\beta_{4}| + ... $
where $\hat{y}_{i}$ is your estimate of $y_{i}$ and and $\beta_{2}$, $\beta_{3}$, ... are the coefficients of the dummies. This would require tuning of $\lambda$.

b) You could put extra constraints such as $|\beta_{j}-\beta_{j+1}| \leq c$ where $c$ is your limit on the dissimilarity. You would have to determine $c$ for that.

For constraint 2: use a logistic regression.

For constraint 1: if there is truly a positive relationship between $x_1$ and $y$, then the regression should find a positive coefficient. I would not know the benefit of forcing negative regression coefficients to become positive if the data says the opposite.

  • 1
    $\begingroup$ On constraint 1, you can estimate the model without any constraint, and if estimated model violates the constraint, then one can think about whether enforcing it would make sense. $\endgroup$ – Matthew Gunn Apr 10 '17 at 17:14
  • $\begingroup$ I have added more context to the problem since now I realize that it is relevant. Estimating the model without constraint 1, with the hope that the regression model will be monotonically increasing in $x_1$, is not an option in this case. The regression should run without any human intervention, in unseen data. $\endgroup$ – PolBM Apr 11 '17 at 8:02

The only obvious tool based on the constraints is some form of Bayesian logistic regression. The reason is that your constraints would define the prior and the likelihood. For example, by assuming $\partial{f}/\partial{x_1}$ is positive, you are assuming that there is a zero probability that the $\hat{\beta}\le{0}$, in the linear analog problem. The bounding assumes the likelihood is some sigmoid function and the easiest way to express this would be through some function $g(h(x_1,x_2))$ where $g$ is the logistic likelihood. This just requires you to solve the prior for the relationship between $x_1$ and $x_2$ through $h$, and the relationships among the partitions as a prior.

Your “smoothness” requirement can be helped through the rather ugly dropping of the “parallel lines” assumption used in Frequentist logistic regression. To see what I mean, lets drop $x_1$ and focus on the relationship between $y$ and $x_2$.

Under the parallel lines assumption $y$ is mediated through a single parameter, say $\beta$, which is constant for all values. This implies that there is a direct relationship between how people behave when inside a frozen block of water, in cold water, in 90 or 100 degree water and in boiling water. This is obviously not the case since a person frozen in a block of water would have no activity and a person in boiling water might have brief frenetic activity followed by no activity. This is not at all like the behavior at 80-100 degrees.

When you drop the parallel lines assumption you could map $\hat{\beta}_1\dots\hat{\beta}_{10}$ in such a way that it varies in a smoothly increasing or decreasing way. It would, however, create a nightmare of a prior distribution. You would have to restrict the posterior so that, for example, $\beta_{n+1}-\beta{n}\le{k}_n;k_n>0$ for increasing functions to reach “smoothness.” This may be unnecessary, but you should plan for it. I am also assuming no interaction effect between $x_1$ and $x_2$ and that you somehow can construct a clean, proper prior. You would first have to construct this as a series of probability statements, but you may be able to simplify this if you do drop the parallel lines assumption by treating $x_2$ as an ordered partition. You would need a lot of data, because the partitioning would consume a lot of explanatory power.

One other method that may preserve your smoothness requirements is to not treat $x_2$ as a $1-10$ variable, but as a ten bit variable, where each bit has its own constant that is added to or multiplied with $x_1$. For example, if $x_2=3$ then it is coded as $0000000100$ and where that is present then $c_3$ is put in the function as either $g(c_3+h(x_1))$ or $g(c_3h(x_1))$. It could even be both a multiplication and additive constant together inside $g$ but outside $h$ such as $g(c_3+k_3h(x_1))$.

The sticking point will be the probability statements to show the relationship between the variables prior to collecting the data. You could use Bayesian model selection to test competing models.

I had to use something like this in a problem because of the weird restrictions on the variables that I faced.

  • $\begingroup$ I feel like it would be better to code $x_2$ as dummy variables as you and @Timo suggest. Could you add some literature on bayesian logistic regression and/or a some tips (package etc.) in how to implement this?. $\endgroup$ – PolBM Apr 11 '17 at 7:44
  • 1
    $\begingroup$ Although the content on logistic regression in Peter Lee's book "Bayesian Statistics" is short, it is pretty good. You can look for it at ISBN-13: 978-1118332573. If you have never done anything with Bayesian methods before William Bolstad's sophomore level book ISBN-13: 978-1118091562 will let you look at them at a basic level. The most practical book, however, is neither of these. John Krushke's two books that begin with "Doing Bayesian Analysis" covers it extensively. His books are ISBN-13: 978-0124058880 and ISBN-13: 978-0123814852. $\endgroup$ – Dave Harris Apr 11 '17 at 14:24
  • 1
    $\begingroup$ Krushke's books cover the implementation of STAN, BUGS and JAGS in the R programming language. You could also implement it directly in either STAN, BUGS or JAGS. STAN uses a new, high speed, numerical integration method. BUGS and JAGS use plain-vanilla numerical integration. No matter what you use, you will probably need to write some custom software. If you get stuck somewhere, I can help over the summer months. Post something here and I can get you my contact information. There is a brief, very steep, learning curve. Then it is really easy. $\endgroup$ – Dave Harris Apr 11 '17 at 14:32

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.