simulation study using AIC for selecting covariance structre of errors in mixed model (proc mixed)

[fatemeh]

Hi all,

how can i do simulation study for choosing the best covariance structure of errors using AIC and BIC in mixed models (using proc mixed) .appreciate any help.

[Rick_SAS]

Are there any more details you can provide? What covariance structure(s) do you intend to simulate? 

A complete description of the simulation study would be quite lengthy. If you plan to do a lot of simulation, I recommend Simulating Data with SAS, which provides many details and examples of how to efficiently run simulation studies in SAS.

Briefly, I suggest you read "Simulation in SAS: The slow way or the BY way" which outlines the main ideas. You might also want to read "Simulate many samples from a logistic regression model," which provides a nontrivial regression example.

Overall, you design would look like this 

1. Choose a covariance structure for the simulated data

2. Simulate B samples using the BY-group method. (B large, but use smaller B while debugging).

3. Use PROC MIXED with a BY variable to analyze the data.  Save the AIC and BIC statistics in a data set. Do this for each covariance structure in your study; each structure requires its own PROC MIXED step.

4. Examine and compare the distributions of AIC and BIC. Do the AIC (BIC) values for one structure tend to be less than the AIC scores for a different structure?

[fatemeh]

Thanks for your response. For examine and compare the distribution of AIC what procedure should I use? Proc univariate or proc means?

[Rick_SAS · Accepted solution]

Either would work.PROC MEANS is more concise. PROC UNIVARIATE enables you to produce comparative histograms and fit distributions (e.g., normal) to the approximate sampling distributions.

You need to use the procedure(s) that address your research questions. For example, if you are hypothesizing that the mean AIC value is lower for the AR(1) structure than for the CS structure, you might use PROC TTEST to run a two-sample t test for the means.

[fatemeh]

Hello ,

i have a question about the estimating sampling distribution of the variance of random effects of the data(varCS and varR)  that i was simulating using  repeated measures data having fixed effect and random error effects.i used proc univariate to see the sampling distributions but why my estimations  are not close to the variance of CS and residuals ? R is a compound symmetry matrix. bellow is my code

proc iml;
FixedVar = "_X1":"_X48";
use All;
read all var FixedVar into X;
close All;
use SolutionF; read all var {Estimate} into beta;
close;
eta = X*beta;                     
use CovParms; read all var {Estimate} into var; close;                 
varCS=var[1];
varR=var[2];             /* Variance of residuals*/
k=15;
s=17;
b=varCS*j(k,k,1 )+varR*I(k);
R=b;
do i=2 to s;
R=block(R,b);
end;
R = I(s) @ b;*print R;
zer = j(1,k*s,0);*print zer;
call randseed(1234);
eps= j(nrow(x) , 500);
call randgen(eps, "normal" ,zer ,R);
Y = eta +eps; * print Y;
yNames = "y1": ("y"+strip(char(500)));          
create Simdat from Y[c=yNames];  append from Y;  close;
quit;

[Rick_SAS]

I think you want eps to be MVN(0, R), in which case you should call 

eps = RandNormal(1, zer, R);

instead of 

call randgen(eps, "Normal", zer, R);

[fatemeh]

why it gives me error that  Matrix Y has not been set to a value?bellow is the changed code.

proc iml;
FixedVar = "_X1":"_X48";
use All;
read all var FixedVar into X;
close All;
use SolutionF; read all var {Estimate} into beta;
close;
eta = X*beta;                     
use CovParms; read all var {Estimate} into var; close;                 
varCS=var[1];
varR=var[2];             /* Variance of residuals*/
k=15;
s=17;
b=varCS*j(k,k,1 )+varR*I(k);
R=b;
do i=2 to s;
R=block(R,b);
end;
R = I(s) @ b;*print R;
zer = j(1,k*s,0);*print zer;
call randseed(1234);
eps= j(nrow(x) , 500);
eps=randnormal(1 ,zer ,R);
Y = eta +eps; * print Y;
yNames = "y1": ("y"+strip(char(500)));          
create Simdat from Y[c=yNames];  append from Y;  close;
quit;

[Rick_SAS]

Because your dimensions are not compatible.  I assume that eta is a column vector with N elements where N=nrow(X).  The eps vector is a row vector with dimension 1 x 255.

[fatemeh]

thanks your solution solved my problem.

[Rick_SAS]

I don't know. It is very difficult to guess what you might be doing wrong. I suggest you make sure that you are not having convergence problems. Good luck.

[fatemeh]

Hello.DR.Wicklin,

First i appreciate you for all your knowledge that you are sharing with us. I am applying proc mixed with only repeat statement and the covariance structure is sp(power). my question is about sas output. I am confused about the first covparm estimate of the output table, which is sp(pow). what is its relation with rho. To be clear, i attached the output table and my sas code below. Thanks for any help.

[Rick_SAS]

I suggest you open a new thread with an appropriate title to ask about the meaning of the parameter in the mixed model.