How do I find an appropriate test statistic for the derivative of a quadratic regression curve? Suppose I've managed to fit a quadratic regression curve $Y=\beta_0+\beta_1X+\beta_2X^2$ to a dataset. Given some $X=x$, I'm looking for an appropriate test statistic to check if $x$ is an extreme point of the quadratic, i.e., I want to test the null hypothesis $\frac{dY}{dX}=0$ at $X=x$, i.e., $\beta_1+2\beta_2x=0$. What is an appropriate test statistic for this?
 A: The fit produces estimates $\hat\beta = (\hat\beta_0, \hat\beta_1, \hat\beta_2)$ which, because they are functions of the random response $Y,$ are random variables.  These random variables have a covariance matrix.  Least-squares and maximum likelihood procedures routinely make an estimate of this covariance matrix available: let's call it $\hat V.$  OLS also reports "degrees of freedom" for $\hat V.$
Given that you selected $x$ independently of any reaction to the data or the regression results, you can treat it as constant.  Thus, $\hat\beta_1 + 2\hat\beta_2 x$ is a random variable whose properties are determined by the joint distribution of $(\hat\beta_1, \hat\beta_2).$  We readily compute that this is an unbiased estimator of $\beta_1 + 2x\beta_2$ because
$$\mathbb{E}(\hat\beta_1 + 2\hat\beta_2 x) = \mathbb{E}(\hat\beta_1) + 2x\mathbb{E}(\hat\beta_2) = \beta_1 + 2x\beta_2.$$
Its variance equals
$$\operatorname{Var}(\hat\beta_1 + 2\hat\beta_2 x) = \operatorname{Var}(\hat\beta_1) + (2x)^2 \operatorname{Var}(\hat\beta_2) + 2x \operatorname{Cov}(\hat\beta_1, \hat\beta_2).$$
The three variance-covariance terms are estimated from the corresponding entries in $\hat V;$ namely, $\operatorname{Var}(\hat\beta_1) \approx \hat V_{\hat\beta_1,\hat\beta_2}$ and so on.
Thus you can test hypotheses about $\beta_1 + 2x\beta_2$ (or indeed about any constant linear combination of the parameters) exactly as you would test hypotheses about the parameters individually; namely, compare the difference between the estimate and its hypothesized value (under the null) to its standard deviation:


*

*For least squares procedures, the test statistic for a comparison to $0$ is $$t = \frac{\hat\beta_1 + 2x\hat\beta_2 - 0}{\sqrt{\hat V_{\hat\beta_1,\hat\beta_1} + (2x)^2 \hat V_{\hat\beta_2,\hat\beta_2} + 2x \hat V_{\hat\beta_1,\hat\beta_2}}}$$ and it ought to approximately have a Student t distribution with the same degrees of freedom associated with either of $\hat\beta_1$ or $\hat\beta_2$ (which is routinely reported in software output).

*For maximum likelihood procedures, the same statistic is suggestively called $Z$ and it ought to approximately have a standard Normal distribution.
Compute the p-value accordingly.

To illustrate this procedure (and demonstrate it works), the R code below simulates from this model, carries out OLS fit and the foregoing calculations for a range of possible $x$ (although in the application it is valid to use these p-values to make decisions only about a single $x$), and plots the results.

At the left the fit is the dark curve, in contrast to the true underlying model shown in gray.  At the right a p-value of $0.05$ is plotted as a horizontal dashed line for reference.  Note the zoomed horizontal ($x$) scale at the right.
Observe, too, how $\beta_1 + 2x\beta_2 = 0\beta_0 + \beta_1 + 2x\beta_2$ and its estimate are represented in this vector-oriented code as the array of coefficients $(0, 1, 2x)$ (stored in the contrast variable). This is possible because the software always organizes the rows and columns $\hat V$ in this order.
n <- 10              # Sample size
beta <- c(5, 0, -1)  # True parameters
sigma <- 1           # Error SD
#
# Generate a dataset.
#
set.seed(17)
x <- seq(-2, 2, length.out=n)
y.0 <- cbind(1, x, x^2) %*% beta
y <- y.0 + rnorm(n, 0, sigma)
#
# OLS regression.
#
fit <- lm(y ~ x + I(x^2))
#
# Extract estimates.
#
beta.hat <- coef(fit)
V.hat <- vcov(fit)
df <- fit$df.residual
#
# Compute p-values for a range of possible x.0.
#
x.0 <- seq(-0.5, 0.5, by=0.005)
p.values <- sapply(x.0, function(x.0) {
  contrast <- c(0, 1, 2*x.0)
  se <- sqrt(contrast %*% V.hat %*% contrast)
  t.stat <- contrast %*% beta.hat / se
  2*pt(-abs(t.stat), df)
})
#
# Create a data structure for plotting the model and the fit.
#
X <- data.frame(x = seq(min(x), max(x), length.out=101))
X$y <- predict(fit, newdata=X)
X$y.0 <- with(X, cbind(1, x, x^2) %*% beta)
#
# Plot data and the results.
#
par(mfrow=c(1,2))
plot(x, y, main="Data, Model, and Fit")
with(X, lines(x, y, lwd=2))
with(X, lines(x, y.0, col="Gray", lwd=2))

plot(x.0, p.values, type="l", log="y", main="P-values", xlab="x", lwd=2)
abline(h=0.05, col="Red", lty=3)
par(mfrow=c(1,1))

