How exactly does one “control for other variables”? Here is the article that motivated this question: Does impatience make us fat?
I liked this article, and it nicely demonstrates the concept of “controlling for other variables” (IQ, career, income, age, etc) in order to best isolate the true relationship between just the 2 variables in question.  
Can you explain to me how you actually control for variables on a typical data set?    
E.g., if you have 2 people with the same impatience level and BMI, but different incomes, how do you treat these data?  Do you categorize them into different subgroups that do have similar income, patience, and BMI?  But, eventually there are dozens of variables to control for (IQ, career, income, age, etc)  How do you then aggregate these (potentially) 100’s of subgroups?  In fact, I have a feeling this approach is barking up the wrong tree, now that I’ve verbalized it.
Thanks for shedding any light on something I've meant to get to the bottom of for a few years now...!
 A: There is an excellent discussion so far of covariate adjustment as a means of "controlling for other variables". But I think that is only part of the story. In fact, there are many (other) design, model, and machine learning based strategies to address the impact of a number of possible confounding variables. This is a brief survey of some of the most important (non-adjustment) topics. While adjustment is the most widely used means of "controlling" for other variables, I think a good statistician should have an understanding of what it does (and doesn't do) in the context of other processes and procedures.
Matching:
Matching is a method of designing a paired analysis where observations are grouped into sets of 2 who are otherwise similar in their most important aspects. For instance, you might sample two individuals who are concordant in their education, income, professional tenure, age, marital status, (etc. etc.) but who are discordant in terms of their impatience. For binary exposures, the simple paired-t test suffices to test for a mean difference in their BMI controlling for all the matching features. If you are modeling a continuous exposure, an analogous measure would be a regression model through the origin for the differences. See Carlin 2005
$$E[Y_1 - Y_2] = \beta_0 (X_1 - X_2)$$
Weighting
Weighting is yet another univariate analysis which models the association between a continuous or binary predictor $X$ and an outcome $Y$ so that the distribution of exposure levels is homogeneous between groups. These results are typically reported as standardized such as age-standardized mortality for two countries or several hospitals. Indirect standardization calculates an expected outcome distribution from the rates obtained in a "control" or "healthy" population that is projected to the distribution of strata in the referent population. Direct standardization goes the other way. These methods are typically used for a binary outcome. Propensity score weighting accounts of the probability of a binary exposure and controls for those variables in that regard. It is similar to direct standardization for an exposure. See Rothman, Modern Epidemiology 3rd edition.
Randomization and Quasirandomization
It's a subtle point, but if you are actually able to randomize people to a certain experimental condition, then the impact of other variables is mitigated.
It's a remarkably stronger condition, because you do not even need to know what those other variables are.
In that sense, you have "controlled" for their influence. This is not possible in observational research, but it turns out that propensity score methods create a simple probabilistic measure for exposure that allows one to weight, adjust, or match participants so that they can be analyzed in the same fashion as a quasi-randomized study. See Rosenbaum, Rubin 1983.
Microsimulation
Another way of simulating data that might have been obtained from a randomized study is to perform microsimulation. Here, one can actually turn their attention to larger and more sophisticated, machine learning like models. A term which Judea Pearl has coined that I like is "Oracle Models": complex networks which are capable of generating predictions and forecast for a number of features and outcomes. It turns out one can "fold down" the information of such an oracle model to simulate outcomes in a balanced cohort of people who represent a randomized cohort, balanced in their "control variable" distribution, and using simple t-test routines to assess the magnitude and precision of possible differences. See Rutter , Zaslavsky, and Feuer 2012
Matching, weighting, and covariate adjustment in a regression model all estimate the same associations, and thus all can be claimed to be ways of "controlling" for other variables.
A: Of course some math will be involved, but it's not much: Euclid would have understood it well.  All you really need to know is how to add and rescale vectors.  Although this goes by the name of "linear algebra" nowadays, you only need to visualize it in two dimensions.  This enables us to avoid the matrix machinery of linear algebra and focus on the concepts.

A Geometric Story
In the first figure, $y$ is the sum of  $y_{\cdot 1}$ and $\alpha x_1$.  (A vector $x_1$ scaled by a numeric factor $\alpha$; Greek letters $\alpha$ (alpha), $\beta$ (beta), and $\gamma$ (gamma) will refer to such numerical scale factors.)

This figure actually began with the original vectors (shown as solid lines) $x_1$ and $y$.  The least-squares "match" of $y$ to $x_1$ is found by taking the multiple of $x_1$ that comes closest to $y$ in the plane of the figure.  That's how $\alpha$ was found.  Taking this match away from $y$ left $y_{\cdot 1}$, the residual of $y$ with respect to $x_1$.  ( The dot "$\cdot$" will consistently indicate which vectors have been "matched," "taken out," or "controlled for.")
We can match other vectors to $x_1$.  Here is a picture where $x_2$ was matched to $x_1$, expressing it as a multiple $\beta$ of $x_1$ plus its residual $x_{2\cdot 1}$:

(It does not matter that the plane containing $x_1$ and $x_2$ could differ from the plane containing $x_1$ and $y$: these two figures are obtained independently of each other.  All they are guaranteed to have in common is the vector $x_1$.)  Similarly, any number of vectors $x_3, x_4, \ldots$ can be matched to $x_1$.
Now consider the plane containing the two residuals $y_{\cdot 1}$ and $x_{2 \cdot 1}$.  I will orient the picture to make $x_{2\cdot 1}$ horizontal, just as I oriented the previous pictures to make $x_1$ horizontal, because this time $x_{2\cdot 1}$ will play the role of matcher:

Observe that in each of the three cases, the residual is perpendicular to the match.  (If it were not, we could adjust the match to get it even closer to $y$, $x_2$, or $y_{\cdot 1}$.)
The key idea is that by the time we get to the last figure, both vectors involved ($x_{2\cdot 1}$ and $y_{\cdot 1}$) are already perpendicular to $x_1$, by construction.  Thus any subsequent adjustment to $y_{\cdot 1}$ involves changes that are all perpendicular to $x_1$.  As a result, the new match $\gamma x_{2\cdot 1}$ and the new residual $y_{\cdot 12}$ remain perpendicular to $x_1$.
(If other vectors are involved, we would proceed in the same way to match their residuals $x_{3\cdot 1}, x_{4\cdot 1}, \ldots$ to $x_2$.)
There is one more important point to make.  This construction has produced a residual $y_{\cdot 12}$ which is perpendicular to both $x_1$ and $x_2$.  This means that $y_{\cdot 12}$ is also the residual in the space (three-dimensional Euclidean realm) spanned by $x_1, x_2,$ and $y$.  That is, this two-step process of matching and taking residuals must have found the location in the $x_1, x_2$ plane that is closest to $y$.  Since in this geometric description it does not matter which of $x_1$ and $x_2$ came first, we conclude that if the process had been done in the other order, starting with $x_2$ as the matcher and then using $x_1$, the result would have been the same.
(If there are additional vectors, we would continue this "take out a matcher" process until each of those vectors had had its turn to be the matcher.  In every case the operations would be the same as shown here and would always occur in a plane.)

Application to Multiple Regression
This geometric process has a direct multiple regression interpretation, because columns of numbers act exactly like geometric vectors.  They have all the properties we require of vectors (axiomatically) and therefore can be thought of and manipulated in the same way with perfect mathematical accuracy and rigor.  In a multiple regression setting with variables $X_1$, $X_2, \ldots$, and $Y$, the objective is to find a combination of $X_1$ and $X_2$ (etc) that comes closest to $Y$.  Geometrically, all such combinations of $X_1$ and $X_2$ (etc) correspond to points in the $X_1, X_2, \ldots$ space.  Fitting multiple regression coefficients is nothing more than projecting ("matching") vectors. The geometric argument has shown that

*

*Matching can be done sequentially and


*The order in which matching is done does not matter.
The process of "taking out" a matcher by replacing all other vectors by their residuals is often referred to as "controlling" for the matcher.  As we saw in the figures, once a matcher has been controlled for, all subsequent calculations make adjustments that are perpendicular to that matcher.  If you like, you may think of "controlling" as "accounting (in the least square sense) for the contribution/influence/effect/association of a matcher on all the other variables."

References
You can see all this in action with data and working code in the answer at https://stats.stackexchange.com/a/46508.  That answer might appeal more to people who prefer arithmetic over plane pictures.  (The arithmetic to adjust the coefficients as matchers are sequentially brought in is straightforward nonetheless.)  The language of matching is from Fred Mosteller and John Tukey.
A: There are many ways to control for variables.
The easiest, and one you came up with, is to stratify your data so you have sub-groups with similar characteristics - there are then methods to pool those results together to get a single "answer". This works if you have a very small number of variables you want to control for, but as you've rightly discovered, this rapidly falls apart as you split your data into smaller and smaller chunks.
A more common approach is to include the variables you want to control for in a regression model. For example, if you have a regression model that can be conceptually described as: 
BMI = Impatience + Race + Gender + Socioeconomic Status + IQ

The estimate you will get for Impatience will be the effect of Impatience within levels of the other covariates - regression allows you to essentially smooth over places where you don't have much data (the problem with the stratification approach), though this should be done with caution.
There are yet more sophisticated ways of controlling for other variables, but odds are when someone says "controlled for other variables", they mean they were included in a regression model.
Alright, you've asked for an example you can work on, to see how this goes. I'll walk you through it step by step. All you need is a copy of R installed.
First, we need some data. Cut and paste the following chunks of code into R. Keep in mind this is a contrived example I made up on the spot, but it shows the process.
covariate <- sample(0:1, 100, replace=TRUE)
exposure  <- runif(100,0,1)+(0.3*covariate)
outcome   <- 2.0+(0.5*exposure)+(0.25*covariate)

That's your data. Note that we already know the relationship between the outcome, the exposure, and the covariate - that's the point of many simulation studies (of which this is an extremely basic example. You start with a structure you know, and you make sure your method can get you the right answer.
Now then, onto the regression model. Type the following:
lm(outcome~exposure)

Did you get an Intercept = 2.0 and an exposure = 0.6766? Or something close to it, given there will be some random variation in the data? Good - this answer is wrong. We know it's wrong. Why is it wrong? We have failed to control for a variable that effects the outcome and the exposure. It's a binary variable, make it anything you please - gender, smoker/non-smoker, etc.
Now run this model:
lm(outcome~exposure+covariate)

This time you should get coefficients of Intercept = 2.00, exposure = 0.50 and a covariate of 0.25. This, as we know, is the right answer. You've controlled for other variables.
Now, what happens when we don't know if we've taken care of all of the variables that we need to (we never really do)? This is called residual confounding, and its a concern in most observational studies - that we have controlled imperfectly, and our answer, while close to right, isn't exact. Does that help more?
A: The software doesn't literally control for variables. If you're familiar with matrix notation of regression $Y=X\beta+\varepsilon$, then you may remember that least squares solution is $b=(X^TX)^{-1}X^TY$. So, the software evaluates this expression numerically using computational linear algebra methods.
