1
$\begingroup$

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?

$\endgroup$
8
  • 1
    $\begingroup$ 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$? $\endgroup$
    – Alex J
    Commented Feb 29 at 0:02
  • $\begingroup$ @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. $\endgroup$
    – Sorfosh
    Commented Feb 29 at 0:06
  • $\begingroup$ Edited, my apologies. I think it is indeed confidence bands. $\endgroup$
    – Sorfosh
    Commented Feb 29 at 0:12
  • $\begingroup$ 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 $\endgroup$
    – Alex J
    Commented Feb 29 at 0:14
  • $\begingroup$ 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 $\endgroup$
    – Alex J
    Commented Feb 29 at 0:16

1 Answer 1

1
$\begingroup$

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
                                             , type = "link"
                                             , 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.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.