# Constructing block matrices with applications to mixed models

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); end; 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!

1. Caitlyn
Posted June 14, 2013 at 2:48 pm | Permalink

I recently had a similar situation where I used the Kronecker product to generate some data. Specifically, I needed to generate a series of bivariate normal datasets. Rather than specifying them individually, I used the Kroenecker product to create a large variance-covariance matrix with the block diagonals corresponding to my bivariate matrix. However, some of this efficiency has been lost because I have to use the BY statement in Proc Mixed when I analyze the datasets. Could you recommend a way to specify the repeated statement that would allow me to process the matrix in a single step(i.e. without the BY statement)? My original code looked like this:

proc mixed data=SIMDATA method=reml ;
class SubjectID Visit Trt Sim_ID;
by Sim_ID;
model Mu = Visit Trt Visit*Trt / s ;
repeated / type=ar(1) sub=zSub_StudyID ;
lsmeans Visit*Tx/diff=control('2' '1') ;
run;

I've tried changing type=Un@ar(1) but this doesn't quite give me what I'm looking for.

Rick Wicklin, PhD, is a distinguished researcher in computational statistics at SAS and is a principal developer of PROC IML and SAS/IML Studio. This blog focuses on statistical programming. It discusses statistical and computational algorithms, statistical graphics, simulation, efficiency, and data analysis. Rick is author of the books Statistical Programming with SAS/IML Software and Simulating Data with SAS.