# Confidence intervals around functions of estimated parameters

Suppose I estimate the following model, $$y = \hat{\beta_0} + \hat{\beta_1} X_1 + \hat{\beta_2} X_2\ ,$$ and am additionally interested in the quantity $$\gamma = \sqrt{\beta_0^2 - 4\beta_1\beta_2}\ .$$ An estimate of the above, $$\hat{\gamma}$$, is clearly given by inserting the estimated $$\hat\beta_i$$'s in the above formula.

But how does one go about computing the confidence interval around such a "compound" quantity that is not directly estimated?

My attempt.

My "empirical" approach is as follows. I generate a multivariate normal distribution of size $$n$$ using the estimated means $$\hat\beta_i$$'s and the covariance matrix of the fit. I consequently have a distribution of size $$n$$ for the $$\gamma$$. The mean of this distribution $$\bar\gamma$$ is naturally then close to the $$\hat\gamma$$ computed above, and I now additionally have a standard deviation of this distribution $$\sigma$$.

Using $$\sigma$$ to generate a confidence interval.

If I use the standard definition $$\bar\gamma \pm t(\alpha, n-1) \sigma / \sqrt{n}$$, where $$t$$ is the PPF of the $$t$$-distribution, I run into the following quandary: since I pick $$n$$ myself, I can in fact make the standard error as small as I like. This is not what I want: I am looking for a "natural" confidence interval which depends solely on the uncertainty in my estimates of the $$\beta$$'s. Therefore, my idea is the following: $$\sigma$$ stays positive, and seems to "converge" to its true value for large $$n$$: hence, 95% of the values in my $$\gamma$$ distribution lie within the interval $$\bar\gamma \pm 2\sigma$$.

Is this approach valid?

• Is your model a quadratic function? $$y = \hat{\beta_0} + \hat{\beta_1} x + \hat{\beta_2} x^2$$ What for do you need to estimate $\gamma$? Commented Nov 9, 2021 at 8:50
• @SextusEmpiricus No, my model isn't a quadratic function. I attempted to frame the question as a "minimal working example"--my particular use case is in fact slightly more complicated, but I didn't want the discussion to get bogged down in details. Commented Nov 9, 2021 at 12:46

Two common approaches for this problem are to calculate the non-linear combination of the coefficients directly from the regression or to bootstrap it.

The variance in the former is based on the "delta method", an approximation appropriate in large samples. This was suggested in the other answer, but statistics software can make the calculation a whole lot easier.

The variance for the latter comes from resampling the data in memory with replacement, fitting the model, calculating the coefficient combination, and then using the sampled distribution to get the confidence interval.

Here's an example of both using Stata:

. sysuse auto, clear
(1978 automobile data)

. set seed 11082021

. regress price mpg foreign

Source |       SS           df       MS      Number of obs   =        74
-------------+----------------------------------   F(2, 71)        =     14.07
Model |   180261702         2  90130850.8   Prob > F        =    0.0000
Residual |   454803695        71  6405685.84   R-squared       =    0.2838
-------------+----------------------------------   Adj R-squared   =    0.2637
Total |   635065396        73  8699525.97   Root MSE        =    2530.9

------------------------------------------------------------------------------
price | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
mpg |  -294.1955   55.69172    -5.28   0.000    -405.2417   -183.1494
foreign |   1767.292    700.158     2.52   0.014     371.2169    3163.368
_cons |   11905.42   1158.634    10.28   0.000     9595.164    14215.67
------------------------------------------------------------------------------

. nlcom (gamma_dm:sqrt(_b[_cons] - 4*_b[mpg]*_b[foreign]))

gamma_dm: sqrt(_b[_cons] - 4*_b[mpg]*_b[foreign])

------------------------------------------------------------------------------
price | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
gamma_dm |   1446.245   361.0078     4.01   0.000     738.6823    2153.807
------------------------------------------------------------------------------


The 95% CI using the delta method is [738.6823, 2153.807].

Boostrapping yields [740.5149, 2151.974], which is fairly similar:

. bootstrap (gamma_bs:sqrt(_b[_cons] - 4*_b[mpg]*_b[foreign])), reps(500) nodots: regress price mpg foreign

