title: Estimating a single model changing multiple variables
body: <p>I have a panel dataset and I am trying to estimate multiple models using different key independent variables. 
  My dataset looks something like this</p>
  
  <pre><code>df &lt;- data.frame(y = rnorm(9),id= c(1,1,1,2,2,2,3,3,3), time=c(1,2,3,1,2,3,1,2,3), x1 = rnorm(9), x2 = rnorm(9), x3 = rnorm(9), c1 = rnorm(9),  c2 = rnorm(9))
  
  df &lt;- pdata.frame(df, index=c("id","time"))
  </code></pre>
  
  <p>Where <em>y</em> is the key dependent variable, <em>id</em> is the unit identifier and <em>time</em> is the time identifier. <em>c1 and c2</em> are controls, and <em>x1,x2,x3</em> are key independent variables of interest. 
  I know that I can estimate these model separately like this:</p>
  
  <pre><code>library(plm)
  m1&lt;- plm(y ~  lag(x1, k = 1)+ lag(y,k=1) + lag(c1, k = 1) + 
    lag(c2, k = 1), data = df, effect = "twoways", 
  model = "within")
  
  m2&lt;- plm(y ~  lag(x2, k = 1)+ lag(y,k=1) + lag(c1, k = 1) + 
    lag(c2, k = 1), data = df, effect = "twoways", 
  model = "within")
  </code></pre>
  
  <p>And so on...</p>
  
  <p>However, I was looking for a more concise way to write and examine different models since I have 16 independent variables of interest and multiple dependent variables to examine. 
  I would like to:</p>
  
  <p>1) Find a way to "group" all of the independent variables. So that I can write the model only once (the controls and dv remain the same) and r computes all of the models for each independent variable.</p>
  
  <p>2) Ideally, the command should also immediately display all of the results in an orderly fashion (stargazer?), so that I can quickly examine which model performs better and which one performs worse. </p>
  
  <p>Thank you a lot for your help</p>