# Confidence box for coefficients of linear regression?

I am learning linear regression and I am trying to create a visualisation. Say I want to estimate a power model $$y=ax^b$$ using linear regression. I take the logarithm to get

$$\ln(y)=\ln(a)+b\ln(x)$$ Now I can run linear regression on log log of the data. Say I run the regression and I get that $$\hat{y}=\hat{\beta_0}+\hat{\beta_1}x$$. I understand how to compute the confidence intervals for each $$\beta_0$$ and $$\beta_1$$ individually from well known formulas. However, I would like to figure out the confidence "box" for $$(\beta_0,\beta_1)$$. I assume this box is a little smaller than the cross products of individual confidence boxes for $$\beta_0,\beta_1$$ in general, but it may not be. To my understanding, my question is related to the confidence band. The confidence band gives us a "region" within which with 95% chance our true line lies. The prediction band is: $$\hat{y}=\hat{\beta_0}+\hat{\beta_1}x \pm \sqrt{\text{Quadratic stuff}}$$ This is fine, but how do I now find the confidence bands for $$ax^b$$? My naive approach is to take the prediction band above, replace all $$x$$ with $$\ln(x)$$ and raise both sides to the power $$e$$. This of course will no longer be power functions, but will certainly resemble them. Are those the new prediction bands for the power model?

• To clarify: are you asking for the confidence interval for $\bar y$ are some given $x$? Or for the prediction interval for $y$ given some $x$? Commented Feb 29 at 0:02
• @AlexJ I would like the confidence "box" (2 dimensional interval), s.t, with probability of at least 95%, the true parameters of the regression lay in that box. I think that would correspond to confidence bands. Commented Feb 29 at 0:06
• Edited, my apologies. I think it is indeed confidence bands. Commented Feb 29 at 0:12
• This is only half a solution, because I don't have the knowledge off the top of my head. But $(\hat\beta_0, \hat\beta_1)$ should be (approximately) distributed as a 2-dimensional Gaussian with expectation $(\beta_0, \beta_1)$ and variance-covariance matrix $\Sigma$ (which you will get from your model fitting). You should be able to plot this in 2D. That's half of it. The half that's missing is getting the 95% confidence region (it will not be a square box, it will be some vaguely elliptical shape), which idk off the top of my head, sorry Commented Feb 29 at 0:14
• I still find your question a bit unclear, btw. For example, the question appears to be about a 2D confidence region for $(\beta_0, \beta_1)$. But then you're asking about "but how do I now find the confidence bands for $ax^b$" which is related but different Commented Feb 29 at 0:16

As I see it, there are two options. The first is to use the log-transformed formulation you are using. The second is to use a generalised linear model formulation. I think the second way is better, but I will illustrata both.

In both, we start with the power equation $$y=ax^b$$ that we're interesting in estimating for a dataset of $$x$$ and $$y$$ variables.

## Option 1 - log transform

This is the approach you propose in your question body.

Taking logarithms on both sides gives $$\log y= \log a + b \log x$$.

By using ordinary least squares, we are actually modelling the mean of the LOG: $$E[\log y] = \log a + b \log x$$. The parameter estimate $$\hat\beta_0$$ is estimating $$\log a$$, and $$\hat\beta_1$$ is estimating $$b$$. Note $$(\hat\beta_0, \hat\beta_1)$$ is a 2D Gaussian with variance-covariance matrix estimated during the modelling process.

If you want to get a confidence interval on $$\log y$$ (e.g. for plotting on a line chart), that's easy. For a sequence of $$x$$ values, use some standard formulae for variances:

\begin{aligned} \text{var}(\log y_i) &= \text{var}(\hat\beta_0 + \hat\beta_1 \log x_i) \\ &= \text{var}(\hat\beta_0) + (\log x_i)^2 \text{var}(\hat\beta_1) + (2\times 1\times \log x_i)\text{cov}(\hat\beta_0, \hat\beta_1) \end{aligned}

You can get all of these values from the output of your model.

BUT, there is an issue here. Exponentiating $$E[\log y]$$, i.e. $$\exp(E[\log y])$$, doesn't return $$y$$. Which brings me to the second option...

## Option 2 - log-link GLM

Again we start with $$y=ax^b$$. Now we'll model instead $$\log(E[y])=\gamma_0 + \gamma_1 \log x$$. This can be fit with most statistical packages using a GLM with a log link. Again, $$(\hat\gamma_0, \hat\gamma_1)$$ is a (approximately) 2D Gaussian with variance-covariance matrix estimated during the modelling process.

Here's the benefit of this way. When we go to estimate the variance

\begin{aligned} \text{var}(\log(E[y_i])) &= \text{var}(\hat\gamma_0 + \hat\gamma_1 \log x_i) \\ &= \text{var}(\hat\gamma_0) + (\log x_i)^2 \text{var}(\hat\gamma_1) + (2\times 1\times \log x_i)\text{cov}(\hat\gamma_0, \hat\gamma_1) \end{aligned}

Now if we exponentiate both sides it's more useful.

• $$\exp(\log(E[y]))=E[y]$$, which is what we want
• We can calculate a confidence interval on $$\log(E[y])$$ in log-space, using the above equation. Let's call that interval $$[L,U]$$. Exponentiating that interval gives us an (asymmetric) confidence interval on $$E[y]$$

Here's some R code with an example

library(ggplot2)

f <- function(x, a, b) a * x^b
A <- 0.5 # set some parameter values
B <- 1.9

set.seed(12345)

df <- data.frame(
x = 1:50
)
df$$y <- f(df$$x, a = A, b = B) +
rnorm(n = nrow(df), mean = 0, sd = 20) # just adding some noise

ggplot(data = df) + geom_point(aes(x = x, y = y))

mdl <- glm(data = df
, formula = y ~ log(x)
, family = gaussian(link = "log")
, start = c(0.6, 2))
summary(mdl)

exp(coef(mdl)["(Intercept)"]) # approx = A
coef(mdl)["log(x)"] # approx = B

# predictions of log(E[y])
df[, c("pred_log", "se_log")] <- predict.glm(object = mdl
, newdata = df
, se.fit = TRUE)[c("fit", "se.fit")]
# pred and se are on log scale
df$$pred_real <- exp(df$$pred_log)
df$$L_log <- df$$pred_log - 1.96 * df$$se_log df$$U_log <- df$$pred_log + 1.96 * df$$se_log
df$$L_real <- exp(df$$L_log)
df$$U_real <- exp(df$$U_log)

ggplot(data = df) +
geom_point(aes(x = x, y = y)) + # raw data
geom_ribbon(aes(x = x, ymin = L_real, ymax = U_real) # confidence bands
, alpha = 0.2
, fill = "blue") +
geom_line(aes(x = x, y = pred_real)) + # predicted curve
geom_line(aes(x = x # actual curve
, y = f(x = x, a = A, b = B))
, colour = "green")


Black points: actual data. Shaded blue interval is confidence interval on the mean for various $$x$$. Black line is the fitted curve, green line is the true curve.