SAS/IML: how to use individual variance components in RANDNORMAL

Question

This is a programming question, but I'll give you a little of the stats background first. This question refers to part of a data sim for a mixed-effects location scale model (i.e., heterogeneous variances). I'm trying to simulate two MVN variance components using the RANDNORMAL function in IML. Because both variance components are heterogeneous, the variances used by RANDNORMAL will differ across people. Thus, I need IML to select the specific row (e.g., row 1 = person 1) and use the RANDNORMAL function before moving onto the next row, and so on.

My example code below is for 2 people. I use DO to loop through each person's specific variance components (VC1 and VC2). I get the error: "Module RANDNORMAL called again before exit from prior call." I am assuming I need some kind of BREAK or EXIT function in the DO loop, but none I have tried work.

PROC IML;
    ColNames = {"ID" "VC1" "VC2"};
    A = {1 2 3, 
         2 8 9};
    PRINT A[COLNAME=ColNames];
    /*Set men of each variance component to 0*/ 
    MeanVector = {0, 0};
    /*Loop through each person's data using THEIR OWN variances*/
    DO i = 1 TO 2;
        VC1 = A[i,2];
        VC2 = A[i,3];
        CovMatrix = {VC1 0, 
                     0   VC2};
        CALL RANDSEED(1);  
        U = RANDNORMAL(2, MeanVector, CovMatrix);
    END;
QUIT;

Any help is appreciated. Oh, and I'm using SAS 9.4.

You're generating 2 rows total, or 2 rows per person? (Also, if you're using IML, use the right tag - Rick answers the [tag:sas-iml] questions but I don't know if he looks at just [tag:sas].) — Joe, Sep 18 '14 at 21:29
Each row is one person. Oh, and they wouldn't let use a sole "IML" tag because I don't have 1500 reputation... — Ryan W., Sep 19 '14 at 13:46

score 1 · Accepted Answer · answered Sep 18 '14 at 21:42

1

You want to move some things around, but mostly you don't want to rewrite U twice: you need to write U's 1st row, then U's 2nd row, if I understand what you're trying to do. The below is a bit more efficient also, since I j() the U and _cv matrices rather than constructing then de novo every time through the loop (which is slow).

proc iml;

  a = {1 2 3,2 8 9};
  print(a);
  _mv = {0,0};
  U = J(2,2);
  _cv = J(2,2,0);
  CALL RANDSEED(1);   

  do i = 1 to 2;
    _cv[1,1] = a[i,2];
    _cv[2,2] = a[i,3];
    U[i,] = randnormal(1,_mv, _cv);
  end;
  print(u);
quit;

answered Sep 18 '14 at 21:42

Joe

62,789
6
49
67

This answer is completely correct and worked great. With that said, my example code was programmed similar to what I think the steps in IML should look like, so Rick's code fix below kept that consistent for me (less mental muscle required). Your solution will inevitably be very useful for others. Thanks! – Ryan W. Sep 19 '14 at 14:53
I highly recommend reading Rick's other blog posts - he has some good examples of how certain things should be done, including why to use the `j()` how I do above versus regenerating them like you do each time through. (But, I'm not really an IML programmer, so largely my answer is based on that reading!) – Joe Sep 19 '14 at 15:38
After comparing both solutions and ramping up to 1000 people, I found your code to be much more efficient. This is not Rick's fault, it's my own crappy programming. It just took me some time to figure out exactly what you're doing. I switch between answer check marks because I want to select both responses as they get me to the same place. Anyways, I really appreciate it! Thanks again for your patience. – Ryan W. Sep 19 '14 at 15:53

score 1 · Answer 2 · answered Sep 19 '14 at 14:02

1

Your mistake is the line

CovMatrix = {VC1 0, 0   VC2}; /* wrong */

which is not valid SAS/IML syntax. Instead, use @Joe's approach or use

 CovMatrix = (VC1 || 0) // (0 || VC2);

For details, see the article "How to build matrices from expressions."

You might also be interested in this article that describes how to carry out this simulation with a block-diagonal matrix: "Constructing block matrices with applications to mixed models."

answered Sep 19 '14 at 14:02

Rick

1,210
6
11

This is great, thanks Rick! Your first linked blog post does a great job at explaining what I did wrong. Honestly, I just couldn't articulate my question well enough in any Google search to get to that blog post. I had only used curly brackets for matrices in the past, but I had no idea that curly brackets were only good for fixed values (er, literals). – Ryan W. Sep 19 '14 at 14:49

SAS/IML: how to use individual variance components in RANDNORMAL

2 Answers2