1

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.

Sandro
  • 101
  • 7
  • You can't pass a formula like that to lm. – user2974951 Dec 07 '21 at 09:55
  • @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 Dec 07 '21 at 10:08
  • 1
    Just put the `j` outside the quotation like : `model <- lm(paste(j, "~", i[[1]]))` ` – maydin Dec 07 '21 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 Dec 07 '21 at 10:18

1 Answers1

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
maydin
  • 3,715
  • 3
  • 10
  • 27
  • 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 Dec 07 '21 at 11:26
  • 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? – Sandro Dec 07 '21 at 11:30
  • 1
    @Sandro See my edit please. – maydin Dec 07 '21 at 12:38
  • thank you again, you've been most helpful! @maydin – Sandro Dec 07 '21 at 20:11