Linear regression                                          Number of obs =  74
Replications  = 499

Command: regress price mpg foreign
[gamma_bs]_bs_1: sqrt(_b[_cons] - 4*_b[mpg]*_b[foreign])

------------------------------------------------------------------------------
|   Observed   Bootstrap                         Normal-based
| coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
gamma_bs     |
_bs_1 |   1446.245   360.0728     4.02   0.000     740.5149    2151.974
------------------------------------------------------------------------------
Note: One or more parameters could not be estimated in 1 bootstrap replicate;
standard-error estimates include only complete replications.


Your proposed solution would work if you have lots of data, but here it does not do so well with only 74 observations:

. quietly regress price mpg foreign

. corr2data b_mpg b_foreign b_cons, n(500) means(e(b)) cov(e(V)) clear
(obs 500)

. gen gamma_sim = sqrt(b_cons - 4*b_mpg*b_foreign)
(3 missing values generated)

. sum gamma_sim

Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
gamma_sim |        497    1426.183    366.6408   197.5594   2397.263

. display "[" %-9.4f r(mean) + invttail(r(N)-1,.975)*r(sd) ", " r(mean) + invttail(r(N)-1,.025)*r(sd) "]"
[705.8224 , 2146.5434]


The CI here is [705.8224 , 2146.5434], which is noticeably different from the two CIs above.

My thought is that if you are going to simulate, you might as well bootstrap and not rely on the normal approximation that is only valid asymptotically. If you have lots of data, the difference between bootstrapping and sampling from MVN parameterized by estimated coefficients and variance should not be noticeable.

• You mention 'with replacement', so is that bootstrapping resampling by using the residuals and shuffling them? This would not be necessary if we know that the distribution of the coefficients is multivariate normal (which is approximately true since it is a sum of many variables) and we can sample from that multivariate normal. With small samples this would give a better outcome. Commented Nov 9, 2021 at 6:36
• @Anthony the concept is the same, compute the distribution of the non-linear combination $\sqrt{b^2-4ac}$, by simulation rather than an analytic computation. The only difference is that you make use of the assumption that the sample distribution for the coefficients follows a particular distribution, a multivariate normal distribution. So you use that distribution for the simulation. With bootstrapping you use an empirical distribution (you use the observed residuals). Bootstrapping is preferred when the assumption is not accurate (low sample size and non-normal error distribution). Commented Nov 9, 2021 at 8:45
• I'm marking this as the accepted answer since it is the most general. Thanks for all the responses and discussion. Commented Nov 9, 2021 at 8:49
• @Anthony I added some thoughts about your MVN sampling approach at the bottom of my answer. Commented Nov 9, 2021 at 23:32
• I'm actually surprised the third method doesn't align better, given that it sits right inbetween delta and bootstrap in terms of approximations (MVN approximation of original parameters, but exact transformation of MVN approximation). I noticed that you have 3 "missing" values in the proposed method; I suspect these are actually square roots of negative values? I'm curious if you used the empirical CIs (i.e. 2.5th, 97.5th percentile of simulated sample) if the values would align better? Commented Nov 10, 2021 at 4:30

Usually we take normality assumption for linear regression models. That is, $$y_i\sim N(\beta^Tx_i,\sigma^2)$$. From this assumption we derive the asymptotic distribution of $$\hat{\beta}$$, which is also normal: $$\hat{\beta}\overset{.}{\sim}N(\beta,\sigma^2(X^TX)^{-1})$$. Let us denote the covariance matrix of $$\hat{\beta}$$ as $$Var(\hat{\beta})=V=\sigma^2(X^TX)^{-1}$$.

Next, we take $$\gamma=g(\beta)=\sqrt{\beta_0^2-4\beta_1\beta_2}$$ so $$\hat{\gamma}=g(\hat{\beta})$$. The function $$g$$ is differentiable and let us assume it is non-zero. Equipped with these, we can apply the delta method to get the following estimator:

$$\hat{\gamma}\sim N\left(\gamma,\nabla g(\hat{\beta})^T\cdot V\cdot \nabla g(\hat{\beta})\right)$$

Where $$\nabla g(\hat{\beta})$$ is the gradient vector of $$g$$ :

