Run nested function for multiple regressions on independent variables for each outcome variable and plot coefficients graphically

I am trying to run a set of nested functions in order to perform multiple regressions and graph the results. For the regression, I have a list of outcome variables, and a list of independent variables. I want to run each independent variable as an interaction term with in separate regressions, and do this for each outcome variable. My intention is to perform the above, store the results as dataframes, then create ggplot graphs for each outcome variable, with each independent variable as a different predicted line.

My code as of right now would be something like the following (modelled using the mtcars dataset)

library(ggplot2)
library(tidyverse)

outcome_var_list = mtcars %>% select(qsec,hp) %>% names()
var_list = mtcars %>% select(cyl,wt) %>% names()

iterate <- sapply(outcome_var_list,function(x){(one_df = sapply(var_list,function(k){

mod <- lm(data=mtcars, paste(x," ~ ",k,":factor(drat)",sep = ""))

mod_df <- as.data.frame(coef(mod))

})
  
#using rownames and coefficient as placeholder names for the actual column names
mod_df %>% ggplot(aes(x = factor(mod_df), y= coefficient, color = variable)) + geom_line()

})

The above does not successfully run. Problems I encounter are that the coefficients do not seem to store as dataframes within the sapply function, and only take on a matrix. Is it perhaps not the right case in which to use nested functions? Or, is there something else I must change to execute this properly? I would like to store each line plot the function produces, and title it with the corresponding outcome variable. Is this possible, or is there another approach to writing a function that I should consider?

EDIT

Here is the code to generate a single plot for only 1 outcome variable and only 2 independent variables (manually)

mod <- lm(data= mtcars, qsec ~ cyl:factor(drat))
mod_df <- as.data.frame(coef(mod))
mod_df <- tibble::rownames_to_column(mod_df)
#rename rowname
mod_df <- mod_df %>% mutate(rowname = str_replace(mod_df$rowname,pattern = "cyl:factor[(]drat[)]",replacement = ""))
#rename 
names(mod_df) <- c("xvar","yvar")
#remove intercept
mod_df <- mod_df[2:nrow(mod_df),]
#add label for independent variable
mod_df$color = "cyl"
  
#repeat for second independent variable

mod2 <- lm(data= mtcars, qsec ~ wt:factor(drat))
mod2_df <- as.data.frame(coef(mod2))
mod2_df <- tibble::rownames_to_column(mod2_df)
#rename rowname
mod2_df <- mod2_df %>% mutate(rowname = str_replace(mod2_df$rowname,pattern = "wt:factor[(]drat[)]",replacement = ""))
#rename 
names(mod2_df) <- c("xvar","yvar")
#remove intercept
mod2_df <- mod2_df[2:nrow(mod2_df),]
#add label for independent variable
mod2_df$color = "wt"

combined_df <- rbind(mod_df,mod2_df)
combined_df %>% ggplot(aes(x=as.numeric(xvar), y=yvar, color = color)) + geom_line() + ggtitle("qsec") + theme_light() + theme(plot.title = element_text(hjust = 0.5))

The resulting plot is the following:

Is there a way to automate this over multiple outcome variables and independent variables for each one?

The code to obtain all the necessary coefficients will be:

mtcars %>%
  select(all_of(c(outcome_var_list, var_list)), drat)%>%
  pivot_longer(all_of(var_list))%>%
  reframe(broom::tidy(lm(as.matrix(pick(all_of(outcome_var_list)))~
               name/factor(drat):value + 0, data = .)))%>%
  select(response, term, estimate)%>%
  mutate(color = str_extract(term, "(?<=name)[^:]+"),
         xvar = as.numeric(str_extract(term, "[0-9.]+")))%>%
  drop_na()

Check all the values results generated against your values. They should match.

eg if you only need the coefficients regarding qsec you could filter the dataframe.

base R code:

fm <- sprintf("cbind(%s)~time/factor(drat):%s + 0", 
              toString(outcome_var_list), var_list[1])

cfs <- subset(mtcars, select = c(outcome_var_list, var_list, "drat"))|>
  reshape(varying = list(var_list), dir = 'long', times = var_list) |>
  lm(fm, data=_) |>
  coef()

final_df <- na.omit(data.frame(xvar = as.numeric(trimws(rownames(cfs),,"\\D")),
      yvar = c(cfs), color = sub("time(.*?):.*", '\', rownames(cfs)),
      response = colnames(cfs)[col(cfs)]))

  xvar       yvar color response
3  2.76  -1.286761   cyl     qsec
4  2.76  -2.701420    wt     qsec
5  2.93  -1.189608   cyl     qsec
6  2.93  -1.900810    wt     qsec
7  3.00  -1.209608   cyl     qsec
8  3.00  -1.869331    wt     qsec
9  3.07  -1.228775   cyl     qsec
10 3.07  -2.664112    wt     qsec
11 3.08  -1.319161   cyl     qsec
12 3.08  -2.760143    wt     qsec
13 3.15  -1.292108   cyl     qsec
14 3.15  -3.141629    wt     qsec
15 3.21  -1.457108   cyl     qsec
16 3.21  -3.394749    wt     qsec
17 3.23  -1.259608   cyl     qsec
18 3.23  -1.971797    wt     qsec
19 3.54  -1.612108   cyl     qsec

Compare the results to your combinded df.

Note that the method provided is generic