I am estimating many linear models with interactions between two factor variables drawn from a large set of factor variables as predictors. Factors may differ in their number of levels and manually counting the number of levels is time-intensive. From these models, I am trying to generate a data frame containing a row for each coefficient estimate of interaction terms only. The coefficient estimates in the lm object are stored in a named vector but I store them in a df.

Currently, I have designed the function to generate a data frame containing these terms for each model. I could do this, save the results, then eventually read in and append/join the data frames but I know this is very slow and inefficient.

Does anyone have an idea for how to calculate the number of interaction terms that will be in the model then store the results in a data frame I've generated to fit the number of results?

Below is a minimal working example where I print the data frames after generating them. I could have saved them as well.

```
# 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)
  }
}

```