The HPLMIXED Procedure

Example 8.1 Computing BLUPs for a Large Number of Subjects

Suppose you are using health measurements on patients treated by each medical center to monitor the performance of those centers. Different measurements within each patient are correlated, and there is enough data to fit the parameters of an unstructured covariance model for this correlation. In fact, long experience with historical data provides you with values for the covariance model that are essentially known, and the task is to apply these known values in order to compute best linear unbiased predictors (BLUPs) of the random effect of medical center. You can use these BLUPs to determine the best and worst performing medical centers, adjusting for other factors, on a weekly basis. Another reason why you want to do this with fixed values for the covariance parameters is to make the week-to-week BLUPs more comparable.

Although you cannot use the REPEATED statement in PROC HPLMIXED to fit models in this release, you can use it to compute BLUPs for such models with known values of the variance parameters. For illustration, the following statements create a simulated data set of a given week’s worth of patient health measurements across 100 different medical centers. Measurements at three different times are simulated for each patient, and each center has about 50 patients. The simulated model includes a fixed gender effect, a random effect due to center, and covariance between different measurements on the same patient.

   %let NCenter  = 100;
   %let NPatient = %eval(&NCenter*50);
   %let NTime    = 3;
   %let SigmaC   = 2.0;
   %let SigmaP   = 4.0;
   %let SigmaE   = 8.0;
   %let Seed     = 12345;

   data WeekSim;
      keep Gender Center Patient Time Measurement;
      array PGender{&NPatient};
      array PCenter{&NPatient};
      array PEffect{&NPatient};
      array CEffect{&NCenter};
      array GEffect{2};

      do Center = 1 to &NCenter;
         CEffect{Center} = sqrt(&SigmaC)*rannor(&Seed);
         end;

      GEffect{1} = 10*ranuni(&Seed);
      GEffect{2} = 10*ranuni(&Seed);

      do Patient = 1 to &NPatient;
         PGender{Patient} = 1 + int(2       *ranuni(&Seed));
         PCenter{Patient} = 1 + int(&NCenter*ranuni(&Seed));
         PEffect{Patient} = sqrt(&SigmaP)*rannor(&Seed);
         end;

      do Patient = 1 to &NPatient;
         Gender = PGender{Patient};
         Center = PCenter{Patient};
         Mean = 1 + GEffect{Gender} + CEffect{Center} + PEffect{Patient};
         do Time = 1 to &nTime;
            Measurement = Mean + sqrt(&SigmaE)*rannor(&Seed);
         output;
         end;
      end;
   run;

Suppose that the known values for the covariance parameters are

\begin{eqnarray*}  \mbox{Var(Center)} &  = &  1.7564 \\ \mbox{Cov(Patient)} &  = &  \left[ \begin{array}{ccc} 11.4555 &  3.6883 &  4.5951 \\ 3.6883 &  11.2071 &  3.6311 \\ 4.5951 &  3.6311 &  12.1050 \end{array} \right] \end{eqnarray*}

Incidentally, these are not precisely the same estimates you would get if you estimated these parameters based on the preceding data (for example, with the HPLMIXED procedure).

The following statements use PROC HPLMIXED to compute the BLUPs for the random medical center effect. Instead of simply displaying them (as PROC HPMIXED does), PROC HPLMIXED sorts them and displays the five highest and lowest values. To run these statements successfully, you need to set the macro variables GRIDHOST and GRIDINSTALLLOC to resolve to appropriate values, or you can replace the references to macro variables with appropriate values.

ods listing close;
proc hplmixed data=WeekSim blup;
   performance host="&GRIDHOST" install="&GRIDINSTALLLOC" nodes=20;
   class Gender Center Patient Time;
   model Measurement = Gender;
   random   Center / s;
   repeated Time   / sub=Patient type=un;
   parms   1.7564
          11.4555
           3.6883 11.2071
           4.5951  3.6311 12.1050;
   ods output SolutionR=BLUPs;
run;
ods listing;

proc sort data=BLUPs;
   by Estimate;
run;

data BLUPs; set BLUPs;
   Rank = _N_;
run;
proc print data=BLUPs;
   where ((Rank <= 5) | (Rank >= 96));
   var Center Estimate;
run;

Three parts of the PROC HPLMIXED syntax are required in order to compute BLUPs for this model: the BLUP option in the HPLMIXED statement, the REPEATED statement, and the PARMS statement with fixed values for all parameters. The resulting values of the best and worst performing medical centers for this week are shown in Output 8.1.1. Apparently, for this week’s data, medical center 54 had the most decreasing effect, and medical center 48 the most increasing effect, on patient measurements overall.

Output 8.1.1: Highest and Lowest Medical Center BLUPs

Obs Center Estimate
1 54 -2.9369
2 7 -2.4614
3 50 -2.2467
4 51 -2.2281
5 93 -2.1644
96 26 2.1603
97 99 2.2718
98 44 2.4222
99 60 2.6089
100 48 2.6443