Run Multiple Regression Models in for-Loop in R (Example)

 

In this article, I’ll show how to estimate multiple regression models in a for-loop in the R programming language.

Table of contents:

If you want to know more about these topics, keep reading…

 

Introducing Example Data

The following data is used as basement for this R programming tutorial:

set.seed(98274)                          # Creating example data
y <- rnorm(1000)
x1 <- rnorm(1000) + 0.2 * y
x2 <- rnorm(1000) + 0.2 * x1 + 0.1 * y
x3 <- rnorm(1000) - 0.1 * x1 + 0.3 * x2 - 0.3 * y
data <- data.frame(y, x1, x2, x3)
head(data)                               # Head of data
#            y         x1         x2          x3
# 1  0.5587036 -0.3779533 -0.5320515 -0.92069263
# 2  0.8422515 -1.3835572  1.2782521  0.87967960
# 3 -0.5395343 -0.9729798 -0.1515273 -0.05973894
# 4 -0.3522260  1.2977564 -0.3512013 -0.77239810
# 5  1.5848675 -1.3152806 -2.3644414 -1.14651812
# 6  0.2207957  1.8860636  0.1967851 -0.04963894

As you can see based on the previous RStudio console output, our example data consists of four numeric columns. The first variable is our regression outcome and the three other variables are our predictors.

 

Example: Running Multiple Linear Regression Models in for-Loop

In this Example, I’ll show how to run three regression models within a for-loop in R. In each for-loop iteration, we are increasing the complexity of our model by adding another predictor variable to the model.

First, we have to create a list in which we will store the outputs of our for-loop iterations:

mod_summaries <- list()                  # Create empty list

Now, we can write a for-loop that runs multiple linear regression models as shown below:

for(i in 2:ncol(data)) {                 # Head of for-loop
 
  predictors_i <- colnames(data)[2:i]    # Create vector of predictor names
  mod_summaries[[i - 1]] <- summary(     # Store regression model summary in list
    lm(y ~ ., data[ , c("y", predictors_i)]))
 
}

Let’s have a look at the output of our previously executed for-loop:

mod_summaries                            # Return summaries of all models

 

multiple linear regressions executed in a for-loop in r

 

As you can see in Figure 1, we have created a list containing three different summary statistics of three different linear regressions.

 

Video, Further Resources & Summary

If you need further explanations on the content of this tutorial, I can recommend having a look at the following video that I have published on my YouTube channel.

In the video, I’m explaining how top loop and repeat the estimation of multiple regression models using the R programming syntax of this tutorial in RStudio.

 

The YouTube video will be added soon.

 

In addition, you may want to read the other R programming tutorials of my homepage.

 

Summary: At this point you should know how to write a for-loop executing several linear regressions in R programming. Please let me know in the comments below, in case you have further questions. Furthermore, please subscribe to my email newsletter to receive updates on new articles.

 

Subscribe to the Statistics Globe Newsletter

Get regular updates on the latest tutorials, offers & news at Statistics Globe.
I hate spam & you may opt out anytime: Privacy Policy.


