0
$\begingroup$

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!

$\endgroup$
  • 1
    $\begingroup$ Write sane loops? $\endgroup$ – Gregor --reinstate Monica-- Jul 18 '17 at 0:16
  • 2
    $\begingroup$ Why are you avoiding a for loop? This seems to call for one. $\endgroup$ – Matthew Drury Jul 18 '17 at 0:31
  • 5
    $\begingroup$ Don't worry about efficiency - just write something hacky that works. If it takes hours, then worry about efficiency. $\endgroup$ – Gregor --reinstate Monica-- Jul 18 '17 at 0:33
  • 4
    $\begingroup$ for loops get a bad rap; many times they are perfectly appropriate. $\endgroup$ – Mark White Jul 18 '17 at 2:13
  • 2
    $\begingroup$ 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. $\endgroup$ – Björn Jul 18 '17 at 6:00
4
$\begingroup$

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)]
#   r.squared adj.r.squared   p.value       AIC
# 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.

$\endgroup$
1
$\begingroup$

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
#         formula_string r.squared    adj.r.squared     sigma
#         <chr>          <dbl>        <dbl>             <dbl>
#   1     Y ~ X * Z1     0.2150010   -0.17749851        0.2024755
#   2     Y ~ X * Z2     0.3910735    0.08661019        0.1783283
$\endgroup$
  • $\begingroup$ extremely sweet $\endgroup$ – ℕʘʘḆḽḘ Jul 18 '17 at 2:33
  • $\begingroup$ is it possible to combine this purr solution with lm + coeftest? I am interested in the output of coeftest actually $\endgroup$ – ℕʘʘḆḽḘ Jul 18 '17 at 2:35
  • 1
    $\begingroup$ Yea, absolutely. What's your sample code? $\endgroup$ – Rahul Jul 18 '17 at 2:37
  • $\begingroup$ also, in your original solution it is hard to actually know which model is which. how can we keep track of this? can you store the model formula in the dataframe? $\endgroup$ – ℕʘʘḆḽḘ Jul 18 '17 at 2:38
  • $\begingroup$ Done. Here's sample code for you. $\endgroup$ – Rahul Jul 18 '17 at 2:45

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