Rao Blackwellized Gibbs

Introduction

The procedure described in this section is still experimental. Users are requested to share their experience with the developers.

Models which are too large to fit with a standard model may be fitted using a specialized Gibbs sampling approach which cycles over submodels. As an introductory example, an analysis of the model y ~ mu A !r B on 500,000 records when A and B had 25,000 levels each was taking around 20 minutes per iteration and needs 10-20 iterations to converge. Using the code

 !RBGIBBS 20
 y ~ mu !SM 1 A !r !SM 2 B

performed 220 Gibbs cycles in under a minute. Depending on the user requirement, the user could use the variance components obtained from the Gibbs run as start values for a standard run which would then require just 3 iterations.

The !RBGIBBS n qualifier sets the number of BURNIN cycles. !MAXIT defaults to n*11. So in the example, parameters and effects reported were the average of 200 cycles (after discarding the first 20).

The essence of the approach is to partition the model into submodels which individually are easy to fit. So in the example, each cycle consists of two model fits, y-Zβ ~ μ +Xα and y-Xα ~ μ +Zβ. Sampling variation is added to the effects before adjusting the response to remove bias.

A larger example analysing 1 million records from 260,000 animals took 2 hours for the first iteration when fitted conventially and 1 minute per cycle when fitted as 3 parts using:

 !RBG 20
 LR ~ mu CNV NP EDAD LD REB !r !sm 2 ID  !F !sm 1 ref(FP) +
        !r !sm 3 giv(PedID,1)

The average parameter values from the RBG run were 0.2527 ( ID) and 0.1977 ( giv(PedID,1)) and fitted in the standard way, the LogLs were

    1 LogL=-4954.28     S2=  1292.4       995814 df
    2 LogL=-4954.02     S2=  1292.6       995814 df

yielding updated parameter values of 0.2509 and 0.1997 respectively.

When one of the submodels contains fixed effects, an inefficiency is introduced unless care is taken. Fitting the Sparse fixed effects sweeps out the mean ( mu) but adding noise to the effects results in the mean of the adjusted data not being Zero. However, the effect mu will be flagged as aliased and not fitted in the other submodels. The best solution is to formally add a constraint to the fixed factor in the submodel, and the most efficient constraint is to set one level to Zero. E.g., in our motivating example, use ref(A,1) in the model instead of A.

This analysis uses ideas from Gibbs sampling (Thompson 1994) as described by Gilmour and Thompson (2003) building on work by Clayton and Rasbash (1999), Garcia-Cortes and Sorensen (2001) and Sorensen and Gianola (2002). Further information is available in a document available from the author.

Return to index