Linear regression with changing variance

Question

I want to perform linear regression on some data. For every value of x, the data values are distributed normally across y, around some mean. However, the variance increases linearly as x increases. I made this example graph:

Blue is the regression line, red are data points, black shows the normal distribution, and green visualizes the variance increasing.

How can I calculate a regression for the change in variance, while also performing a linear regression of the data? The data is heteroscedastic, and I've read up on methods for doing linear regression on such data. However, I haven't found anything on estimating the actual change in variance of the data.

I haven't studied stats rigorously, so any simple explanations or resources I could look at further would be appreciated.

More Details:

The original dataset follows $y = a/x + b$. The variance as $x$ changes follows a similar model $s^2 = c/x + d$. I transformed the data using $x' = 1/x$ to make the data linear (just to simplify the problem). Here is a sample graph (left is transformed, right is original):

Answer 1

This sounds like a special case of heteroscedasticity.

There are two issues:

1. What estimator should you use in the presence of heteroscedasticity?
2. How should you calculate your standard errors?

The most straightforward thing to do is run a regular regression but use heteroscedastic robust standard errors. As @Glen_b suggests in the comments though, you probably can do better than this by efficiently exploiting known structure on your problem.

What estimator to use?

• You could just run a normal regression.

  • In the presence of heteroscedasticity, the regular ordinary least squares (OLS) estimator is still consistent. In layman's terms, OLS still works given enough data. But OLS is not efficient.

• You could run weighted least squares, an application of generalized least squares. The loose idea is to give more weight to observations with low variance error terms.

  • Since you probably don't know ex-ante how the variance of the error term varies with $x$, you probably have to do something like feasible gls.

If you run a regular OLS regression, you should not use the usual standard errors based upon assumptions of homoscedasticity. Instead you should use heteroscedastic robust standard errors. Any stats package can do this.

Answer 2

Your data violate the assumption of homoscedasticity. You can use a regression method that produces standard errors that are robust to heteroscedasticity. What software are you using to run your regression? If you are using R, you can use the sandwich package to estimate robust standard errors.

Answer 3

Machine learning guy here -- as much as I love the stats folks ML usually wins the day in real-world applications, and heteroscedasticity is one such common occurrence. A more general solution that is non-parametric and works for this and other (even highly) nonlinear regression problems is to just use quantile regression with a bagged decision tree. You can follow some links at Matlab to learn more (their examples are easily transferable to Python or R).

Answer 4

To address the question of how to "estimate the actual change in variance of the data", let me reiterate my answers to this related question: nothing stops you from modelling the variance of your model with something (slightly) more complicated than a constant $\sigma^{2}$.

It helps to take a step back in think about the logic of regression: you're replacing a constant mean with a function of $x$, such that $y_{i} \sim N(\mu,\sigma^{2})$ becomes $y_{i} \sim N(\alpha + \beta x_{i}, \sigma^{2})$, where I'm assuming normality to keep it simple. You can take it a step further an assume $y_{i} \sim N \left(\alpha + \beta x_{i}, \sigma_{i}^{2} \right)$, and model the heteroscedasticity, for example, as $\sigma_{i}^{2} = \sigma^{2} x_{i}^{\delta}$. Your Gaussian log-likelihood function then becomes $$ \ell(\alpha, \beta, \sigma^{2}, \delta) = -\frac{n}{2} \log ( 2 \pi) - \frac{1}{2} \sum_{i=1}^{n}\log \left( \log \sigma^{2} + \delta \log x_{i} \right) -\frac{1}{2}\sum_{i=1}^{n} \left( \frac{y_{i} - \alpha - \beta x_{i}}{\sigma x_{i}^{\delta/2}} \right)^{2} $$ This is essentially the model discussed in Kmenta, Elements of Econometrics, pp. 279-292.