Confidence interval for GLM or the maximum of a function?

Question

Imagine I have a set of (xi,yi) measures. I can show it on a scatter plot

I want to choose the value of x that maximizes y, or I could fit a function and find the values of the parameters that maximize that function.

To put it simple I decide to fit a general linear regression model (quadratic on the variable) Y=a+bx+cx^2. I can do it with R, using the glm() funtcion. And I would get the values of a, b, c their standard errors SE(a), SE(b), SE(c), R2 and more things.

And the max would be located at x=-b/2c

My question is, How can I calculate the confidence interval for that x, the predictor? Is it just the typical formula x+-Z·S/sqrt(n)? and using some result provided by R for that S? or do I need to calculate that standard using a more complex methods, such as the square root of the sum of derivatives squared...? Or it's something completly different?

How would you do it with a simulation with x ? maybe calculate "y" for different values of x with perturbation. then choose the 95% max values of y (??) and calculate the confidence interval for the associated x's ? How would you do it?

Answer 1

In the comments there are good proposals, like the traditional delta method. A newer method which might be more exact, is constructing a profile likelihood confidence interval, which today is quite easy.

I will show below with an simulated example, in R. First write the linear regression (polynomial) model as $$ Y_i = a + b x_i + c x_i^2 + \epsilon_i $$ Now, assuming $c<0$ so as to have a maximum ... (If $c>0$ the same methods will give a confidence interval for the minimum), the maximum is at $$ M=-\frac{b}{2c} $$ $M$ is our interest (or focus) parameter, so reparametrize the regression model as $$ Y_i = a + b x_i\left( 1-\frac{x_i}{2M}\right) +\epsilon_i $$ This is now no longer linear in the parameters, but we can fit it as a nonlinear regression model:

a <- 1
b <- 1
c <- -1/2
x <- seq(from=-5,  to=5,  by=1/3)

set.seed(7*11*13) # My public seed
Y <- a + b*x  +  c*x^2  +  rnorm(length(x), 0, 2)

mydata <- data.frame(Y, x)

mod0 <- nls(Y ~ a  + b*x*(1-x/(2*M)), data=mydata,
            start=list(a=0, b=0.5, M=2))  

 summary(mod0)

Formula: Y ~ a + b * x * (1 - x/(2 * M))

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   1.2778     0.6284   2.033   0.0516 .  
b   0.9671     0.1404   6.888 1.74e-07 ***
M   0.8975     0.1571   5.711 3.99e-06 ***
  ...  

Now we call the confint function, which calls the profile function , doing likelihood profiling, and constructs the confidence interval from there:

 confint(mod0, 3)
Waiting for profiling to be done...
     2.5%     97.5% 
0.6044179 1.2658407 

We can also plot the likelihood profile for M:

Now, as an exercise readers can compare this with the delta method ...