We develop a new approximative estimation method for conditional Shapley values obtained using a linear regression model. We develop a new estimation method and outperform existing methodology and implementations. Compared to the sequential method in the shapr-package (i.e fit one and one model), our method runs in minutes and not in hours. Compared to the iterative method in the shapr-package, we obtain better estimates in less than or almost the same amount of time. When the number of covariates becomes too large, one can still fit thousands of regression models at once using our method. We focus on a linear regression model, but one can easily extend the method to accommodate several types of splines that can be estimated using multivariate linear regression due to linearity in the parameters.