# Conditionally heteroskedastic linear regression: How can I model variance from given predictors?

As a motivating example, consider this hypothesis: A friend and I were discussing how it might be the case that, the more one person is given control in the process of making a movie, the higher the variance of the quality of that film. The more hands that go into making a film, the safer of a bet it is—but the less likely it will be that it is great.

Here are some simulated data that could characterize this hypothesis:

set.seed(1839)
x <- runif(1000, 0, 100)
y <- rnorm(1000, 0, x)
dat <- data.frame(x, y)

library(ggplot2)
ggplot(dat, aes(x = x, y = y)) +
geom_jitter() +
theme_light() +
labs(y = "Quality of Film",
x = "Amount of Artistic Control in One Person")


Note that here, I am not interested in predicting the expected value, but the expected variance of $y$ from $x$.

I am familiar with how one can create models for both location and dispersion in beta regression and how one can create models for both location and scale in conditionally heteroskedastic truncated regression.

However, is there a way to model conditional heteroskedasticity in an ordinary least squares regression like this? Are there regression models that exist where one can model the variance of $y$ as a function of $x$? When I try looking for these models, I usually end up with people wanting to use sandwich estimators or weighted least squares or something.

I don't want to treat heteroskedasticity as a problem with the data, but something to be explicitly modeled as a function of predictors $X$. Is there a way to do this with linear regression?

