We propose a multithreshold change plane regression model which naturally partitions the observed subjects into subgroups with different covariate effects. The underlying grouping variable is a linear function of observed covariates and thus multiple thresholds produce change planes in the covariate space. We contribute a novel two‐stage estimation approach to determine the number of subgroups, the location of thresholds, and all other regression parameters. In the first stage we adopt a group selection principle to consistently identify the number of subgroups, while in the second stage change point locations and model parameter estimates are refined by a penalized induced smoothing technique. Our procedure allows sparse solutions for relatively moderate‐ or high‐dimensional covariates. We further establish the asymptotic properties of our proposed estimators under appropriate technical conditions. We evaluate the performance of the proposed methods by simulation studies and provide illustrations using two medical data examples. Our proposal for subgroup identification may lead to an immediate application in personalized medicine.