# 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...!

-
Anyone? I'm not sure how I can do a practice example? Should I make some dummy data, and see how you guys would control for the other variables? –  JackOfAll Oct 21 '11 at 21:33
A large proportion of the questions tagged "regression" (of which there are currently about 600) will provide explicit examples. –  whuber Oct 21 '11 at 21:39
Epi & Bernd, Thanks so much for trying to answer this. Unfortunately, these answers are a big leap from my question, and are over my head. Maybe it' b/c I don't have experience with R, and just a basic Statistics 101 foundation. Just as feedback to your teaching, once you abstracted away from BMI, age, impatience, etc to "covariate" et al, you totally lost me. Auto-generating pseudo-data also was not helpful in clarifying the concepts. In fact, it made it worse. It's hard to learn on dummy data with no inherent meaning, unless you already know the principle being explained (ie: Teacher knows i –  JackOfAll Oct 23 '11 at 23:20
This is not a complete answer or anything, but I think it's worthwhile to read "Let's Put Garbage-Can Regressions and Garbage-Can Probits Where They Belong" by Chris Achen. (PDF link: http://qssi.psu.edu/files/Achen_GarbageCan.pdf) This applies to both Bayesian and Frequentist approaches equally. Just throwing terms into your set-up is not sufficient to "control" for effects, but sadly this is what passes for control in a lot of the literature. –  EMS Oct 5 '12 at 19:41

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?

-
Thanks. Anyone know a simple example regression based example online or in a textbook that I can work through? –  JackOfAll Oct 20 '11 at 21:47
@JackOfAll There are likely hundreds of such examples - what areas/types of questions are you interested in, and what software packages can you use? –  Fomite Oct 20 '11 at 21:49
Well, any academic/contrived example is fine by me. I have Excel, which can do a multi-variable regression, correct? Or do I need something like R to do this? –  JackOfAll Oct 20 '11 at 22:25
@JackOfAll I'll see if I can find/rig up a contrived example for you :) Honestly, I don't know for Excel, but I'll include R code as well. The example won't exactly be tricky, so you should be fine. –  Fomite Oct 21 '11 at 16:24
+1 For answering this without the negativity that I would use. :) In typical parlance, controlling for other variables means the authors threw them into the regression. It doesn't really mean what they think it means if they have not validated that the variables are relatively independent and that the entire model structure (usually some kind of GLM) is well-founded. In short, my view is that whenever someone uses this phrase, it means they have very little clue about statistics, and one should re-calculate the results using the stratification method you offered. –  Iterator Oct 22 '11 at 12:49

1. Introduction

I like @EpiGrad's answer (+1) but let me take a different perspective. In the following I am referring to this PDF document: "Multiple Regression Analysis: Estimation", which has a section on "A 'Partialling Out' Interpretation of Multiple Regression" (p. 83f.). Unfortunately, I have no idea who is the author of this chapter and I will refer to it as REGCHAPTER. A similar explanation can be found in Kohler/Kreuter (2009) "Data Analysis Using Stata", chapter 8.2.3 "What does 'under control' mean?".

I will use @EpiGrad's example to explain this approach. R code and results can be found in the Appendix.

It also should be noted that "controling for other variables" does only make sense when the explanatory variables are moderately correlated (collinearity). In the aforementioned example, the Product-Moment correlation between exposure and covariate is 0.50, i.e.

> cor(covariate, exposure)
[1] 0.5036915


2. Residuals

I assume that you have a basic understanding of the concept of residuals in regression analysis. Here is the Wikipedia explanation: " If one runs a regression on some data, then the deviations of the dependent variable observations from the fitted function are the residuals".

3. What does 'under control' mean?

Controlling for the variable covariate, the effect (regression weight) of exposure on outcome can be described as follows (I am sloppy and skip most indices and all hats, please refer to the above mentioned text for a precise description):

$\beta_1=(\sum resid_{i1} \cdot y_i)/(\sum resid^2_{i1})$

$resid_{i1}$ are the residuals when we regress exposure on covariate, i.e.

$exposure = const. + \beta_{covariate} \cdot covariate + resid$.

The "residuals [..] are the part of $x_{i1}$ that is uncorrelated with $x_{i2}$. [...] Thus, $\hat{\beta}_1$ measures the sample relationship between $y$ and $x_1$ after $x_2$ has been partialled out" (REGCHAPTER 84). "Partialled out" means "controlled for".

I will demonstrate this idea using @EpiGrad's example data. First, I will regress exposure on covariate. Since I am only interested in the residualslmEC.resid, I omit the output.

summary(lmEC <- lm(exposure ~ covariate))
lmEC.resid <- residuals(lmEC)


The next step is to regress outcome on these residuals (lmEC.resid):

