Statistical approaches for testing hypotheses of heterogeneity in fitness functions are needed to accommodate studies of phenotypic selection with repeated sampling across study units, populations or years. In this study, we tested directly for among-population variation in complex fitness functions and demonstrate a new approach for locating the region of the trait distribution where variation in fitness and traits is greatest. We modelled heterogeneity in fitness functions among populations by treating regression coefficients of fitness on traits as random variates. We applied random regression using two model specifications, (i) spline-based curve and (ii) stepwise, to a 2-year study of selection among 16 populations of the gall wasp, Belonocnema treatae. Log-likelihood ratio tests of variance components and 10-fold cross-validation were used to assess the evidence that selection varied among populations. Ten-fold cross-validation prediction error sums of squares (PSS) indicated that spline-based fitness functions were population specific and that the strength of evidence for heterogeneity in selection differed between years. Hypothesis testing of variance components from both models was consistent with the PSS results. Both the stepwise model and the local prediction error estimates of spline-based fitness functions identified the region(s) of the phenotype distribution harbouring the greatest heterogeneity among populations. The adopted framework advances our understanding of phenotypic selection in natural populations by extending the analysis of spline-based fitness functions to testing for heterogeneity among study units and isolating the regions of the phenotypic distribution where this variation is most pronounced.