Assume we are having observational data. There is a continuous outcome variable and a grouping variable. I would like to plot the distribution of the outcome variable by group but I want to adjust them for potential confounders, so that the final distributions are net of these other variables. In other words, I want to plot the distributions of the three groups if the adjusting variables were highly similar in all groups. What is the best way to do this? My first idea was an OLS regression and then predict values and plot them, like:

Code:
reg outcome i.group conf1 conf2 conf3
predict adjusted
histogram adjusted if group == 1 ...
Is this sufficient? Would I need interaction terms in the regression model? I thought about balancing approaches as well but this will only work for 2 groups, right?