Nested loop with simple linear regression

Question

I am aware of several posts on this topic, see here Print and export loop from simple linear regression and here How to Loop/Repeat a Linear Regression in R.

However, I am interested in looping through not only a set of predictors but also a set of outcomes. Please see my code attempt below.

set.seed(42)
n <- 100
age <- runif(n, min=45, max=85)
sex <-  factor(sample(c('Male','Female'), n, rep=TRUE, prob=c(.6, .4)))
smoke <- factor(sample(c("Never", 'Former', 'Current'), n, rep=TRUE, prob=c(.25, .6, .15)))
bmi <- runif(n, min=16, max=45)
outcome1 <- rbinom(n = 100, size = 1, prob = 0.3)
outcome2 <- rbinom(n = 100, size = 1, prob = 0.13)
predictorlist<- list("age","sex", "bmi", "smoke")
outcomelist <- list("outcome1","outcome2")

for (i in predictorlist){ 
  for (j in outcomelist){ 
    model <- lm(paste("j ~", i[[1]])) 
    print(summary(model)) 
  } 
}

In doing so I receive the following error message:

Error in model.frame.default(formula = paste("j ~", i[[1]]), drop.unused.levels = TRUE) : 
  variable lengths differ (found for 'age')

A quick check of their lengths reveals everything to be in order

> length(age)
[1] 100
> length(sex)
[1] 100
> length(bmi)
[1] 100
> length(smoke)
[1] 100
> length(outcome1)
[1] 100
> length(outcome2)
[1] 100

Any assistance is, naturally, tremendously appreciated.

@user2974951 thanks, any assistance? This seems to work for (i in predictorlist){ model <- lm(paste("outcome1 ~", i)) print(summary(model)) } but am interested in also looping through the outcomes. — Sandro
– Sandro, Commented Dec 7, 2021 at 10:08
Just put the j outside the quotation like : model <- lm(paste(j, "~", i[[1]])) ` — maydin
– maydin, Commented Dec 7, 2021 at 10:08
@maydin thank you worked a treat! Any advice on something I could add to my loop to enable quick export of the coefficients and SE for each model? — Sandro
– Sandro, Commented Dec 7, 2021 at 10:18

maydin · Accepted Answer · 2021-12-07 12:38:12Z

1

my_models <- list()
n <- 1
for (i in 1:length(predictorlist)){ 
  for (j in 1:length(outcomelist)){ 
    formula <- paste(outcomelist[[j]], "~", predictorlist[[i]])
    out <- coefficients(summary(model))
    Coeff <- rownames(out)  
    rownames(out) <- NULL 
    model <- lm(formula)
    my_models[[n]]<-data.frame(Model=formula ,Coeff = Coeff,out )

    n <- n+1
  } 
}

do.call(rbind,my_models)

gives,

              Model       Coeff      Estimate  Std..Error       t.value     Pr...t..
1    outcome1 ~ age (Intercept)  5.551115e-16 0.085820144  6.468313e-15 1.0000000000
2    outcome1 ~ age smokeFormer  1.756757e-01 0.094117066  1.866566e+00 0.0649830089
3    outcome1 ~ age  smokeNever -1.638898e-16 0.131940938 -1.242145e-15 1.0000000000
4    outcome2 ~ age (Intercept)  4.997536e-01 0.239278415  2.088586e+00 0.0393369696
5    outcome2 ~ age         age -3.936905e-03 0.003567878 -1.103430e+00 0.2725424813
6    outcome1 ~ sex (Intercept)  3.600995e-01 0.188106966  1.914334e+00 0.0584945208
7    outcome1 ~ sex         age -3.487458e-03 0.002804861 -1.243362e+00 0.2167006313
8    outcome2 ~ sex (Intercept)  2.173913e-01 0.063533308  3.421690e+00 0.0009091606
9    outcome2 ~ sex     sexMale  4.186795e-02 0.086457881  4.842584e-01 0.6292829854
10   outcome1 ~ bmi (Intercept)  1.304348e-01 0.050088617  2.604080e+00 0.0106444776
11   outcome1 ~ bmi     sexMale -8.051530e-04 0.068161974 -1.181235e-02 0.9905993425
12   outcome2 ~ bmi (Intercept)  3.919620e-01 0.159188779  2.462246e+00 0.0155528798
13   outcome2 ~ bmi         bmi -4.967139e-03 0.005010602 -9.913259e-01 0.3239678274
14 outcome1 ~ smoke (Intercept)  1.085051e-01 0.125958776  8.614332e-01 0.3911023372
15 outcome1 ~ smoke         bmi  7.025988e-04 0.003964659  1.772154e-01 0.8597049249
16 outcome2 ~ smoke (Intercept)  3.333333e-01 0.111424058  2.991574e+00 0.0035189594
17 outcome2 ~ smoke smokeFormer -1.036036e-01 0.122196316 -8.478456e-01 0.3986112406
18 outcome2 ~ smoke  smokeNever -1.515152e-01 0.171304709 -8.844774e-01 0.3786256699

edited Dec 7, 2021 at 12:38

answered Dec 7, 2021 at 10:37

maydin

3,7533 gold badges14 silver badges31 bronze badges

Sign up to request clarification or add additional context in comments.

4 Comments

Sandro Over a year ago

very much appreciated! I am trying to adapt the ModelNo such that it reports the outcome name rather than a number. This will help when there are more than two outcomes.

Sandro Over a year ago

If I replace ModelNo=n' with ModelNo=j` it will report it with either a 1 or 2 given outcomelist is of length=2. Any way to get the actual character string of outcome1 and outcome2 in there?

maydin Over a year ago

@Sandro See my edit please.

Sandro Over a year ago

thank you again, you've been most helpful! @maydin

Collectives™ on Stack Overflow

Nested loop with simple linear regression

1 Answer 1

4 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

4 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related