Hello, I am doing some regression analysis where I want to graph two regressions, one for con19 against retired and one for con19 against ret_eng, where the first regression line should show con19 = _cons + retired*x (replacing _cons and retired with their coefficients shown in the first regression output below) and the second regression line should show con19 = (_cons+england) + (retired+ret_eng)*x (replacing _cons, retired and ret_eng with their coefficients shown in the first regression output below).
My problem is that for the regressions I want to graph, the first one ends up illustrating con19 against retired, as if it takes the coefficients used in the second regression output, and the second one illustrates con19 against ret_eng, as if it takes the coefficients used in the third regression output. Neither regressions take into account the impact of the other variables, the code I am using to plot these graphs is:
Code:
twoway lfit con19 retired || lfit con19 ret_eng
I don't know what I am meant to code to tell STATA to plot a graph based on the output of the first regression, I know at the moment I am essentially telling STATA to regress con19 against retired and just plot that, and then regress con19 against ret_eng and just plot that, whereas what I actually want is to tell STATA to regress con19 against retired, england and ret_eng and then fit a line of con19 against retired, and fit another line of con19 against ret_eng, using the correct coefficients of the first regression.

So, my question is, how to I plot a line-of-best-fit graph where I have multiple independent variables on STATA?

Code:
. regress con19 retired england ret_eng

Source | SS df MS Number of obs = 632
-------------+---------------------------------- F(3, 628) = 193.42
Model | 85582.293 3 28527.431 Prob > F = 0.0000
Residual | 92621.868 628 147.487051 R-squared = 0.4802
-------------+---------------------------------- Adj R-squared = 0.4778
Total | 178204.161 631 282.415469 Root MSE = 12.144

------------------------------------------------------------------------------
con19 | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
retired | 2.225232 .4406403 5.05 0.000 1.359925 3.090539
england | 16.28971 7.233304 2.25 0.025 2.085322 30.49411
ret_eng | .3576843 .4605856 0.78 0.438 -.5467902 1.262159
_cons | -5.867193 6.963177 -0.84 0.400 -19.54112 7.806736
------------------------------------------------------------------------------

. regress con19 retired

Source | SS df MS Number of obs = 632
-------------+---------------------------------- F(1, 630) = 222.97
Model | 46583.2751 1 46583.2751 Prob > F = 0.0000
Residual | 131620.886 630 208.922041 R-squared = 0.2614
-------------+---------------------------------- Adj R-squared = 0.2602
Total | 178204.161 631 282.415469 Root MSE = 14.454

------------------------------------------------------------------------------
con19 | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
retired | 2.256395 .1511097 14.93 0.000 1.959656 2.553135
_cons | 11.66918 2.236417 5.22 0.000 7.27745 16.06092
------------------------------------------------------------------------------

. regress con19 ret_eng

Source | SS df MS Number of obs = 632
-------------+---------------------------------- F(1, 630) = 432.25
Model | 72514.5752 1 72514.5752 Prob > F = 0.0000
Residual | 105689.586 630 167.761247 R-squared = 0.4069
-------------+---------------------------------- Adj R-squared = 0.4060
Total | 178204.161 631 282.415469 Root MSE = 12.952

------------------------------------------------------------------------------
con19 | Coef. Std. Err. t P>|t| [95% Conf. Interval]
-------------+----------------------------------------------------------------
ret_eng | 1.712313 .08236 20.79 0.000 1.55058 1.874047
_cons | 23.62391 1.104736 21.38 0.000 21.4545 25.79332
------------------------------------------------------------------------------

.