• What's the problem you see with weighted least squares? I don't see how it implies "a problem with the data"; if your points have different variance of a known form (up to an unknown scale), but the other assumptions of regression hold, it's the optimal way to estimate the parameters of the model. – Glen_b Jan 9 '18 at 4:02
• @Glen_b How does a weighted least squares model the conditional variance? I don't follow how weighting certain observations can do this. In the motivating example above, I want a coefficient that will tell me, "For every unit increase in $x$, the variance in $y$ increases by $b$," or something similar (like you can get in beta regression with the dispersion submodel, for example). – Mark White Jan 9 '18 at 4:08
• Consider for example the gamma GLM. I know the form of the relationship -- variance $\propto$ mean$^2$ (spread proportional to mean). If I compute my model's estimate of the conditional mean, and I have a suitable estimate of the dispersion parameter (the scaling constant), I also have an estimate of the conditional variance. Indeed if both dispersion and the mean can be estimated by maximum likelihood (as they can for a Gamma) then the variance estimate is ML also. – Glen_b Jan 9 '18 at 5:08
• You can set up a relationship between the parameters but you have to be careful how you do it. For example, if you're prepared to have your mean line go through the intercept, it is easy to do. Let $\sigma\propto\mu$, say and then $Y= bx + e$ implies $Y = bx + bx\eta$ for $\eta$ having constant variance. Now divide through both sides by $x$, and you have $Y^* = b + \xi$ where now $\xi$ has constant variance; you can estimate $b$ (compute the mean of $Y^*$) and the sd of $\xi$ (e.g. compute the sd of $Y^*$) and hence get an estimate of the variance of $Y$ as a function of $x$. ...ctd – Glen_b Jan 9 '18 at 5:35
• ctd... (This all requires both $x$ and $b$ be non-negative.) ... So sometimes you can make such a model work, but too often I have seen people set up models where the implied variance (or implied standard deviation, depending on the model) of the data can be negative somewhere where you can observe data. It must be done with some thought. [Note that with this model if the observed $Y$ can't/won't be negative, a normal error doesn't really make sense; a construction with positively distributed variables may do better.] – Glen_b Jan 9 '18 at 5:36

## 3 Answers

Simultaneous modelling of mean and variance using double generalized linear models

The emphasis of gamlss is obviously on generalized additive models (GAMs). The general of idea of simultaneously modelling the mean and variance using generalized linear models (GLMs) rather than GAMs has been around much longer, and you might find the dglm package (Double GLMs) a simpler entry point for this sort of analysis.

I came up with the idea of double generalized linear models as part of my 1985 PhD thesis. In the normal case, the first publication with the same idea is Aitkin (1987).

For your example data we have:

> set.seed(1839)
> x <- runif(1000, 0, 100)
> y <- rnorm(1000, 0, x)
> library(dglm)
> fit <- dglm(y~x, dformula = ~x)
> summary(fit)

Call: dglm(formula = y ~ x, dformula = ~x)

Mean Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.3273     1.2319  -0.266    0.791
x            -0.0336     0.0467  -0.719    0.472
(Dispersion Parameters for gaussian family estimated as below )

Scaled Null Deviance: 1001 on 999 degrees of freedom
Scaled Residual Deviance: 1000 on 998 degrees of freedom

Dispersion Coefficients:
Estimate Std. Error z value  Pr(>|z|)
(Intercept)   4.9319    0.08996    54.8  0.00e+00
x             0.0508    0.00157    32.4 7.99e-230
(Dispersion parameter for Gamma family taken to be 2 )

Scaled Null Deviance: 2082 on 999 degrees of freedom
Scaled Residual Deviance: 1419 on 998 degrees of freedom

Minus Twice the Log-Likelihood: 10296
Number of Alternating Iterations: 6


You can see from the summary results that there is no significant trend for the means but a very highly significant increasing trend for the variances. The dispersion coefficients are exactly twice the sigma coefficients returned by gamlss (in Stefan's answer) because dglm is modeling the variance rather than the standard deviation.

We can do even better by using $\log x$ as the predictor for the variance. You simulated the data to have sd($y_i$)$=\sigma_i=x_i$ so, on the log-scale, the true variance model is $\log\sigma^2_i=2\log x_i$. We can fit this by:

> fit <- dglm(y~x, ~log(x))
> summary(fit)

Call: dglm(formula = y ~ x, dformula = ~log(x))

Mean Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.0189     0.0410  -0.461    0.645
x            -0.0346     0.0323  -1.069    0.285
(Dispersion Parameters for gaussian family estimated as below )

Scaled Null Deviance: 1001 on 999 degrees of freedom
Scaled Residual Deviance: 1000 on 998 degrees of freedom

Dispersion Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   0.0843     0.1694   0.498    0.619
log(x)        1.9850     0.0453  43.806    0.000
(Dispersion parameter for Gamma family taken to be 2 )

Scaled Null Deviance: 2068 on 999 degrees of freedom
Scaled Residual Deviance: 1188 on 998 degrees of freedom

Minus Twice the Log-Likelihood: 10079
Number of Alternating Iterations: 4


You can see that the variance trend coefficient is estimated to be 1.985 when the true value is 2. The variance intercept is estimated to be 0.0843 when the true value is 0.

How it works

The basic idea is quite simple. We first fit a linear regression to the data:

> fit.m <- lm(y~x)


Then we extract the squared deviance residuals (aka unit deviances) from the fit and fit a log-linear chisquared glm to them:

> d <- residuals(fit.m)^2
> fit.d <- glm(d~x, family=Gamma(link="log"))


Then we take the fitted values from the variance (dispersion) model and use them as weights for the mean model:

> w <- 1/fitted(fit.d)
> fit.m <- lm(y~x, weights=w)


The dglm() function simply iterates this process to convergence. It can be shown that this computes MLEs of all the parameters.

To examine the mean-model fit, you should use:

> summary(fit.m, dispersion=1)


The dispersion should be set to 1 because the variances have already been incorporated into the weights.

To examine the variance model, you should use

> summary(fit.d, dispersion=2)


The dispersion should be set to 2 because the unit deviances follow scaled chisquared distributions on 1 df, and chisquared on 1 df has squared coefficient of variation equal to 2.

The above example is for normal data, but the dglm() function works for any GLM family.

REML

The dglm package has the additional ability to use a REML method to return approximately unbiased variance estimators, unlike ML estimation which will return variance estimators that are systematically too small. For this data it makes little difference, because the number of mean parameters (2) is small compared to the number of observations (1000):

> fit <- dglm(y~x, ~x, method="reml")
> summary(fit)

Call: dglm(formula = y ~ x, dformula = ~x, method = "reml")

Mean Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.3297     1.2353  -0.267    0.790
x            -0.0335     0.0467  -0.716    0.474
(Dispersion Parameters for gaussian family estimated as below )

Scaled Null Deviance: 998 on 999 degrees of freedom
Scaled Residual Deviance: 997 on 998 degrees of freedom

Dispersion Coefficients:
Estimate Std. Error z value  Pr(>|z|)
(Intercept)   4.9428    0.09029    54.7  0.00e+00
x             0.0507    0.00157    32.2 3.96e-227
(Dispersion parameter for Gamma family taken to be 2 )

Scaled Null Deviance: 2082 on 999 degrees of freedom
Scaled Residual Deviance: 1419 on 998 degrees of freedom

Minus Twice the Log-Likelihood: 10296
Number of Alternating Iterations: 7


References

Aitkin, M. (1987). Modelling variance heterogeneity in normal regression using GLIM. Applied Statistics 36(3), 332-339.

Smyth, G. K. (1989). Generalized linear models with varying dispersion. J. R. Statist. Soc. B, 51, 47–60.

Smyth, G. K., and Verbyla, A. P. (1999). Adjusted likelihood methods for modelling dispersion in generalized linear models. Environmetrics, 10, 696-709. http://www.statsci.org/smyth/pubs/Ties98-Preprint.pdf

Smyth, G. K., and Verbyla, A. P. (1999). Double generalized linear models: approximate REML and diagnostics. In Statistical Modelling: Proceedings of the 14th International Workshop on Statistical Modelling, Graz, Austria, July 19-23, 1999, H. Friedl, A. Berghold, G. Kauermann (eds.), Technical University, Graz, Austria, pages 66-80. http://www.statsci.org/smyth/pubs/iwsm99-Preprint.pdf

Smyth, GK, and Verbyla, AP (2009). Leverage adjustments for dispersion modelling in generalized nonlinear models. Australian and New Zealand Journal of Statistics 51, 433-448.

• brilliant! And it looks like Smyth, G. K. (1989). Generalized linear models with varying dispersion. J. R. Statist. Soc. B, 51, 47–60. Might be a good place to start, reading-wise? Or is there a review paper/book you would suggest for covering these double GLMs? – Mark White Jan 9 '18 at 6:39
• @MarkWhite You might start with Murray Aitkin's 1987 article or my 1999 article in Environmetrics. – Gordon Smyth Jan 9 '18 at 7:03
• This is indeed a simpler entry point and nicely explained! Until now I have not heard about the dglm package. Great answer! – Stefan Jan 11 '18 at 1:00

If I understand you correctly, this can be done using the gamlss() function in the gamlss package.

Given your example, you can model the scale parameter $\sigma$ as follows:

fit <- gamlss(y~x, sigma.formula = ~x, data = dat, family = NO)
> summary(fit)
******************************************************************
Family:  c("NO", "Normal")

Call:  gamlss(formula = y ~ x, sigma.formula = ~x, family = NO, data = dat)

Fitting method: RS()

------------------------------------------------------------------
Mu link function:  identity
Mu Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.32765    1.23130  -0.266    0.790
x           -0.03359    0.04671  -0.719    0.472

------------------------------------------------------------------
Sigma link function:  log
Sigma Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.4663263  0.0533045   46.27   <2e-16 ***
x           0.0254100  0.0009709   26.17   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

------------------------------------------------------------------
No. of observations in the fit:  1000
Degrees of Freedom for the fit:  4
Residual Deg. of Freedom:  996
at cycle:  3

Global Deviance:     10296.29
AIC:     10304.29
SBC:     10323.92
******************************************************************


See also Table 1 in the link above for all other available distributions.

• Ah, marvelous. I should have been aware of this package before now. Looks like they cite anything else I should need, too. For anyone that stumbles upon this in the future, it looks like the authors have a nice book on this package, too: amazon.com/Flexible-Regression-Smoothing-GAMLSS-Chapman/dp/… – Mark White Jan 9 '18 at 5:09
• @MarkWhite Indeed the book is great! I bought it myself because I have to deal a lot with zero inflated proportion data (zero inflated beta). At the same time it is also quite dense and complex. They also have a great website but for some reason it's down. – Stefan Jan 9 '18 at 5:12

@Stefan gave rich family of models with the use of gamlss so my response here is more of an extended comment.

For some models with non-constant variance one can just do a simple "rearrangement" of the model and use simpler functions such as lm.

Suppose the model is

$$y=a+b x + x \epsilon$$

with $\epsilon \sim N(0,\sigma^2)$. We can just divide everything by $x$ and end up with a "standard" linear model with constant variance:

$${y \over x} = {a \over x} + b + \epsilon$$

We now have a typical linear model with $y/x$ as the dependent variable, $1/x$ is the independent variable with "slope" $a$ and intercept $b$ and the variance is constant.

• I like this answer! I never thought about it that way! – Stefan Jan 11 '18 at 1:03