Replicating R model in SAS

Question

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.

Answer 1

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.