regression with constraints

Question

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.

Answer 1

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]$

Then

$$\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

Answer 2

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.

Answer 3

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.