Lecture 19: A Few R Demos, Interactions and Collinearity

BIOS 600 - Spring 2026

Authors

Affiliation

Announcements

Exam Corrections due Thursday 4/16 at 11:59pm,but email me if you need extra time.
Lab 8 due this Friday at 11:59pm.

Overview

Some R demos, helpful project things
Interpreting categorical variables, plus releveling factors
Interactions & how to interpret interaction terms
Checking for multicollinearity

Computational setup

library(tidyverse)
library(tidymodels)
library(kableExtra)
library(car) # for VIF

R Demo: merging datasets

Let’s say you have two data sets..

# Create two sample data frames
df1 <- data.frame(
  ID = c(1, 2, 3, 4),
  Name = c("Alice", "Bob", "Charlie", "David")
)

df2 <- data.frame(
  ID = c(1, 3, 5),
  City = c("New York", "London", "Paris")
)

df1

  ID    Name
1  1   Alice
2  2     Bob
3  3 Charlie
4  4   David

df2

  ID     City
1  1 New York
2  3   London
3  5    Paris

R Demo: Merging datasets

As long as you have an identifying column (e.g. studyid, patient number, observation number etc.), you can merge using left_join():

# The 'by' argument specifies the 
# common column(s) to join on.
merged_df <- df1 |>
  left_join(df2, by = "ID")

merged_df

  ID    Name     City
1  1   Alice New York
2  2     Bob     <NA>
3  3 Charlie   London
4  4   David     <NA>

R Demo: Dropping NA values

We haven’t talked about missing data much, but it may be simplest to conduct a “complete case analysis” where you only perform an analysis on complete cases:

merged_df |>
  drop_na()

  ID    Name     City
1  1   Alice New York
2  3 Charlie   London

R Demo: `write.csv()`

You can save your cleaned/merged dataset at any point using write.csv()

write.csv(merged_df, file = "data/data_clean.csv", 
          row.names = F) 
# do not create a column of row names

Note that you may need to include row.names = F so that R does not create an additional column of row numbers.

Categorical variables from last time…

Consider the model

\[Obesity_i = \beta_0 + \beta_1(HDI == Middle)_i + \]

\[ \beta_2(HDI == BottomTen)_i + \epsilon_i\]

The parameter interpretations are below.

\(\beta_0\) represents the expected obesity percentage for a state with 0 for the two dummy variables. That is, in the top ten HDI
\(\beta_1\) represents the expected difference in obesity percentage for a state in the middle HDI category, compared to the top ten
\(\beta_2\) represents the expected difference in obesity percentage for a state in the bottom ten HDI category, compared to the top ten

Fitting the model

fit2 <- lm(Obesity ~ HDI, data = cdc)
tidy(fit2) |>
  kable(digits = 3)

term	estimate	std.error	statistic
(Intercept)	34.430	0.816	42.195
HDIMiddle	-4.797	0.942	-5.091
HDITop ten	-8.080	1.154	-7.002

Wait a minute! R automatically assigned the lowest level of HDI as the reference level. This is because R assigns factor levels to characters alphabetically, so HDIBottom Ten was before HDIMiddle and HDITop Ten.
Let’s see how we can change the reference level. But before doing that…

Writing out the model

Write out the new model that has been fit. (i.e. \(Obesity_i = \ldots\)). Do we include an error term? Why or why not?

Changing the reference level

cdc_new <- cdc |>
  mutate(HDI_relevel = fct_relevel(HDI, "Top ten"))

This reorders the factor levels so that “Top ten” becomes the first (reference) level.
Now, when we run the regression, the model will automatically use Top ten as the baseline.

Re-running the model

fit2_relevel <- lm(Obesity ~ HDI_relevel, 
                   data = cdc_new)
tidy(fit2_relevel) |>
  kable(digits = 3)

term	estimate	std.error	statistic	p.value
(Intercept)	26.350	0.816	32.293	0.000
HDI_relevelBottom ten	8.080	1.154	7.002	0.000
HDI_relevelMiddle	3.283	0.942	3.485	0.001

Interpretation and Model

What is the interpretation of \(\hat\beta_2 = 3.283\) from the output above?

Interactions

Sometimes, the relationship between one predictor and the outcome depends on the value of another predictor variable. For example, the effect of exercise % on obesity may be different for smokers vs. non-smokers. To model such a relationship, we create an interaction term.
This is created simply by multiplying two predictors \(x_1\) and \(x_2\) to create a new predictor, \(x_1x_2\). When interaction terms are in a model, interpretations can become tricky.

Interactions

Let’s consider a model with main effects of exercise and smoking and an interaction term between them:

\[Obesity_i = \beta_0 + \beta_1 Exercise + \beta_2 Smoking\]

\[+ \beta_3 (Exercise_i \cdot Smoking_i) + \epsilon_i\]

