Complex non-linear systems biology models comprise relevant knowledge on processes of pharmacological interest. They are, however, too complex to be used in inferential settings, for example, to allow for the estimation of patient-specific parameters for individual dose optimisation. Thus, there is a need for simple models with interpretable components to infer the drug effect in a clinical setting. In particular, it is essential to accurately quantify and simulate the interindividual variability in the drug response in order to account for covariates like body weight, age and genetic disposition. To this end, non-linear model order reduction and simplification methods can be used if they maintain model interpretability during reduction and consider an entire population rather than just a single reference individual. We present a sample-based approach for robust model order reduction and propose two improvements for efficiency. In particular, we introduce a new sampling method to generate the virtual population based on transformed latin hypercube sampling. Thereby, the sample is stratified in the relevant parameter-space directions, which are identified using empirical observability Gramians. We illustrate our approach in application to a blood coagulation pathway model, where we reduce the complexity from a 62-dimensional highly non-linear to a six-dimensional and a nine-dimensional system of ordinary differential equations for two scenarios, respectively.