# Generate fake data
a <- as.factor(sample(0:1, 20, replace = TRUE))
b <- as.factor(sample(c("a","b","c","d","e","f"), 20, replace = TRUE))
c <- as.factor(sample(0:10, 20, replace = TRUE))
d <- as.factor(sample(0:12, 20, replace = TRUE))
y <- rnorm(20)
df <- data.frame(y,a,b,c,d)

# The factor variable names are:
vars <- c("a","b","c", "d")

# Loop through all the factors
for (i in 1:(length(vars) - 1)){
  for (j in (i+1):length(vars)){

    # Generate the right-hand side of the formula using
    # the fact that (x+y+z)^2 expands in the lm() formula
    # to all main and interaction terms for all two-way
    # interactions: (x + y + z + x:y + x:z + z:y)
    rhs <- c(vars[i], vars[j]) %>%
      paste(., collapse = "+") %>%
      paste0("(", ., ")", "^2")

    # Generate left-hand side
    lhs <- paste0("y", " ~ ")

    # Generate the model formula
    my_mod <- paste0(lhs, rhs) %>%
      formula()

    # Fit the model, save coefficients
    mod_sum <- lm(my_mod, data = df)
    mod_coef <- mod_sum$coefficients

    # Identify interaction coefficients by the ":"
    # symbol and keep only the interaction terms in a df
    int_coefs_df <- mod_coef %>%
      names() %>%
      grep(":",.) %>%
      mod_coef[.] %>%
      data.frame(estimate = .)

    print(int_coefs_df)
  }
}