How exactly does one “control for other variables”?

Question

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

Answer 1

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?

Answer 2

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
>