What is the expected change in obesity % given a one percentage point increase in exercise %?

Interactions model

fit_int <- lm(Obesity ~ Exercise + 
                Smoking + 
                Exercise*Smoking, 
                   data = cdc)
tidy(fit_int) |>
  kable(digits = 3)

term	estimate	std.error	statistic	p.value
(Intercept)	38.577	13.867	2.782	0.008
Exercise	-0.358	0.270	-1.325	0.192
Smoking	0.408	0.729	0.559	0.579
Exercise:Smoking	0.003	0.015	0.186	0.853

Interpretation: If exercise increases by 1 percentage point, the expected change in obesity % depends on smoking %, and is: -.358 + .003(Smoking%)

More interpretation

So,

When Smoking = 0%, the slope of Exercise is -0.358 (exercise is assoc. w/ lower obesity).
When Smoking = 50%, slope = -0.358 + 0.003(50) = -0.208 (still negative association, but smaller in magnitude).

However: The interaction term is not significant. So we should remove the interaction term and just use the model with main effects only.

The math behind it

How did we get that interpretation? Here’s the math behind it…

For a given exercise (ex) and smoking (sm) percentage, our predicted obesity (ob) percentage is

\[\hat {ob}_i = \hat \beta_0 + \hat \beta_1 ex_i + \hat \beta_2 sm_i + \hat \beta_3 ex_i sm_i\]

and for a state at the same smoking % but 1 percentage point higher in exercise %, the predicted obesity percentage is

\[\hat {ob}_{i'} = \hat \beta_0 + \hat \beta_1 (ex_i + 1) + \hat \beta_2 sm_i + \hat \beta_3 (ex_i + 1) sm_i\]

\[ = \hat \beta_0 + \hat \beta_1 ex_i + \hat \beta_1 + \hat \beta_2 sm_i + \hat \beta_3 ex_i sm_i + \hat \beta_3 sm_i\]

The math behind it

Subtracting, we have \(\hat ob_{i'} - \hat ob_i= \hat \beta_1 + \hat \beta_3 sm_i\), which is the expected change in obesity % for a 1% change in exercise %.
Takeaway: In the interactions model, the relationship of exercise and obesity depends on the level of smoking in that state.

Quick review

Take a few minutes to review the math in small groups.

Interactions

Luckily, interpretation of interaction terms with categorical predictors is easier than with continuous predictors.
Since categorical predictors are based on dummy variables, they can only take on the values of 0 or 1 in the model.
Again, an interaction effect implies that the regression coefficient for an explanatory variable would change depending on the value of another predictor (for instance, the relationship between exercise and obesity might depend on whether a state is in High, Medium, or Low HDI groups).

Multicollinearity

One common problem in multiple regression is multicollinearity, which occurs when multiple highly correlated variables are used as predictors.
Collinearity simply refers to two variables behind highly correlated.
In this case, the model can become unstable (often seen as standard errors that get huge and lead to huge confidence interval estimates), and it can be difficult to assess the relationships of the predictors.

Diagnosing multicollinearity

If nothing is significant when you “expect” something to be, we have some clues:

Individual predictors are significant in simple linear regression models,
but standard errors and interval estimates are huge,
and the overall F test is significant

A significant overall F test with no significant individual variable test is a typical sign of collinearity. We can check out the correlations among the three predictors.

Checking VIF

The common way to check for multicollinearity is through VIF values. VIF values indicate the extent of multicollinearity.
A VIF value of 1 suggests no correlation with other predictors.
Values between 1 and 5 generally indicate moderate correlation.
Values above 5 (or sometimes 10, depending on the context) are often considered problematic, indicating severe multicollinearity and potential issues with coefficient stability and interpretation.

VIF

Let’s say we had the following model:

fit4 <- lm(Obesity ~ Exercise + Smoking + Over65, 
           data = cdc)
vif_values <- vif(fit4) # car package
vif_values

Exercise  Smoking   Over65 
1.420894 1.541720 1.123272

In this model, none of the VIFs are above 5 (in fact, they’re much closer to 1), so multicollinearity is not an issue here.
Advice on what to do differs! Sometimes, removing a variable is more problematic than multicollinearity, sometimes, multicollinearity is worse.

Recap

Some R demos: merging, dropping NA, write.csv(), re-leveling a factor
Interactions, interpretations, and the math behind it all
Collinearity

Next class

Logistic regression!

Announcements

Overview

Computational setup

R Demo: merging datasets

R Demo: Merging datasets

R Demo: Dropping NA values

R Demo: write.csv()

Categorical variables from last time…

Fitting the model

Writing out the model

Changing the reference level

Re-running the model

Interactions

Interactions

Interactions model

More interpretation

The math behind it

The math behind it

Interactions

Multicollinearity

Diagnosing multicollinearity

Checking VIF

VIF

Recap

Next class

R Demo: `write.csv()`