# R: how to chose best interaction variable in an OLS regression? [closed]

I have a simple regression setup:

$Y$ = $\alpha$ + $\beta X$ + $\gamma Z$ + $\delta XZ$ + $\epsilon$

As you can see, the variable X is interacted with Z and the interaction coefficient of interest is $\delta$

In my dataframe, I have many possible variables for Z and I would like to choose automatically the best one.

That is, in each regression $X$ is the same but pick the $Z$ whose regression R-square is highest or the one where the t-stat on $\delta$ is highest among all possible regressions.

Here is a simple example

library(dplyr)
set.seed(10)
dataframe = data_frame(Y = runif(10), X= runif(10),  Z1=runif(10), Z2 =runif(10))

> dataframe
# A tibble: 10 x 4
Y          X        Z1         Z2
<dbl>      <dbl>     <dbl>      <dbl>
1 0.50747820 0.65165567 0.8647212 0.53559704
2 0.30676851 0.56773775 0.6153524 0.09308813
3 0.42690767 0.11350898 0.7751099 0.16980304
4 0.69310208 0.59592531 0.3555687 0.89983245
5 0.08513597 0.35804998 0.4058500 0.42263761
6 0.22543662 0.42880942 0.7066469 0.74774647
7 0.27453052 0.05190332 0.8382877 0.82265258
8 0.27230507 0.26417767 0.2395891 0.95465365
9 0.61582931 0.39879073 0.7707715 0.68544451
10 0.42967153 0.83613414 0.3558977 0.50050323


How can I do that in R without writing crazy loops? Thanks!

• Write sane loops? Jul 18, 2017 at 0:16
• Why are you avoiding a for loop? This seems to call for one. Jul 18, 2017 at 0:31
• Don't worry about efficiency - just write something hacky that works. If it takes hours, then worry about efficiency. Jul 18, 2017 at 0:33
• for loops get a bad rap; many times they are perfectly appropriate. Jul 18, 2017 at 2:13
• You probably do not want to pick a single model on this way - it is essentially certain to cause (severe) overfitting, any p-values become onvalid etc. With the same set-up you can probably do model averaging. Jul 18, 2017 at 6:00

Something like this:

library(dplyr)
library(broom)
z_options = c("Z1", "Z2")
formula_string = paste("Y ~ X *", z_options)
model_list = list()
for(i in seq_along(formula_string)) {
model_list[[i]] = lm(as.formula(formula_string[i]), dataframe)
}

model_glance = bind_rows(lapply(model_list, glance))
model_glance[, c(1, 2, 5, 8)]
# 1 0.2150010   -0.17749851 0.6677345  1.327784
# 2 0.3910735    0.08661019 0.3620319 -1.212065


Use whatever metric you'd like to pick the best.

I tend to not like loops either. I'm a big fan of Hadley's purrr package. Clean, elegant, readable code.

Here's how to use to use it, built off of @Gregor's solution.

library(dplyr)
library(purrr)
library(broom)

set.seed(10)
dataframe = data_frame(
Y = runif(10),
X = runif(10),
Z1 = runif(10),
Z2 = runif(10)
)

z_options = c("Z1", "Z2")
formula_string = paste("Y ~ X *", z_options)

tibble(formula_string) %>%
bind_cols(
map_df(.x = formula_string,
.f = ~ glance(lm(as.formula(.x),
dataframe)))
) %>%
select(1:4)

# A tibble: 2 x 4

• is it possible to combine this purr solution with lm + coeftest? I am interested in the output of coeftest actually Jul 18, 2017 at 2:35