28 Comments. Leave new

  • Fernando Ramires
    March 1, 2021 1:46 pm

    the procedure can follow the stepwise line….

    Reply
  • Farieda Hassanen
    July 24, 2021 5:22 pm

    Very good

    Reply
  • Thanks for the excellent explanation. If you can provide steps on how to use PCA (principal component analysis) on series of images, that would be great!

    Reply
  • How can I modify this to have models for each variable by itself? For example I want the the second model to just have x2 and the third model to just have x3 so on, rather than them building onto each other. Also could I change this to do and OddsRatio with each predictor?

    Reply
    • Hey Helen,

      To change the predictor variables as you want, you only have to change

      predictors_i <- colnames(data)[2:i]

      to

      predictors_i <- colnames(data)[i]

      (i.e. remove the “2:”).

      I don’t have a tutorial on Odds Ratios yet. However, I think this article explains it nicely: https://jarrettmeyer.com/2019/07/23/odds-ratio-in-r

      Regards

      Joachim

      Reply
  • Hi Mr Joachim
    I need your help because till now no one answered my question and I hope to find the solution with you,

    How can find the range (minimum, maximum) values for Intercept, slope and R-square in r programming?

    Best Regards
    Hussein Ali

    Reply
  • Simon von Zahn
    October 28, 2021 7:31 pm

    Hello Joachim,
    will this example provide me knowing which of the tested models to use for a maximised adj. R squared?
    And which packages are required to run it?
    Best Regards
    Simon

    Reply
  • Hi Joachim,

    I’m trying to use this method to run regressions on “n” independent variables, so up to every combination available from my columns. Is there a way I can do this with this code? Thanks for your help!

    Reply
  • Hi Joachim,
    thank you for sharing this guide. I run several linear regression models with it. How can I extract all p-values and estimates of the list of models I created? I would like to store them in one dataframe.

    Reply
  • How can we extract the predicted value from this linear models that stored in the list

    Reply
  • Hii Joachim,

    Thanks for the nice example. How can I expand this code to extract the coeff, SD and P-value and put the in a separate list?

    Reply
  • Hey Joachim,
    Can I modify the code and use for multiple ‘target variables’ instead of multiple ‘predictors’? If so, could you examplify the code.
    Thanks

    Reply
    • Hey Sal,

      Yes, this is possible. Please have a look at the modified example code below:

      set.seed(98274)                          # Creating example data
      y <- rnorm(1000)
      x1 <- rnorm(1000) + 0.2 * y
      x2 <- rnorm(1000) + 0.2 * x1 + 0.1 * y
      x3 <- rnorm(1000) - 0.1 * x1 + 0.3 * x2 - 0.3 * y
      data <- data.frame(y1 = y,
                         y2 = y + rnorm(1000),
                         y3 = y + 10 * rnorm(1000),
                         x1, x2, x3)
      head(data)                               # Head of data
      #           y1          y2         y3         x1         x2          x3
      # 1  0.5587036  1.53788165  0.1669674 -0.3779533 -0.5320515 -0.92069263
      # 2  0.8422515  0.08528937  0.6847646 -1.3835572  1.2782521  0.87967960
      # 3 -0.5395343 -0.98086394 -1.2222470 -0.9729798 -0.1515273 -0.05973894
      # 4 -0.3522260  0.30012623 -1.7663313  1.2977564 -0.3512013 -0.77239810
      # 5  1.5848675 -0.48901240 -1.3411279 -1.3152806 -2.3644414 -1.14651812
      # 6  0.2207957  0.11100091 -0.2485206  1.8860636  0.1967851 -0.04963894
       
      mod_summaries <- list()                  # Create empty list
       
      out_names <- c("y1", "y2", "y3")         # Names of all target variables
       
      for(i in 1:length(out_names)) {          # Estimate models
       
        mod_summaries[[i]] <- summary(lm(get(out_names[i]) ~ x1 + x2 + x3, data))
      }
       
      mod_summaries                            # Return summaries of all models
      # [[1]]
      # 
      # Call:
      # lm(formula = get(out_names[i]) ~ x1 + x2 + x3, data = data)
      # 
      # Residuals:
      #      Min       1Q   Median       3Q      Max 
      # -2.67110 -0.65829 -0.03136  0.66088  3.15418 
      # 
      # Coefficients:
      #              Estimate Std. Error t value Pr(>|t|)    
      # (Intercept)  0.003017   0.030502   0.099    0.921    
      # x1           0.131145   0.031817   4.122 4.07e-05 ***
      # x2           0.173592   0.031297   5.547 3.73e-08 ***
      # x3          -0.211426   0.029224  -7.235 9.31e-13 ***
      # ---
      # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      # 
      # Residual standard error: 0.9631 on 996 degrees of freedom
      # Multiple R-squared:  0.09867,	Adjusted R-squared:  0.09596 
      # F-statistic: 36.35 on 3 and 996 DF,  p-value: < 2.2e-16
      # 
      # 
      # [[2]]
      # 
      # Call:
      # lm(formula = get(out_names[i]) ~ x1 + x2 + x3, data = data)
      # 
      # Residuals:
      #     Min      1Q  Median      3Q     Max 
      # -4.3833 -0.9343  0.0150  0.8776  5.3319 
      # 
      # Coefficients:
      #             Estimate Std. Error t value Pr(>|t|)    
      # (Intercept) -0.03054    0.04400  -0.694  0.48784    
      # x1           0.11453    0.04590   2.495  0.01274 *  
      # x2           0.13429    0.04515   2.974  0.00301 ** 
      # x3          -0.21587    0.04216  -5.121 3.66e-07 ***
      # ---
      # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      # 
      # Residual standard error: 1.389 on 996 degrees of freedom
      # Multiple R-squared:  0.04237,	Adjusted R-squared:  0.03948 
      # F-statistic: 14.69 on 3 and 996 DF,  p-value: 2.301e-09
      # 
      # 
      # [[3]]
      # 
      # Call:
      # lm(formula = get(out_names[i]) ~ x1 + x2 + x3, data = data)
      # 
      # Residuals:
      #     Min      1Q  Median      3Q     Max 
      # -34.264  -6.389  -0.109   6.519  34.821 
      # 
      # Coefficients:
      #             Estimate Std. Error t value Pr(>|t|)   
      # (Intercept) -0.44531    0.31494  -1.414  0.15769   
      # x1           0.33912    0.32852   1.032  0.30220   
      # x2           0.85905    0.32315   2.658  0.00798 **
      # x3          -0.09715    0.30175  -0.322  0.74755   
      # ---
      # Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
      # 
      # Residual standard error: 9.944 on 996 degrees of freedom
      # Multiple R-squared:  0.0105,	Adjusted R-squared:  0.007517 
      # F-statistic: 3.522 on 3 and 996 DF,  p-value: 0.01465

      I hope this helps!

      Joachim

      Reply
  • The last example is really useful for running multiple regressions in a for loop. What if I wanted to extract the regression coefficients for each i in out_names (Y1, Y2, Y3), and then transform those? For example, for each target variable (Y1, Y2, Y3) if I were interested in understanding the fraction of X1*x/Yi (such that X1*x/Yi +X2*x/Yi+X3*x/Yi = 1). Could I use a loop to extract the regression coefficient and store it so I can use it in a new function?

    Reply

Leave a Reply

Your email address will not be published. Required fields are marked *

Fill out this field
Fill out this field
Please enter a valid email address.

Menu
Top