I have 30 pairs of observations $(x_i, y_i)$ and hypothesize that in the observed range of $x$ a non-monotone relationship $y_i = f(x_i) + \epsilon$ may hold. I am willing to assume that $\text{Var}(\epsilon | x_i) =\sigma^2$.

I am interested in

  1. Testing if $f(x_i)$ is non-monotone, and if I reject monotonicity then
  2. Estimating $ x^* = \text{argmin} f(x)$.

Approach so far

My approach has been to use linear regression with a quadratic basis expansion of $x$ thus $$ E(Y|x) = a + bx + cx^2, $$ and then conduct inference on $\hat c$ and $\hat x^* = -\hat{b}/(2\hat{c})$ . If $\hat c>0$, and $\hat x^* = -\hat{b}/(2\hat{c})$ lies in the range of the observed data with 95% confidence, I conclude that $f$ is concave and that a sign change in the derivative of $f$ occurred. For $\hat x^*$ I am using the bootstrap, since the delta-method doesn't seem reliable for a ratio of estimators. In the figure below, see an example of such a quadratic fit and the 95% confidence interval on $x^*$, a one-sided CI on $c$, and a 95% CI on $b$, which gives the slope at $x=0$ if the quadratic model holds.

enter image description here


  1. I am looking for references to a principled treatment of this problem (maybe in response surface modeling?) I am nervous about inventing a wheel here.

  2. Although I doubt I can do more than fit a quadratic since I only have 30 data points, there is something alarming about assuming a polynomial holds. I am vaguely aware of non-parametric convex regression, could it be used somehow here?

  • $\begingroup$ I would estimate the first derivative of f(assuming its regular enough). And you can use any base to do your estimation (here you use polynomial, it could be any functional base, as long as you have few parameter to estimate this is not a problem). I don't have time to dig further, but i can say one last important thing. Its always better to directly estimate what we want to estimate and not do multiple estimation step, its less precise and hide many problems. $\endgroup$
    – PauZen
    Dec 17, 2019 at 9:50
  • $\begingroup$ From the plot, it appears that using values where AUC34 = 0 is not contributing to the overall model. If you separately model AUC34 equal to zero and AUC34 greater than zero, combining those separate models, would this be useful? $\endgroup$ Dec 17, 2019 at 13:29
  • $\begingroup$ @JamesPhillips In my view, the AUC34 = 0 are very important contributors, and eliminating them makes it much harder to estimate the functional relationship, since it dramatically reduces the variation in $x$. You might worry that there could a threshold effect, as do I. However, we address that in a separate analysis. $\endgroup$
    – Andrew M
    Dec 17, 2019 at 20:38
  • $\begingroup$ @PauZen: Thanks for the comment. Under the quadratic model $E (f'(x)) = b + 2cx$ and indeed that is the quantity I am using for inference. However, I don't see how I can avoid an intersection test, and inference on a non-linear function of $(b, c)$, as I need to verify that both $\exists x \in \text{range}(X)_i: f'(x) = 0$ and $c>0$, since just finding a zero doesn't insure I found a local minimum, nor does the fact that $c>0$ mean that a minimum occurs. So an ideal model would provide inference on the the roots of $f'(x)$ for $f$ unimodal. $\endgroup$
    – Andrew M
    Dec 17, 2019 at 20:49
  • $\begingroup$ By your model (the quadratic fit) you assume by assumption a lot about the form of your function (especially about global minimum). $\endgroup$
    – PauZen
    Dec 18, 2019 at 16:03

1 Answer 1


Unimodal regression generically estimates $f(x)$

Isotonic regression provides non-parametric estimates of a function $f(x)$ such that $f(x_i) \leq f(x_{i'})$ for $x_i < x_{i'}$ (with the obvious changes in the inequalities for a decreasing function). The solution can be shown to be a step function, with knots at a subset of the observed values of $x$.

Frisen (1986) describes a modification to isotonic regression that permits estimation of unimodal functions, that essentially boils down to searching points $x$ as candidate turning points, and then picking the turning point that minimizes the residual sum of squares.

This partially addresses my question, in that the isotonic model is nested within the unimodal model, so there should be a way to compare model classes. However I remain unsure about how I can conduct inference on the "turning point" of the unimodal regression. I will have to experiment.


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.