Bayesian analysis via Gibbs sampling, restricted maximum likelihood (REML), and Method R were used to estimate variance components for several models of simulated data. Four simulated data sets that included direct genetic effects and different combinations of maternal, permanent environmental, and dominance effects were used. Parents were selected randomly, on phenotype across or within contemporary groups, or on BLUP of genetic value. Estimates by Bayesian analysis and REML were always empirically unbiased in large data sets. Estimates by Method R were biased only with phenotypic selection across contemporary groups; estimates of the additive variance were biased upward, and all the other estimates were biased downward. No empirical bias was observed for Method R under selection within contemporary groups or in data without contemporary group effects. The bias of Method R estimates in small data sets was evaluated using a simple direct additive model. Method R gave biased estimates in small data sets in all types of selection except BLUP. In populations where the selection is based on BLUP of genetic value or where phenotypic selection is practiced mostly within contemporary groups, estimates by Method R are likely to be unbiased. In this case, Method R is an alternative to single-trait REML and Bayesian analysis for analyses of large data sets when the other methods are too expensive to apply.