Estimating the mean of a variable by controlling for age etc.: Why using median values? just read this thread here (How do I estimate the mean of a variable for different groups by controlling for age, gender, education...?) and found it very helpful! First of all, thanks for that to the authors. 
I have a question that is related to the question in the thread linked above: 
In a paper I am reading it is said: 
"The means for each Household-quintile were estimated from the regression intercept and coefficients controlling for all socio-demographic characteristics (i.e. median age (5 years), female, median birth order (2nd), mothers with primary education, crowding (3.6–
5 people ⁄ room), rural setting, and the most represented province (Toliary)), and adjusting for clustering at the household and
community level."
Having read the post above, I do understand that they were doing a regression and using the intercept and the coefficient as the estimated means. 
BUT 1) why are they using median age and median birth order and only the most represented province instead of a set of province indicator variables as in the regressions before?
AND 2) how can I cluster at the household and the community level at the same time? 
I appreciate your help very much since I have been thinking about those questions quite a while already. 
THANKS!
PS: I am using Stata in case this information is needed. 
 A: *

*From a quick skim of the paper, it looks like those are merely the covariate values used to predict the scores in the graphs after the model was estimated. The authors allow HH wealth and maternal education to vary, keeping the other covariates fixed at the median/modal values. These numbers seem to correspond to median/modal values from Table 1 (summary stats). If you allow province to vary, you would get a set of very similar graphs with the bars all shifted up or down. Most readers care about the wealth and education effects, so the authors focus on how scores change when those two variables increase, for a kind of representative student (a five-year-old female, etc.). This also allows them to put the predicted scores and the coefficients in the same graph with two axes.


You can do something like this in Stata with (setting mpg to the median and foreign dummy to the base of zero and allowing weight quintile to vary):
sysuse auto, clear
xtile wq = weight, nq(5)
lab var wq "Weight Quintile"

reg price i.foreign i.wq c.mpg
margins, dydx(wq) post
estimates store betas

reg price i.foreign i.wq c.mpg
margins wq, at((median) mpg (base) foreign) post
estimates store yhats

capture ssc install coefplot
coefplot yhats betas, baselevels xlab(#10) xline(0)

This yields predictions, which are just $$\hat y = \hat \beta_{q}+(4999.846 -79.08586\cdot 20)=\hat \beta_{q}+3418.1288:$$

You can get coefplot to give you two x-axes and bars, but I think it makes for a worse chart, so I did not do it.


*For two-way clustering, you can use cluster2, cgmreg, or
(xt)ivreg2 from SSC. The xt* is for panel data models, the others for OLS.

