5
$\begingroup$

I've got a model that I've developed in R, but also need to express in SAS. It's a double GLM, that is, I fit both the mean and (log-)variance as linear combinations of the predictors:

$E(Y) = X_1'b_1$

$\log V(Y) = X_2'b_2$

where Y has a normal distribution, $X_1$ and $X_2$ are the vectors of independent variables, and $b_1$ and $b_2$ are the coefficients to be estimated. $X_1$ and $X_2$ can be the same, but need not be.

I can fit this in R using gls() and the varComb and varIdent functions. I've also written a custom function that maximises the likelihood using optim/nlminb, and verified that it returns the same output as gls.

I would now like to translate this into SAS. I know that I can use PROC MIXED:

proc mixed;
  class x2;
  model y = x1;
  repeated /group = x2;
run;

However, this only gives me what I want if I have 1 variable in the /GROUP option. If I enter 2 or more variables, MIXED can only handle this by treating each individual combination of levels as a distinct group (that is, it takes the cartesian product). For example, if I have 2 variables in $X_2$, with 3 and 4 levels respectively, MIXED will fit 12 parameters for the variance. What I want is for the log-variance to be additive in the variables specified, ie 6 parameters.

Is there a way of doing this in MIXED or any other proc? I could probably code something in NLP, but I'd really prefer not to.

$\endgroup$
4
$\begingroup$

Rather than the code you present, I assume you're doing something like class x21 x22 followed by the repeated clause group=x21*x22 to end up with 12 parameters. This is the only option I'm aware of within SAS, i.e. I don't think you can get a straightforward stratification of the variation across the combined levels of two variables. But, unless I'm missing something, you should be able re-code your variable into a combined x3 to achieve this.

Also, just to check -- I assumed that your x2 variable is prepared such that it holds the appropriate exponential values for your log-linear variance model. Otherwise, you may want to check out the local option to the repeated statement. For example, the following model in R:

wgt <- varExp(form =~ x)
fit <- gls(y ~ x, weights=wgt, data=dat)

corresponds to something like the following in SAS:

proc mixed data=dat;
  model y = x / solution;
  repeated / local=exp(x);
run;

The covariance parameter estimates will differ by a factor of 2 between R and SAS. This is because the R version models $var(\epsilon_{ij}) = \sigma^2 e^{ 2bx_{i} }$ (see ?varExp). I'm guessing you're processing your input and not using varExp because this isn't easy to accomplish in R if you have factor variables, but I thought I'd mention it in case it's helpful.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.