 <p>You can do this easily with <code>dplyr</code>. Assuming your data frame is named <code>df</code>, run:</p> <pre><code>library(dplyr) library(forecast) model_fits &lt;- group_by(df, Region) %&gt;% do(fit=auto.arima(.$Sales)) </code></pre> <p>The result is a data frame containing the model fits for each region:</p> <pre><code>&gt; head(model_fits) Source: local data frame [6 x 2] Groups: &lt;by row&gt; Region fit 1 A &lt;S3:Arima&gt; 2 B &lt;S3:Arima&gt; 3 C &lt;S3:Arima&gt; 4 D &lt;S3:Arima&gt; 5 E &lt;S3:Arima&gt; 6 F &lt;S3:Arima&gt; </code></pre> <p>You can get a list with each model fit like so:</p> <pre><code>&gt; model_fits$fit [[1]] Series: .$Sales ARIMA(0,0,0) with non-zero mean Coefficients: intercept 196.0000 s.e. 14.4486 sigma^2 estimated as 2088: log likelihood=-52.41 AIC=108.82 AICc=110.53 BIC=109.42 [[2]] Series: .$Sales ARIMA(0,0,0) with non-zero mean Coefficients: intercept 179.2000 s.e. 14.3561 sigma^2 estimated as 2061: log likelihood=-52.34 AIC=108.69 AICc=110.4 BIC=109.29 </code></pre>