I have multiple dataframes and I would like to evaluate (multiple) different models on each. MWE

```
df1 <- data.frame(A3 = c(-5, 5, 1), 
                  B3 = c(0, 10, 1))
df2 <- data.frame(A4 = c(5, 15, 1))
                  B4 = c(10, 20, 1))

myfun <- function(arg1, arg2){ # arg1 =1 or 2
  if (arg2 == 1){
    eqn <- paste0("A", arg1+2) ~ paste0("B", arg1+2) + I(as.name(paste0("B", arg1+2))^2)
  } else {
    eqn <- paste0("A", arg1+2) ~ paste0("B", arg1+2) + I(as.name(paste0("B", arg1+2))^2) +I(as.name(paste0("B", arg1+2))^3)
  }

  return (lm(formula = eqn, data = eval(as.name(paste0("df", arg1)))
    )
  )
}

```

For example if I run `myfun(1,2)` I would like to get `lm(A4 ~ B4 + I(B4^2) + I(B4^3), data = df2)`. But whatever I try I get the following error message `Error in (paste0("B", arg1 +2))^2 : non-numeric argument to binary operator`. From what I read in `?I`, I imagine this is because R isolates whatever is passed into `I()`, so it doesn't realize I am trying to transform a variable: is that what is going on, and is it something I can fix? Also, is there a better way to estimate multiple models quickly? All the similar questions I found used the same `data.frame` across models, while I have to account for the response (and predictor) variables coming from different dataframes for different models.