Constructing block matrices with applications to mixed models

November 7, 2012

(This article was originally published at The DO Loop, and syndicated at StatsBlogs.)

The other day I was constructing covariance matrices for simulating data for a mixed model with repeated measurements. I was using the SAS/IML BLOCK function to build up the "R-side" covariance matrix from smaller blocks. The matrix I was constructing was block-diagonal and looked like this:

The matrix represents a covariance matrix for four individuals where each individual has three repeated measurements and where the measurements for each individual exhibit a compound symmetric covariance structure.

To construct the block-diagonal matrix, the first step is to construct the little 3x3 compound symmetric matrix, as follows:

proc iml;
k=3;                   /* number of measurements per individual */
b = j(k,k,1) + 2*I(k); /* residual cov=1; common cov=2 */
print b[label="cov"];

To simulate data for multiple subjects, it is common to build R, a big block-diagonal matrix such that each block equals b. This matrix is shown at the top of this post. I initially used the BLOCK function to constructs the covariance matrix, as follows:

s=4;                   /* number of individuals */
R = b;                 /* block diagonal: one kxk block for each indiv */
do i = 2 to s; 
   R = block(R,b); 
print R;

That works, and the BLOCK function is good to know about, but I recently realized that there is a more efficient way to construct a block-diagonal matrix. Calling the BLOCK function in a DO loop is unnecessary! The SAS/IML language supports the direct product (Kronecker product) operator, @. This is not an operator that I use every day (or even every year!), but it turns out that this operator makes it super-easy to construct block-diagonal matrices. The following statement computes exactly the same block-diagonal matrix:

/* block-diagonal matrix with s blocks, each block equals b */
R = I(s) @ b;

I always think it is cool to discover new ways to use the SAS/IML language to avoid writing a DO loop. In any matrix-vector language (SAS/IML, R, MATLAB,...) that is a key to efficient performance. This trick to construct a block-diagonal matrix efficiently is the latest entry on my ever-growing list!

tags: Efficiency, Matrix Computations, Statistical Programming

Please comment on the article here: The DO Loop

Tags: , , ,