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

Question

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.

Answer 1

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)

Answer 2

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.

Answer 3

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.