Much of the phenotypic variation observed between even closely related species may be driven by differences in gene expression levels. The current availability of reliable techniques like RNA-Seq, which can quantify expression levels across species, has enabled comparative studies. Ornstein-Uhlenbeck (OU) processes have been proposed to model gene expression evolution as they model both random drift and stabilizing selection and can be extended to model changes in selection regimes. The OU models provide a statistical framework that allows comparisons of specific hypotheses of selective regimes, including random drift, constrained drift, and expression level shifts. In this way, inferences may be made about the mode of selection acting on the expression level of a gene. We augment this model to include within-species expression variance, allowing for modeling of nonevolutionary expression variance that could be caused by individual genetic, environmental, or technical variation. Through simulations, we explore the reliability of parameter estimates and the extent to which different selective regimes can be distinguished using phylogenies of varying size using both the typical OU model and our extended model. We find that if individual variation is not accounted for, nonevolutionary expression variation is often mistaken for strong stabilizing selection. The methods presented in this article are increasingly relevant as comparative expression data becomes more available and researchers turn to expression as a primary evolving phenotype.