<p>You can do this by combining <a href="http://cran.r-project.org/web/packages/dplyr/index.html" rel="nofollow">dplyr</a>, <a href="http://cran.r-project.org/web/packages/tidyr/index.html" rel="nofollow">tidyr</a> and my <a href="https://github.com/dgrtwo/broom" rel="nofollow">broom</a> package (you can install them with <code>install.packages</code>). First you need to gather all the numeric columns into a single column:</p>

<pre><code>library(dplyr)
library(tidyr)
tidied &lt;- myDat %&gt;%
    gather(column, value, -X, -Recipe, -Step, -Stage, -Prod)
</code></pre>

<p>To understand what this does, you can read up on <a href="http://blog.rstudio.org/2014/07/22/introducing-tidyr/" rel="nofollow">tidyr's gather operation</a>. (This assumes that all columns besides X, Recipe, Step, Stage, and Prod are numeric and therefore should be predicted in your regression. If that's not the case, you need to remove them beforehand. You'll need to produce a reproducible example of the problem if you need a more customized solution).</p>

<p>Then perform each regression, while grouping by the column and the four grouping variables.</p>

<pre><code>library(broom)

regressions &lt;- tidied %&gt;%
    group_by(column, Recipe, Step, Stage, Prod) %&gt;%
    do(mod = lm(value ~ X))

glances &lt;- regressions %&gt;% glance(mod)
</code></pre>

<p>The resulting <code>glances</code> data frame will have one row for each combination of column, Recipe, Step, Stage, and Prod, along with an <code>r.squared</code> column containing the R-squared from each model. (It will also contain <code>adj.r.squared</code>, along with other columns such as F-test p-value: see <a href="https://github.com/dgrtwo/broom#tidying-functions" rel="nofollow">here</a> for more). Running <code>coefs &lt;- regressions %&gt;% tidy(mod)</code> will probably also be useful for you, as it will get the coefficient estimates and p-values from each regression.</p>

<p>A similar use case is described in the <a href="http://cran.r-project.org/web/packages/broom/vignettes/broom_and_dplyr.html" rel="nofollow">"broom and dplyr" vignette</a>, and in Section 3.1 of <a href="http://arxiv.org/abs/1412.3565" rel="nofollow">the broom manuscript</a>.</p>
Regression of a Data Frame with multiple factor groupings
<p>I am working on a regression script.
I have a data.frame with roughly 130 columns, of which I need to do a regression for one column (lets call it X column) against all the other ~100 numeric columns.  </p>

<p>Before the regression is calculated, I need to group the data by 4 factors: <code>myDat$Recipe</code>, <code>myDat$Step</code>, <code>myDat$Stage</code>, and <code>myDat$Prod</code> while still keeping the other ~100 columns and row data attached for the regression. Then I need to do a regression of each column ~ X column and print out the R^2 value with the column name. This is what I've tried so far but it is getting overly complicated and I know there's got to be a better way.</p>

<pre><code> rm(list=ls())
 myDat &lt;- read.csv(file="C:/Users/Documents/myDat.csv",              header=TRUE, sep=",")

for(j in myDat$Recipe)
{
  myDatj &lt;- subset(myDat, myDat$Recipe == j) 
  for(k in myDatj$Step)
  {
    myDatk &lt;- subset(myDatj, myDatj$Step == k) 
    for(i in myDatk$Stage)
    {
      myDati &lt;- subset(myDatk, myDatk$Stage == i)
      for(m in myDati$Prod)
      {
        myDatm &lt;- subset(myDati, myDati$Prod == m)
          if(is.numeric(myDatm[3,i]))  
          {     
          fit &lt;- lm(myDatk[,i] ~ X, data=myDatm) 
          rsq &lt;- summary(fit)$r.squared
            {
              writeLines(paste(rsq,i,"\n"))
           }  
         }
      }
    }
  }  
}      
</code></pre>