Variance components models, also known as mixed effects model, are a central theme in statistics. When there are a large number of variance components, one wants to select a subset of those that are associated with response. Existing methods are limited to finding random components at individual level or within one variance component. We propose selection of variance components based on a penalized log-likelihood with adaptive penalty. This is solved by a majorization-minimization (MM) algorithm, which is simple, numerically stable, and globally convergent. Performance of the proposed methodology is evaluated empirically through simulation studies and real data analysis. In theory, we establish a non-asymptotic error bound for the iterates from the algorithm and characterize the region in which the MM iterates converge to a global optimum of the population likelihood. This result provides a theoretical guideline in terminating MM iterations.