[output omitted]

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.45074    0.02058 119.095  < 2e-16 ***
lmEC.resid   0.50000    0.07612   6.569 2.45e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

[output omitted]


As you can see, the regression weight for lmEC.resid (see column Estimate, $\beta_{lmEC.resid}=0.50$) in this simple regression is equal to the multiple regression weight for covariate, which also is $0.50$ (see @EpiGrad's answer or the R output below).

Appendix

R Code

set.seed(1)
covariate <- sample(0:1, 100, replace=TRUE)
exposure <- runif(100,0,1)+(0.3*covariate)
outcome <- 2.0+(0.5*exposure)+(0.25*covariate)

## Simple regression analysis
summary(lm(outcome ~ exposure))

## Multiple regression analysis
summary(lm(outcome ~ exposure + covariate))

## Correlation between covariate and exposure
cor(covariate, exposure)

## "Partialling-out" approach
## Regress exposure on covariate
summary(lmEC <- lm(exposure ~ covariate))
## Save residuals
lmEC.resid <- residuals(lmEC)
## Regress outcome on residuals
summary(lm(outcome ~ lmEC.resid))

## Check formula
sum(lmEC.resid*outcome)/(sum(lmEC.resid^2))


R Output

> set.seed(1)
> covariate <- sample(0:1, 100, replace=TRUE)
> exposure <- runif(100,0,1)+(0.3*covariate)
> outcome <- 2.0+(0.5*exposure)+(0.25*covariate)
>
> ## Simple regression analysis
> summary(lm(outcome ~ exposure))

Call:
lm(formula = outcome ~ exposure)

Residuals:
Min        1Q    Median        3Q       Max
-0.183265 -0.090531  0.001628  0.085434  0.187535

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.98702    0.02549   77.96   <2e-16 ***
exposure     0.70103    0.03483   20.13   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.109 on 98 degrees of freedom
Multiple R-squared: 0.8052,     Adjusted R-squared: 0.8032
F-statistic: 405.1 on 1 and 98 DF,  p-value: < 2.2e-16

>
> ## Multiple regression analysis
> summary(lm(outcome ~ exposure + covariate))

Call:
lm(formula = outcome ~ exposure + covariate)

Residuals:
Min         1Q     Median         3Q        Max
-7.765e-16 -7.450e-18  4.630e-18  1.553e-17  4.895e-16

Coefficients:
Estimate Std. Error   t value Pr(>|t|)
(Intercept) 2.000e+00  2.221e-17 9.006e+16   <2e-16 ***
exposure    5.000e-01  3.508e-17 1.425e+16   <2e-16 ***
covariate   2.500e-01  2.198e-17 1.138e+16   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.485e-17 on 97 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:     1
F-statistic: 3.322e+32 on 2 and 97 DF,  p-value: < 2.2e-16

>
> ## Correlation between covariate and exposure
> cor(covariate, exposure)
[1] 0.5036915
>
> ## "Partialling-out" approach
> ## Regress exposure on covariate
> summary(lmEC <- lm(exposure ~ covariate))

Call:
lm(formula = exposure ~ covariate)

Residuals:
Min       1Q   Median       3Q      Max
-0.49695 -0.24113  0.00857  0.21629  0.46715

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.51003    0.03787  13.468  < 2e-16 ***
covariate    0.31550    0.05466   5.772  9.2e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2731 on 98 degrees of freedom
Multiple R-squared: 0.2537,     Adjusted R-squared: 0.2461
F-statistic: 33.32 on 1 and 98 DF,  p-value: 9.198e-08

> ## Save residuals
> lmEC.resid <- residuals(lmEC)
> ## Regress outcome on residuals
> summary(lm(outcome ~ lmEC.resid))

Call:
lm(formula = outcome ~ lmEC.resid)

Residuals:
Min      1Q  Median      3Q     Max
-0.1957 -0.1957 -0.1957  0.2120  0.2120

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.45074    0.02058 119.095  < 2e-16 ***
lmEC.resid   0.50000    0.07612   6.569 2.45e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2058 on 98 degrees of freedom
Multiple R-squared: 0.3057,     Adjusted R-squared: 0.2986
F-statistic: 43.15 on 1 and 98 DF,  p-value: 2.45e-09

>
> ## Check formula
> sum(lmEC.resid*outcome)/(sum(lmEC.resid^2))
[1] 0.5
>

-
+1 for being an excellent treatment of the question. –  Fomite Oct 22 '11 at 20:29
That chapter looks like Baby Wooldridge (aka Introductory Econometrics: A Modern Approach by Jeffrey M. Wooldridge) –  Dimitriy V. Masterov Oct 6 '12 at 0:29