Quadratic polynomials in glm and confidence interval of extremum

It is frequent that a variable is modeled using quadratic equation in glm if we suspect that there is a minimum or maximum. However it is less easy to characterize where is located the extremum. One solution is given here but it does not allow to direct fit the data:
Stimson JA, Carmines EG, Zeller RA (1978) Interpreting polynomial regression. Sociological Methods & Research 6 (4):515-615

Let do an example:

First generate some data:

xobs <- 1:100

# If a2 positive, the parabole has a minimum

# If a2 negative, the parabole has a maximum

# The value of a2 is related to the curvature at the point x=b

a2 <- 0.1

# The b value is the position of the maximum or minimum

b <- 40

# The c value is a shift parameter. The "0" will be at c !

c <- 100

yobs <- a2*(xobs-b)^2+c

plot(xobs, yobs, ylim=c(0, 400), type="l", bty="n", las=1)

This second order equation is not of the classical form a2*x^2+a1*x+a0 and thus it cannot be fitted directly using glm.

It must be converted in a second step:

# It can be changed into a classical second order equation using:

a1 <- -2*a2*b

a0 <- a2*b^2+c

yobs <- a2*xobs^2+a1*xobs+a0

points(xobs, yobs, col="red", pch=3)

Now let do some optimization using glm:

yobs <- yobs + rnorm(100)*20

plot(xobs, yobs, ylim=c(0, 500), type="p", bty="n", las=1)

g0 <- glm(formula = yobs ~ xobs + I(xobs^2))

Note the use of I(xobs^2) and not poly(xobs, 2). poly() uses orthogonal polynomials by default and they cannot be used here, so do no forget to use raw=TRUE.

c <- coef(g0)

names(c) <- c("a0", "a1", "a2")

v <- vcov(g0)

rownames(v) <- colnames(v) <- c("a0", "a1", "a2")

a0; a1; a2

(b <- -a1/(2*a2))

(c <- a0 - a2*b^2)

And the confidence interval is...

# How to get the SE and CI of b and c

library(car)

# The confidence interval of the position of the extremum (b)

deltaMethod(c, "-a1/(2*a2)", v)

# The confidence interval of shift (c) 

deltaMethod(c, "a0 - a2*(-a1/(2*a2))^2", v)

# The confidence interval of the curvature at the extremum (a2)

deltaMethod(c, "a2", v)