The LMIXED Procedure

Microarray experiments are an advanced genomic technique that is used in the discovery of new treatments for diseases. Microarray analysis enables tens of thousands of genes to be detected in a single DNA sample. A microarray is a glass slide or membrane that has been spotted or "arrayed" with DNA fragments (oligonucleotides) that represent specific genes. The response of the gene that is detected by a spot is proportional to the intensity of fluorescence associated with that spot. These gene responses can indicate associations with disease conditions, but they can also be affected by systematic biases and different treatments such as sex and genotypes. Statistical models for microarray data attempt to assess the significance and magnitude of gene effects across treatments while adjusting for these systematic biases and to evaluate the significance of differences between treatments.

There are two statistical approaches frequently used in mixed model analysis for microarray data. The first approach is to fit multiple gene-specific models to data that are normalized for systematic biases (Wolfinger et al. 2001; Gibson and Wolfinger 2004). This approach is based on assuming that the biases are independent from the gene effects. If this assumption is untenable, then a second approach fits a single model that combines both the systematic biases and the gene effects (Kerr, Martin, and Churchill 2000; Churchill 2002; Littell et al. 2006). When the number of genes is very large (several hundreds to tens of thousands), this is an analysis for which the sparse matrix approach that the LMIXED procedure implements is well suited.

The following SAS statements simulate a microarray experiment with a so-called loop design structure, which is commonly used in such studies. There are 500 genes, each gene occurs in six arrays, and each array has two dyes. You are interested in the best five and worst five genes based on the best linear unbiased predicted values (BLUP) of these 500 genes.

%let narray  = 6;
%let ndye    = 2;
%let nrow    = 4;
%let ngene   = 500;
%let ntrt    = 6;
%let npin    = 4;
%let ndip    = 4;
%let no      = %eval(&ndye*&nrow*&ngene);
%let tno     = %eval(&narray*&no);

data mycas.microarray;
   keep Gene MArray Dye Trt Pin Dip log2i;
   array PinDist{&tno};
   array DipDist{&tno};
   array GeneDist{&tno};

   array ArrayEffect{&narray};
   array ArrayGeneEffect{%eval(&narray*&ngene)};
   array ArrayDipEffect{%eval(&narray*&ndip)};
   array ArrayPinEffect{%eval(&narray*&npin)};

   do i = 1 to &tno;
      PinDist{i}  = 1 + int(&npin*ranuni(12345));
      DipDist{i}  = 1 + int(&ndip*ranuni(12345));
      GeneDist{i} = 1 + int(&ngene*ranuni(12345));
   end;

   igene = 0;
   idip = 0;
   ipin = 0;
   do i = 1 to &narray;
      ArrayEffect{i} = sqrt(0.014)*rannor(12345);
      do j = 1 to &ngene;
         igene = igene+1;
         ArrayGeneEffect{igene} = sqrt(0.0017)*rannor(12345);
      end;
      do j = 1 to &ndip;
         idip = idip + 1;
         ArrayDipEffect{idip} = sqrt(0.0033)*rannor(12345);
      end;
      do j = 1 to &npin;
         ipin = ipin + 1;
         ArrayPinEffect{ipin} = sqrt(0.037)*rannor(12345);
      end;
   end;

   i = 0;
   do MArray = 1 to &narray;
      do Dye = 1 to &ndye;
         do Row = 1 to &nrow;
            do k = 1 to &ngene;
               if MArray=1 and Dye = 1 then do;
                  Trt = 0;
                  trtc = 0;
               end;
               else do;
                  if trtc >= &no then trtc = 0;
                  if trtc = 0 then do;
                     Trt = Trt + 1;
                     if Trt >= &ntrt then do;
                        Trt = 0;
                        trtc = 0;
                     end;
                  end;
                  trtc = trtc + 1;
               end;
               i = i + 1;
               Pin = PinDist{i};
               Dip = DipDist{i};
               Gene = GeneDist{i};
               a   = ArrayEffect{MArray};
               ag  = ArrayGeneEffect{(MArray-1)*&ngene+Gene};
               ad  = ArrayDipEffect{(MArray-1)*&ndip+Dip};
               ap  = ArrayPinEffect{(MArray-1)*&npin+Pin};
               log2i   = 1 +
                        + Dye
                        + Trt
                        + Gene/1000.0
                        + Dye*Gene/1000.0
                        + Trt*Gene/1000.0
                        + Pin
                        + a
                        + ag
                        + ad
                        + ap
                        + sqrt(0.02)*rannor(12345);
               output;
            end;
         end;
      end;
   end;
run;

A linear mixed model for fitting the log intensity data from such a design is described by Littell et al. (2006) as follows:

You can use the LMIXED procedure with the following statements to fit this model:

ods select Dimensions covParms SolutionR;
proc lmixed data=mycas.microarray dmmethod=sparse;
   class marray dye trt gene pin dip;
   model log2i = dye trt gene dye*gene trt*gene pin;
   random int gene dip pin/subject=marray s;
   ods output solutionr=BLUPs;
run;

By default, because the program has a RANDOM statement with a SUBJECT= effect, PROC LMIXED uses the dense method (see the section Common Subject Effect) in the restricted maximum likelihood (REML) estimation. However, because the design matrices are excessively large with having 4,513 columns and having 3,054 columns (see table Output 10.1.1), the program can potentially exceed your computer’s memory capacity. Here the sparse method (DMMETHOD=SPARSE) is used to improve the handling of these large and sparse design matrices.

Output 10.1.2 shows the "Covariance Parameter Estimation" table, which lists the estimations of covariances across each microarray.

The following program ranks genes on the basis of their BLUP estimates.

ods output off;
proc sort data=BLUPs(firstobs=2 obs=501);
    by descending Estimate;
run;

data BLUPs; set BLUPs;
   Rank = _N_;
run;
ods output on;

proc print data=BLUPs;
   where ((Rank <= 5) | (Rank >= 496));
   var Gene Estimate;
run;

The best five genes and worst five genes with their BLUP estimates are shown in Output 10.1.3.