Holding other predictors constant via simulation in R Imagine predicting salary of some professors from their years of experience (time) controlling for/holding constant their number of publications (pubs).

Question: Is the following regarding the meaning of holding constant their number of pubs correct, and demonstrable via simulation in R?

Imagine we had countless professors, then take a sample of them with the exact same number of pubs (e.g., $1$).

*

*Fit a regression model with only time as predictor, get the regression coef of time.

*Take another sample with pubs of $2$, Fit the regression model again, get the regression coef of time.

*Keep changing pubs to $3, 4,…$ and each time get the regression coef of time.

At the end, average of our regression coefs of time will be a partial regression coefficient that has controlled for the pubs of professors while predicting salary from time.

p.s. Is controlling for a predictor similar to integrating it out?

 A: Yes, if the model is correctly specified.
Suppose your data is generated by
$$
y = \beta_1 x_1 + \beta_2 x_2 + \epsilon, \mbox{ where } E[\epsilon|x_1, x_2] = 0,
$$
i.e.
$$
E[y|x_1, x_2] = \beta_1 x_1 + \beta_2 x_2.
$$
Suppose $x_1$ is the predictor of interest and $x_2$ is control. Conditioning on the control $x_2$ gives
$$
E[y|x_2] = \beta_1 E[x_1|x_2] + \beta_2 x_2. \quad (*)
$$
The empirical counterpart of $(*)$ is the regression you're suggesting---regress $y$ on $x_1$ (with intercept) for a given value of $x_2$.
Note that for any given value of $x_2$, this regression conditional on $x_2$ is already a unbiased estimator of $\beta_1$.
Averaging over $x_2$ makes estimate less noisy. The assumption $E[\epsilon|x_1, x_2] = 0$ implies samples are uncorrelated across $x_2$. Therefore averaging over $x_2$ gives a smaller standard error.
Comment
The statement "the regression conditional on $x_2$ is a unbiased estimator of $\beta_1$" is contingent upon correct specification---correct functional form/no omitted variables/etc. In a real data set, you would have to willing to believe/claim true functional form is linear/no controls are omitted/etc.
If the true population regression function is not linear but $E[\epsilon|x_1, x_2] = 0$ still holds, I would expect averaging the OLS coefficient for $x_1$ from the regression conditional on $x_2$, call it $\hat{\beta}_1|x_2$, over $x_2$ to be close to the OLS coefficient $\hat{\beta}_1$.