$$\nabla g(\beta)=\begin{pmatrix} \frac{\beta_0}{\sqrt{\beta_0^2-4\beta_1\beta_2}} \\ \frac{-2\beta_2}{\sqrt{\beta_0^2-4\beta_1\beta_2}} \\ \frac{-2\beta_1}{\sqrt{\beta_0^2-4\beta_1\beta_2}} \end{pmatrix}$$

Now, you can easily write the value of $$\nabla g$$ at $$\hat{\beta}$$ and then substitute $$\hat{\beta}=(X^TX)^{-1}X^Ty$$ and obtain a nicer form for the variance of $$\hat{\gamma}$$ (this does require some linear algebra work). Finally, you can write down your CI.

• Thanks for the excellent answer! I happened to stumble across this webpage after I posted my question: astro.rug.nl/software/kapteyn/… Could you confirm if equation (100) is the same as would arise from the delta method above? (It sure looks like it to me.) Commented Nov 8, 2021 at 13:26
• (+1) But, as you stated, it is the asymptotic distribution which is normal distributed, and in addition, this Delta method is a 1st order approximation. So, this is only an estimate of the actual distribution when the sample size is finite. For small sample sizes the error might be large and an estimate with the Monte Carlo approach will probably perform better. Commented Nov 8, 2021 at 15:03
• @SextusEmpiricus By a Monte Carlo simulation, do you mean something similar to the approach I described? Would you care to comment? Commented Nov 8, 2021 at 15:35
• @SextusEmpiricus While I do not disagree over the efficiency of the MC, it is a nonparametric method. Why not using parametric methods to the farthest extent possible (under the problem setting)? Also consider (a) the difficulty of sampling WRT $Cov(X)$; (b) generating enough samples so $\hat{\beta}$ converges; (c) given previous two conditions, the difficulty of generating different resampled datasets. For all of the above, $X$ should be large enough - but then again, if it's large enough then we should have no problem using the asymptotic distributions. Commented Nov 8, 2021 at 15:46
• @Spätzle I consider the problem more as computing the distribution of $\sqrt{\beta_0^2 -4\beta_1\beta_2}$ where $\hat\beta_0,\hat\beta_1,\hat\beta_2 \sim MVN(\hat{\mu},\hat\Sigma)$ and $\hat{\mu},\hat\Sigma$ is given (based on estimates from the data). But yes the estimation of the $\hat{\mu},\hat\Sigma$ creates an additional source of error... Commented Nov 8, 2021 at 16:03

It seems like you are estimating the discriminant of a quadratic function, ie. your function is

$$y = \hat{\beta_0} + \hat{\beta_1} X_1 + \hat{\beta_2} X_2 = \hat{\beta_0} + \hat{\beta_1} X + \hat{\beta_2} X^2$$

We can rewrite this in terms of the roots

$$y = \hat{\beta_2} (X-r_1)(X-r_2)$$

And now we have $$\gamma = \beta_2 (r_2-r_1)$$

If you estimate $$\beta_2$$, $$r_1$$, and $$r_2$$ with a non-linear estimation method and obtain an estimate for the approximate multivariate normal distribution of the estimates (based on an estimate of the Fisher information matrix), then you can describe an estimate distribution for $$\gamma$$ as a product of two correlated non central normal distributions.

• how are coefficient variance estimates found in nonlinear estimation? Commented Nov 14, 2021 at 17:16

There's one big problem with your proposed method, although it can be fairly easily remedied. In particular, you compute the estimate of variance of $$\gamma$$ as $$\frac{V(\gamma_i)}{n}$$ (where $$V(\gamma_i)$$ is variance of your samples generated by first sampling from the multi variance normal and then pushing through your function), when the standard error should really just be $$V(\gamma_i)$$.

But furthermore, I would suggest you use the simulated percentiles to create the CI rather than using the mean + se formula.

And one final note, we note that the sqrt function is not defined for negative values. This means if a random draw gives you $$4\beta_1 \beta_2 > \beta_0^2$$, your estimate is not real valued. Assuming that this is not what you want to do, I would say you should consider simply dropping the simulated values when this is true. That may seem adhoc, but it perfectly aligns with defining a restricted parameter space where $$4\beta_1 \beta_2 \le \beta_0^2$$.