I would like to model heteroscedasticity or better, the standard deviation of the mean as a function of x and retrieve the coefficients that says something about how standard deviation behaves with respect to x.

here is some data:

n <- 1000
x <- runif(n)
y <- 1 - 2*x + x*rnorm(n)
r <- lm(y~x)

so the standard deviation of the mean (linear regression line of -2*x ) seems to be a function of x + x^2; how can I retrieve the coefficients of the standard deviations in R?

  • $\begingroup$ You can pass a variance model to the gls function in package nlme. Obviously, only common variance structures are implemented (a quadratic function is not common). $\endgroup$
    – Roland
    Dec 19, 2016 at 7:18
  • $\begingroup$ See also this tag: gamlss $\endgroup$ Jul 8, 2020 at 13:37

1 Answer 1


First, you need to formally define the model you want to estimate. Based on your example, you'll want to set a model like this: $$Y_i = \beta_0 +\beta_1 X_i + \epsilon_i$$ where $\epsilon_i \sim N(0, X_i^2 \sigma^2)$ and each error is iid.

Similar to a Simple Linear Regression model, there are three parameters to estimate, $\beta_0, \beta_1$ and $\sigma^2$. However, the model we've set up is a special case of a Generalized Linear Model. This means that estimating parameters can't be done as easily - at least, there are no simple formulas for them. The most common approach used to estimate the parameters is to emply Maximum Likelihood Estimation. In other words, find the three values which maximize the likelihood function.

This can easily be implemented in R:

# Define the data
n <- 1000
x <- runif(n)
y <- 1 - 2*x + x*rnorm(n)

# Define the negative log-likelihood function
LL <- function(theta) {
  fvec <- dnorm(y - (theta[1] + theta[2] * x), sqrt(x^2 * theta[3]^2))
  loglike <- sum(-log(fvec))


# Define initial parameter guesses and then optimize the parameters
theta0 <- c(0, 0, 0)
opt0 <- optim(theta0, LL)
theta1 <- opt0$par

Unfortunately, this numerical estimator won't have nice statistical properties like unbiasedness or even consistency. But it'll do a passable job of fitting the model and allow you to avoid overthinking things. However, if you want the best possible estimator for your model, check out https://en.wikipedia.org/wiki/Heteroscedasticity-consistent_standard_errors.


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.