# I think standard deviation of y is related to size of x. How do I create a model for this / test this?

I have a sample of data $$(x_i, y_i)$$. I hypothesize that $$y_i$$ is not dependent on $$x_i$$, but the standard deviation of $$y_i$$ depends on $$x_i$$

More concretely, say I assume $$\textrm{Var}(y_i | x_i) = f(x_i)$$. This could be say $$cx_i^2$$ for some value of $$c$$. The question is that how do I visualize this and how do I model this? Right now because I do not have multiple value of $$y_i$$ for same value of $$x_i$$, it is not possible to compute the standard standard deviation.

One way to model this type of relationship is by using generalized least squares (GLS).

GLS allows you to specify variables which have a relationship not just with the mean of the response, but also the variance.

Here is how you can do that in R:

library("nlme")
GLS <- gls(y ~ x, weights = varPower(form = ~ x), data = data)

• I'm not familiar with gls syntax; is varPower() here saying that the variance of y is a function of the mean of x and not its variance (as the question asks for)? If so, can it be changed to accommodate that?
– mkt
Commented Jul 17 at 21:00
• This is a pretty good starting place. Thanks a lot. Does varPower(form=~x) mean the variance is linear or x? Commented Jul 18 at 8:58
• @Lost1 varPower estimates a power-relationship between $x$ and the variance of $y$. There are multiple alternatives described in the help page for varClasses Commented Jul 18 at 9:23

You are asking models of the variance function, or something in the same vein.

There has to be some sort of smoothing in the variance model, because we don't expect anything useful can be done if each $$Var(y_i\mid x_i)$$ has a completely different parameter. The $$c$$ in your example has done that smoothing. If the # of parameters is low compared to sample size, usual estimation techniques (e.g., mle) should work well. Visualization following estimation is then straightforward, e.g., a curve from your estimated variance function overlapped with squared residuals vs $$x$$.

However, another quick nonparametric check is to group $$x_i$$ into bins so that you have multiple $$Y$$ for each bin of $$X$$.

PS: Better to have a good model on $$E(Y_i\mid X_i)$$ first, before spending too much effort in the variance model.

• PS: Be sure to have a good model on $E(Y_i\mid X_i)$ first, before spending too much effort in the variance model. Commented Jul 17 at 16:09
• You can edit your post and add the postscript. Welcome to CV! Commented Jul 17 at 17:41

Generalized additive models for location, scale and shape (GAMLSS) have been developed to do exactly that. As the name implies, they allow not only for explicit modelling of the location but also of scale and shape of the conditional distribution. The gamlss package in R makes this easy.

For a simple conditional normal model with a log-link for the standard deviation you could write:

\begin{align} \mathbf Y &\sim N(\mathbf \mu, \mathbf \sigma) \\ \mathbf \eta_1 &= g_1(\mathbf \mu) = \mathbf \mu = \mathbf X_1\mathbf\beta_1 \\ \mathbf \eta_2 &= g_2(\mathbf \mu) = \log(\mathbf \sigma) = \mathbf X_2\mathbf\beta_2 \\ \end{align} where $$\mathbf Y$$ is the response vector, $$\mathbf X_1$$ and $$\mathbf X_2$$ are the design matrices for the mean and standard deviation and $$\mathbf \beta_1$$ and $$\mathbf \beta_2$$ are the coefficient vectors for the mean and standard deviation.

In R you'd fit such a model like this:

# Load the package
library(gamlss)

# Fit the model
mod <- gamlss(
y~x1 + x2 + x3
, sigma.formula = ~x1 + x2 + x3
, family = "NO"
, data = dat
)

# Inspect fit
plot(mod)

# Look at summary
summary(mod)


In the function gamlss, the first formula is the linear predictor for the mean whereas the formula called sigma.formula specifies the linear predictor for the standard deviation with a log-link.

You can extend this model by adding splines (e.g. pb) or other nonlinear terms such as polynomials to the formulas. If the conditional distribution has more parameters such as the skew $$t$$ distribution for example, you can further model the skewness and kurtosis using tau.formula and nu.formula, respectively.

• I found out about this after asking this question. I ended up using xgboostlss, which can fit to a small subset of gamlss function using regression tree method. This is a good answer. I am using it for forecasting a distribution. Is there a way for gamlss to support this with cross validation in R? Commented Aug 9 at 19:46
• @Lost1 Thank you. The package contains the gamlssCV function that allows for cross-validation. In the reference book, they use this to do model selection. If you have a training and validation dataset, the functions gamlssVGD or getTGD could be used. Commented Aug 9 at 19:54