Adverse pressure gradients, separation and other forms of non-equilibrium flows are often encountered in flows of interest. In these type of flows, the Boussinesq hypothesis does not hold and often leads to erroneous predictions by eddy viscosity models. In an attempt to capture these non-equilibrium effects, lag parameter models introduce a lag parameter, which is derived from an elliptic blending Reynolds stress model. A novel deterministic machine learning algorithm, referred to as Sparse Regression of Turbulent Stress Anisotropy (SpaRTA), has been used with the objective of developing a data-driven turbulence model based on the elliptic blending k − ω lag parameter model and evaluating its performance in terms of generalizability, interpretability and its ability to infer the quantities of interest. Corrective terms are introduced and computed directly from high-fidelity data to account for the model-form error by using the k-corrective-frozen-RANS approach. SpaRTA is then used to infer algebraic stress models for these corrective terms using a Galilean invariant integrity basis. It was shown that the k-corrective-frozen-RANS framework has the ability of representing the mean flow features by propagating the corrective terms through a CFD model in OpenFOAM and comparing its result to high-fidelity data. Cross-validation is used to test the performance of the models on unseen data using three flow cases that involve separation, namely periodic hills (Re=10595), converging-diverging channel (Re=12600) and curved backward-facing step (Re=13700). In order to assess the impact of the additional transport equation of the lag parameter, the same data-driven approach was applied to the conventional two-equation k − ω model. Utilizing an additional transport equation for the lag parameter in this data-driven approach did not result in any significant improvements in terms of predictive capability or generalizability, as both data-driven approaches showed a similar performance, although the data-driven k − ω models were more numerically stable. It was found that corrective terms formulated using a reduced integrity basis yields data-driven models that have a similar predictive capability compared to models that used the full integrity basis to construct the corrective terms. A significant portion of the resulting data-driven models showed an improvement in predictive capability over the standard (non-data-driven) k − ω model. Furthermore, most of the models were able to generalize their predictions to two-dimensional flow cases that had different complexity and showed a significant improvement over the baseline k − ω model.