Calcite | Level 5

## Generalized Least Squares

Greetings all!

I am looking for a good way to fit a generalized least squares model (as I hear it is one of the best ways to deal with heteroskadasticity), but am not sure of the best way to do it.

Searching SAS documentation only turned up proc surveyreg, but I've heard that proc mixed and proc genmod can do it too. In fact, I've heard that proc genmod might be the preferred method for GLS, but I can't seem to figure out how to use it to fit a GLS model.

So here are the questions:
- What is the best method to use to fit a generalized least squares model?
- If proc genmod is the preferred method, how is it done?

Thanks!
Pyrite | Level 9

## Re: Generalized Least Squares

I don't know how you performed your search, but there are a number of SAS procedures which employ generalized least squares estimation methods. Since you don't tell us anything about the problem other than that you want some way to deal with heteroskedasticity, it is a bit difficult to advise on what methods to employ. However, assuming that you are fitting a regression model in which the residual error is porportional to a power of the mean, then I would suggest looking at the MIXED procedure. Such estimation is discussed in the documentation of the REPEATED statement.

Actually, you might find it easier to use the NLMIXED procedure to estimate a model in which the residual variance is proportional to a power of the mean. NLMIXED code that fits a model in which the variance is proportional to the mean can be fit with code like the following:

proc nlmixed data=mydata;
mu = b0 + b1*x1 + b2*x2 + ...;
model y ~ normal(mu, mu*exp(2*log_sig));
estimate "Var(eps) for mu=1" exp(2*log_sig);
run;

There will be problems with the above model if your parameters produce a value of mu which is non-positive (zero or negative). A generalization of the above model which is consistent with the discussion of the POM model in the PROC MIXED documentation would be to fit the model:

proc nlmixed data=mydata;
mu = b0 + b1*x1 + b2*x2 + ...;
model y ~ normal(mu, (mu**theta)*exp(2*log_sig));
estimate "Var(eps) for mu=1" exp(2*log_sig);
run;

I made the assumption above that you were interested in a model in which the heteroskedasticity produced a residual variance which was proportional to a power of the mean. However, you might be interested in a model in which the residual variance differs according to the level of some categorical predictor variable. Such models can be handled quite easily with the MIXED procedure. Again, look at the documentation of the REPEATED statement and pay particular attention to the GROUP= option.

HTH,

Dale
Discussion stats