How to implement a variance function in R?

Question

I want to estimate a multivariate variance function in R. That is, I want to allow the variance (as well as the mean) to vary according to some set of independent variables.

In this particular case, I want to estimate the effects of a set of typical demographic covariates (age, race, education) on the variance of logged wages.

What is a good way to implement this in R? Is there a package that simplifies this?

It may be that this is only a search away - but having searched on the R help pages, Google, Rseek, and StackOverflow, I can't find anything relevant under "variance function" or similar.

Any suggestions gratefully received.

---

Thanks for your responses -- I will try to clarify my question.

I am working in a maximum likelihood framework. I can code this by hand from the log-likelihood, but the real data set has a lot of variables and "optim" is very slow, so I would like to find a package in R that makes this more computationally efficient.

I start with the log-likelihood for a basic OLS regression: $$ \text{ln }L = \sum (-\frac{1}{2} (\text{ln }\sigma^2 - \frac{(y - xB)^2}{\sigma^2})) $$ Then I relax the assumption of constant variance (homoskedasticity) and redefine the variance as: $$ \sigma^2 = exp(Z*\gamma) $$ where $Z$ is the matrix of variables affecting $\sigma^2$. (Exponentiate so that you don't end up with $\sigma^2$ less than zero.) When I substitute the reparameterization of $\sigma^2$ into the original log-likelihood and code the new log-likelihood function in R, I get this:

ll.normal.vary <- function (par, X, Y, Z) {
  beta  <-par[1:ncol(X)]
  gamma <- par[(ncol(X)+1):(ncol(X)+ncol(Z))]   
  -1/2* sum((Z %*% gamma) + ((Y - X %*% beta)^2)/exp(Z %*% gamma))
}

Then I optimize:

v.optim1 <- optim (par = start1, fn=ll.normal.vary, X=x.mat, Y=y.vec, Z=z.mat, 
                   method = "BFGS", hessian = F, control = list(fnscale = -1))
v.optim1$par
v.optim1$value

Here are some sample data if you want to test it:

var1   <- c(0,0,0,1,1,0)
var2   <- c(.28, .07, -.05, .38, .08, -.1)
var3   <- c(-.11, -.17, -.17, -.05, .1, -.01)
x.mat  <- cbind(var1, var2, var3)
y.vec  <- c(.46, .77, .49, .59, .60, .44)
z.mat  <- cbind(var1, var2) 
start1 <- rep(0.1, ncol(x.mat)+ncol(z.mat))

Thanks again for any tips.

Answer 1

Iteratively re-weighted least squares is probably the answer. Start with a set of weights all equal to 1, or equal to the inverse of your best guess of the variance at each point. Use weighted least squares to fit the model with that set of weights. Then re-estimate your weights based on the latest set of parameters; and refit your model with that set of weights. Continue until you get convergence. This is basically the method commonly used for fitting generalized linear models, where the variance at any point depends on the link function and distribution of the response.

There are certainly R functions that do this for you, but if you want to play around with different models of how variance relates to your various parameters you may find it just as easy to do this yourself.

Alternatively, if you have a fairly standard variance model, you may find that just the glm() function, with an appropriate link function and distribution, works for you. The "quasi" family may be particularly useful as it allows you to specify both the link function and the relationship of variance to the response, with minimal assumptions about the other moments of the response.