Estimating 95% confidence interval for vertex of quadratic fit

Question

Can anyone tell me how to obtain the 95% confidence intervals for the x and y coordinates of the vertex (peak in my case) of a quadratic fit? Is there a way to get these CIs in R? Thanks.

Answer 1

There are several ways to approach this problem, but the most straightforward I can think of is this one:

If the population equation is $y = \alpha x^2 + \beta x + \gamma$

then the vertex is at $x_v = -\frac{\beta}{2\alpha}$ (as long as $\alpha\neq 0$)

Equivalently, $2\alpha x_v-\beta=0$.

Correspondingly, the fitted equation is $\hat y = a x^2 + b x + c$

and the estimated vertex is at $\hat{x}_v = -\frac{b}{2a}$ (if $a\neq 0$).

Consequently, a set of values that form a $1-\alpha$ confidence interval (perhaps better called a consonance interval in this case) for $x_v$ are the values, $d$ which would result in acceptance of:

$H_0: 2\alpha d-\beta=0$

at significance level $\alpha$.

Such a set of values might be constructed relatively readily. (Indeed one needn't even be very smart about it, since without doing any algebra we can work our way out from $\hat{x}_v$ until the decision on the hypothesis test changes to rejection, and then perform a binary search, though this assumes the set of values resulting in acceptance form a single interval.)

However, it looks to me like it's reasonably analytically straightforward, since we can, for any given $d$, write a t-interval and at a fixed significance level, for $2\alpha d-\beta$; that will involve $d$ in both the numerator and denominator of the t-statistic, simple manipulations yield a quadratic inequality in $d$, $Q(d)<0$.

The equality is easily solved (involving coefficients in terms of $a$, $b$ and their variances and covariances), and it's then a matter of checking that $Q(\hat x_v)$ is negative and lies between those points of equality, whence you have a simple interval of values of $x_v=d$ that would not lead to rejection and which contains the naive point estimate $Q(\hat x_v)$.

(criticism of this idea is encouraged; I may well be missing some obvious issue)

---

A second approach would be to simulate or bootstrap the distribution of $\hat x_v$.

The bootstrap distribution of $\hat x_v$ could be obtained by resampling (with replacement) $n$ rows of the augmented data matrix $A=[X|y]$ and computing the pseudo-sample point estimate $\hat x_v^*$ each time; there are various techniques for generating an interval from that, including the relatively simple 'take the middle $1-\alpha$ of the resampling distribution of $\hat x_v^*$'.

Another form of bootstrap would be a parametric bootstrap (basically a kind of simulation) - one way to do that: if we condition on the observed $x$, we could simulate from the distribution of $y-\hat y$ to generate new samples, $y^*$ and proceed in similar manner to the nonparametric bootstrap above (but here we're simulating only the noise and condition on $x$, rather than the bivariate $(X,y)$.

Simulation could be performed in a number of other ways. We could for example, simulate the distribution of $-\frac{\hat\beta}{2\hat\alpha}=-\frac{b}{2a}$ under normal assumptions and the same pattern of $x$s, though I believe the shape of that could depend on $\sigma$ (its best if some pivotal quantity is used as the basis for the simulation).

---

Given we're trying to find an interval for a particular $x$, another way to look at the problem is as a kind of inverse regression problem, but there are some details to sort out here.

---

A reference:

Florenz Plassmann and Neha Khanna (2007)
Assessing the Precision of Turning Point Estimates in Polynomial Regression Functions
Econometric Reviews, vol. 26, issue 5, pages 503-528

The 2003 version of paper here

discusses 3 approaches:

a) Taylor expansion/Delta method

b) calculation based on assuming that the coefficient estimates are bivariate normal

c) Finite sample estimates based on MCMC

Answer 2

There are a few different approaches that would be fairly straight forward.

One option would be to fit using non-linear least squares where the vertex is one of the parameters (e.g. y = a + b*(x-c)^2 ). Then you can use profiling to calculate a confidence interval on the parameter of interest.

Another option is to do a Bayesian fit using McMC methods and estimate the x and y coordinates of the vertex from each iteration and use those values to estimate the posterior and credible interval.

Or, if the residuals are assumed to be normal (or at least close enough) then the estimated coefficients from the linear model follow a multivariate normal distribution. You can generate many observations from this distribution, compute the x and y coordinates of the vertex from these values and use those to estimate the distribution of the vertex points and compute the confidence interval from that.