Confidence interval for vertex of quadratic model, transforming parameters estimated by MLE

Question

If I have a model $y=ax^2+bx+c+\epsilon$, and I use maximum likelihood estimation (in R, with nlm function) to estimate $(a,b,c)$ with a Hessian matrix $H$ as a results, can I use this to calculate the (95%) C.I. for the $x$-coordinate of the vertex (i.e. turning point; the minimum or maximum depending on the sign of $a$), $\mu=-\frac{b}{2a}$?

Answer 1

Small Disclaimer

You do not need the nlm package to solve your problem. You have a linear model, it is still linear although your variables are non-linear.

Small Intro

My answer is based on the comment I made to your question. I am going to use the Delta method to construct a confidence interval for you.

The Delta method is based on an approximation of the variance of your new parameter $\mu$. So this is an approximation, and we will verify (somewhat) how good it is with a simulation example.

This is in fact not an easy problem. The parameters in a linear model (your model is linear in the parameters, thus linear) follow a multivariate normal distribution and they are not independent. Thus the exact distribution is not easy to derive, and I do not have an idea of how to find it now. So let the approximation do.

Let's solve this

I will call your new parameter $\mu:= g(a,b)$.

The function you provide is $g(a,b)=-\frac{b}{2a}$. It is worth mentioning that a function of MLEs is an MLE, so the expected value is $$ \text{E}[-\frac{B}{2A}]=-\frac{b}{2a} $$ Where $A$ and $B$ are the RVs corresponding to the estimators for $a$ and $b$ (This may be some abuse of notation). Now the confidence interval is symmetric, so the center will be the expected values.

On the top of page 6 in the document in Delta method link, your problem is solved, i.e. when the function $g$ is the quotient of two RVs. The only difference is that you scale the quotient by $\frac{-1}{2}$. So the variance of $g(A,B)$ is: $$ \begin{align} \text{Var}[g(A,B)] &= \text{Var}[-\frac{B}{2A}]\\ &= \frac{1}{4}\text{Var}[\frac{B}{A}]\\ &\approx \frac{1}{4} (\frac{b}{a})^2(\frac{\text{Var}[B]}{b^2} + \frac{\text{Var}[A]}{a^2} - 2\frac{\text{Cov}[A,B]}{ab}) \end{align} $$ You should have everything you need to calculate this, just plug in the point estimates of $a$ and $b$ and find the corresponding variances and covariances in the covariance matrix of the parameters.

Your 95% confidence interval (I picked $\alpha$ as 5% for the sake of demonstration) thus becomes $$ -\frac{\hat{b}}{2\hat{a}}\pm 1.96\cdot\sqrt{\text{Var}[-\frac{\hat{b}}{2\hat{a}}]} $$

Small Simulation Study

Now let's check how good this estimator really is. This is quite a lot of code, try it out to get a feel for what is happening:

# Start by generating data

a <- 1.5
b <- 2
c <- 3
trueVal <- -0.5*b/a

x <- seq(-5,5,length=100)
y <- a*x^2 + b*x+c

plot(x, y, main = "Some quadratic function")

# Let's remove some of the x's
x <- x[-sample(1:100,80)]
x2 <- x^2
y <- a*x2 + b*x+c

plot(x, y, main = "Removed points")

# Add some error to this
yErr <- y + rnorm(20, mean=0, sd=5)
plot(x, yErr, main = "With error")

# Now let's estimate the parameters
fm <- lm(yErr~x2+x)
params <- coef(fm)

# Look at parameters (They are very close to our defined a,b and c)
summary(fm)

# We can find the covariance matrix by looking at
varMat <- vcov(fm)

# So we can get an estimate for the 95% confidence interval as:
leftVal <- -0.5*params[3]/params[2]-
  1.96*sqrt(0.25*(params[3]/params[2])^2*
    (varMat[3,3]/params[3]^2 +
     varMat[2,2]/params[2]^2 -
     2*varMat[2,3]/(params[2]*params[3]) ))

rightVal <- -0.5*params[3]/params[2]+
  1.96*sqrt(0.25*(params[3]/params[2])^2*
    (varMat[3,3]/params[3]^2 +
     varMat[2,2]/params[2]^2 -
     2*varMat[2,3]/(params[2]*params[3]) ))

# Let's encapsulate this in a function
confQuot <- function(params,varMat){
  leftVal <- -0.5*params[3]/params[2]-
    1.96*sqrt(0.25*(params[3]/params[2])^2*
      (varMat[3,3]/params[3]^2 +
       varMat[2,2]/params[2]^2 -
       2*varMat[2,3]/(params[2]*params[3]) ))
  rightVal <- -0.5*params[3]/params[2]+
    1.96*sqrt(0.25*(params[3]/params[2])^2*
      (varMat[3,3]/params[3]^2 +
       varMat[2,2]/params[2]^2 -
       2*varMat[2,3]/(params[2]*params[3]) ))
  return(c(leftVal,rightVal))
}

# Now let's simulate this and see how often the 
# true value falls in this interval
numSamps <- 10000
counter <- 0
for(i in 1:numSamps){
  yErr <- y + rnorm(20, mean=0, sd=5)
  fm <- lm(yErr~x2+x)
  params <- coef(fm)
  varMat <- vcov(fm)
  confs <- confQuot(params,varMat)
  if(trueVal>confs[1] & trueVal<confs[2]){
    counter <- counter + 1
  }
}

# Simulations give that the true value falls in 
# the interval this proportion of the simulations
(counter/numSamps)

I get the following value from my runs:

> (counter/numSamps)
[1] 0.9352

So it is not exactly 95% but very close indeed. If you are reporting this for something I would at least perform some kind of a simulation like this to get a feel for the true $\alpha$ in your interval. You would also want to try the simulation with different errors and different number of points.

Hope this helped.