3
$\begingroup$

I'll give an example.

In a simple correlation, blood pressure is associated with body fat.

However, I want to adjust %body fat for age and sex. So, in SAS, PROC CORR I use age and sex as partial variables and the association remains significant between blood pressure and body fat(adjusted).

How can I illustrate this in a graph (i.e. blood pressure vs %bodyfat (adjusted))?

My first though was to use a linear regression (PROC REG/PROC GLM) with %bodyfat as the dependent variable and age, sex as the independent/quantitative variables and then output the residuals. Add the mean of %bodyfat to the residuals and then plot this as my y values.

Is this appropriate or is there another way of doing this?

$\endgroup$
  • $\begingroup$ Have you considered a conditional independence network? What would the plot you are suggesting have on the $x$- and $y$-axis? Adding the mean would just offset the location of each residual by the same amount from zero. $\endgroup$ – Frans Rodenburg Jun 28 '19 at 4:08
  • $\begingroup$ Adding the mean would just give some interpretability. On the x I would have blood pressure and on the y i would have adjusted body fat ( for age and sex) $\endgroup$ – mindhabits Jun 28 '19 at 4:16
  • $\begingroup$ If you want to do that, you should add the intercept instead, see my answer. $\endgroup$ – Frans Rodenburg Jun 28 '19 at 5:50
3
$\begingroup$

Since you only have two variables you're adjusting for, one of which is a 2-category factor, the easiest way to visualize the whole thing would just be by using some colouring and symbols:

# simulate data
set.seed(1234)
age <- round(rnorm(30, 40, 10))
sex <- round(runif(30, 0, 1))
bodyfat <- rnorm(30, 10 + 5 * sex + age / 10, 3)
bloodpressure <- rnorm(30, 55 + bodyfat * 2 + age / 10, 5)

# plot original with color and pch set to other variables
cols <- colorRampPalette(c("blue", "black", "red"))(30)
par(mfrow = c(1, 1))
plot(bloodpressure ~ bodyfat, col = cols[order(age)], pch = sex + 17)
legend("bottomright", legend = c("male", "female"), pch = 17:18)

original

You could add a colour scale to your legend for age.


Alternatively, you could adjust bodyfat for age and sex, although this loses its easy interpretation. If you want to adjust as you propose in your question—by regressing body fat on age and sex—you should add the intercept of this model to the residuals, not the mean.

This is equivalent to subtracting $\beta_1 \cdot \text{age}$ and $\beta_2 \cdot \text{sex}$ from the body fat percentage in:

$$\text{body fat percentage} = \beta_0 + \beta_1 \cdot \text{age} + \beta_2 \cdot{sex} + \epsilon$$

That would look like this:

# plot blood pressure vs adjusted bodyfat
LM <- lm(bodyfat ~ age + sex)
adjusted_BF <- (coef(LM)[1] + resid(LM))
par(mfrow = c(1, 2))
plot(bloodpressure ~ bodyfat)
plot(bloodpressure ~ adjusted_BF)

adjusted

As expected, the relationship observed in the right scatter plot is slightly weaker after adjusting for the effects of age and sex.


However, looking at partial correlations is more like adjusting both variables for the remaining variables:

LM2 <- lm(bloodpressure ~ age + sex)
adjusted_BP <- coef(LM2)[1] + resid(LM2)
plot(adjusted_BP ~ adjusted_BF)

partial

But again, be careful not to mislead your readers with this, as these axes no longer represent blood pressure and body fat.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your excellent explanation and detailing of the different methods. I have one follow up question -- I've came across some papers in my field that say this "in linear models adjusted for age and sex, blood pressure was associated with %bodyfat (r=0.45, P<0.01)". How did they obtain a Pearson's R from a linear model? Thanks for your insight. $\endgroup$ – mindhabits Jun 28 '19 at 17:54
  • 1
    $\begingroup$ By linear models adjusted for age and sex, they likely mean $\text{BP} = \beta_0 + \beta_1 \cdot \text{BF} + \beta_2 \cdot \text{age} + \beta_3 \cdot \text{sex} + \epsilon$. Authors often include the coefficient of determination ($R^2$) as a means for the reader to judge the variance in $\text{BP}$ explained by the model. Maybe this is what they meant? You can calculate it as 1 minus the residual sum of squares divided by the total sum of squares: en.wikipedia.org/wiki/Coefficient_of_determination $\endgroup$ – Frans Rodenburg Jun 28 '19 at 22:32
  • $\begingroup$ Understood, and that makes sense. What do you think of this response? $\endgroup$ – mindhabits Jun 29 '19 at 2:26
  • $\begingroup$ I think the responder's second suggestion makes sense. That is, "fix" the effects of the covariates/confounders using the mean and then calculate the Y values with values of my actual predictor of interest $\endgroup$ – mindhabits Jun 29 '19 at 2:32
  • $\begingroup$ Setting the covariates to a suitable value only makes sense if such a value exists. For age you could argue that most people are around a certain age (like the mean), but for sex, would it really make sense to take the mean of the dummy variable? $\endgroup$ – Frans Rodenburg Jun 29 '19 at 8:53
0
$\begingroup$

If you have enough data for each age-group you are interested in, "facets" may be a better way to convey your message than the first graph proposed by @Frans Rodenburg.

Simulate data

set.seed(1234)
age <- round(rnorm(1000, 40, 10))
sex <- round(runif(1000, 0, 1))
bodyfat <- rnorm(1000, 10 + 5 * sex + age / 10, 3)
bloodpressure <- rnorm(1000, 55 + bodyfat * 2 + age / 10, 5)

Data transformation

library(tidyverse)

data <- tibble(bodyfat, bloodpressure, age, sex)
data <- data %>%
          mutate(age_group = case_when(age < 25 ~ "0 to 25",
                                       age >= 25 & age < 40 ~ "25 to 40",
                                       age >= 40 & age < 50 ~ "40 to 50",
                                       age >= 50 & age < 60 ~ "50 to 60",
                                       age >= 60 ~ "60 and over"
                                       )) %>%
          mutate(sex = case_when(sex == 0 ~ "female",
                                 sex == 1 ~ "male"))

Data visualization through facets

ggplot(data, aes(bodyfat, bloodpressure)) +
  geom_point() +
  facet_wrap(age_group ~ sex)

enter image description here

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.