Updating list in nested loop in R

user_12332455 0
Viewed 205

I am re-writting an earlier question I had which was deleted for lack of details.

I would like to extract the p-value of a number of regression that I am running with a loop and have the p-values stored for each of the independant variable as a list. My issue is with storing the new lists correctly. You will find below the code, the sample data and the expected results.

Many thanks!

ind_var <- list(c("Man","Below_Master","Worker","Rural", "Less_than_26"))
depending_var <- list(c("Q1","Q2","Q3","Q4", "Q5", "Q6"))


x <- list()

for (var_p in depending_var) {
 { for (var in ind_var)  
    
x1 <- lm(paste(var, "~", var_p), data = Ech_final_nom_SD, weights = pond)
 
x[[var_p]] <- summary(x1)$coefficients[,4]    #Where shit hits the fan

  }
}

Here is the sample data:

ID Man Below_Master Worker Rural Less than 26 Q1 Q2 Q3 Q4 Q5 Q6
1 1 1 0 0 1 4 3 1 2 2 1
2 0 0 0 1 1 2 2 2 2 3 1
3 0 1 0 0 1 4 3 1 1 3 1
4 0 0 1 0 0 1 4 1 2 3 1
5 0 1 1 0 1 4 3 4 2 3 2
6 1 0 1 1 1 3 2 1 2 2 2
7 0 1 0 0 0 3 3 3 2 1 2
8 1 0 0 0 1 3 3 1 4 1 3
9 1 0 0 1 1 2 3 1 2 2 3
10 1 1 1 0 0 1 1 1 1 2 3

Expected results (with different p-values for all) :

x_Man <- list(c("0.0234","0.9894", "0.0045", "0.9892", "0.0229", "0.0254","0.0894", "0.0047", "0.2392", "0.0029"))

x_Below_Master <- list(c("0.0434","0.9854", "0.0545", "0.0092", "0.0729", "0.0254","0.0894", "0.0047", "0.2392", "0.0029"))

x_Worker <- list(c("0.0234","0.9894", "0.0045", "0.9892", "0.0229", "0.0254","0.0894", "0.0047", "0.2392", "0.0029"))

x_Rural <- list(c("0.0234","0.9894", "0.0045", "0.9892", "0.0229", "0.0254","0.0894", "0.0047", "0.2392", "0.0029"))

x_Less_than_26 <- list(c("0.0234","0.9894", "0.0045", "0.9892", "0.0229", "0.0254","0.0894", "0.0047", "0.2392", "0.0029"))
lm r for-loop list
4 Answers

This should work:

l <- lapply(independent_var, function(x) summary(lm(reformulate(depending_var, x), data = Ech_final_nom_SD))$coefficients[, 4])
names(l) <- independent_var

output

$Man
(Intercept)          Q1          Q2          Q3          Q4          Q5          Q6 
  0.7006954   0.4183649   0.4814779   0.1129608   0.7429505   0.8469599   0.3121451 

$Below_Master
(Intercept)          Q1          Q2          Q3          Q4          Q5          Q6 
  0.4269917   0.3294021   0.8975661   0.4513400   0.1453265   0.4121115   0.9905834 

$Worker
(Intercept)          Q1          Q2          Q3          Q4          Q5          Q6 
  0.9251373   0.6595554   0.9057457   0.7430144   0.9231273   0.5541379   0.6826245 

$Rural
(Intercept)          Q1          Q2          Q3          Q4          Q5          Q6 
  0.9051546   0.8906284   0.5703472   0.8245217   0.6760297   0.6800893   0.9835419 

$Less_than_26
(Intercept)          Q1          Q2          Q3          Q4          Q5          Q6 
 0.23702999  0.05804822  0.37926448  0.26210405  0.16896185  0.12015247  0.47013899

If you want the output to be a dataframe, I'd suggest to apply the same function with map_dfr from purrr:

library(purrr)
library(dplyr)
map_dfr(independent_var, 
        function(x) summary(lm(reformulate(depending_var, x), data = Ech_final_nom_SD))$coefficients[, 4]) %>% 
  mutate(ind_var = independent_var, .before = 1)

Which gives:

# A tibble: 5 × 8
  ind_var      `(Intercept)`     Q1    Q2    Q3    Q4    Q5    Q6
  <chr>                <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Man                  0.701 0.418  0.481 0.113 0.743 0.847 0.312
2 Below_Master         0.427 0.329  0.898 0.451 0.145 0.412 0.991
3 Worker               0.925 0.660  0.906 0.743 0.923 0.554 0.683
4 Rural                0.905 0.891  0.570 0.825 0.676 0.680 0.984
5 Less_than_26         0.237 0.0580 0.379 0.262 0.169 0.120 0.470

Data

Ech_final_nom_SD <- read.table(header = T, text = "ID   Man Below_Master    Worker  Rural   Less_than_26    Q1  Q2  Q3  Q4  Q5  Q6
1   1   1   0   0   1   4   3   1   2   2   1
2   0   0   0   1   1   2   2   2   2   3   1
3   0   1   0   0   1   4   3   1   1   3   1
4   0   0   1   0   0   1   4   1   2   3   1
5   0   1   1   0   1   4   3   4   2   3   2
6   1   0   1   1   1   3   2   1   2   2   2
7   0   1   0   0   0   3   3   3   2   1   2
8   1   0   0   0   1   3   3   1   4   1   3
9   1   0   0   1   1   2   3   1   2   2   3
10  1   1   1   0   0   1   1   1   1   2   3")

independent_var <- c("Man", "Below_Master", "Worker", "Rural", "Less_than_26")
depending_var <- c("Q1","Q2","Q3","Q4", "Q5", "Q6")
user_134606020

Here is a way with the nested for loop in the question. The results object is of a different shape, it's a data.frame.

ind_var <- c("Man","Below_Master","Worker","Rural", "Less_than_26")
depending_var <- c("Q1","Q2","Q3","Q4", "Q5", "Q6")

n_rows <- length(depending_var) * length(ind_var)
x <- data.frame(
  dependent = character(n_rows),
  independent = character(n_rows),
  pval_intercept = numeric(n_rows),
  pval_indep = numeric(n_rows)
)

i <- 0L
for (var_p in depending_var) {
  for (var in ind_var) {
    i <- i + 1L
    fmla <- reformulate(var, response = var_p)
    x1 <- lm(fmla, data = Ech_final_nom_SD, weights = pond)
    x[i, "dependent"] <- var_p
    x[i, "independent"] <- var
    x[i, "pval_intercept"] <- summary(x1)$coefficients[1, 4]
    x[i, "pval_indep"] <- summary(x1)$coefficients[2, 4]
  }
}
x
#>    dependent  independent pval_intercept pval_indep
#> 1         Q1          Man   9.160936e-04 0.80277173
#> 2         Q1 Below_Master   2.026789e-03 0.18690481
#> 3         Q1       Worker   2.231686e-04 0.34553732
#> 4         Q1        Rural   2.331972e-04 0.54473730
#> 5         Q1 Less_than_26   1.768920e-02 0.05845715
#> 6         Q2          Man   3.287724e-05 0.27313922
#> 7         Q2 Below_Master   8.982180e-05 0.72446583
#> 8         Q2       Worker   3.896158e-05 0.56210550
#> 9         Q2        Rural   1.712317e-05 0.38768132
#> 10        Q2 Less_than_26   7.359017e-04 0.93892908
#> 11        Q3          Man   6.976609e-04 0.07359019
#> 12        Q3 Below_Master   3.373039e-02 0.26225475
#> 13        Q3       Worker   1.177118e-02 0.74104616
#> 14        Q3        Rural   3.748440e-03 0.63631632
#> 15        Q3 Less_than_26   3.502726e-02 0.90655491
#> 16        Q4          Man   1.337132e-03 0.47136168
#> 17        Q4 Below_Master   8.923744e-05 0.12647490
#> 18        Q4       Worker   2.199469e-04 0.46193445
#> 19        Q4        Rural   2.862687e-04 1.00000000
#> 20        Q4 Less_than_26   8.392753e-03 0.43036607
#> 21        Q5          Man   3.584230e-05 0.11143429
#> 22        Q5 Below_Master   3.701074e-04 1.00000000
#> 23        Q5       Worker   2.603420e-04 0.35588376
#> 24        Q5        Rural   1.346479e-04 0.74828807
#> 25        Q5 Less_than_26   2.975062e-03 0.62876991
#> 26        Q6          Man   2.911220e-03 0.06558768
#> 27        Q6 Below_Master   1.270703e-03 0.74043945
#> 28        Q6       Worker   1.257419e-03 0.78704595
#> 29        Q6        Rural   7.220649e-04 0.82866767
#> 30        Q6 Less_than_26   5.691830e-03 0.82866767

Created on 2022-08-22 by the reprex package (v2.0.1)

And a solution with a double lapply loop.

res <- lapply(depending_var, \(y) {
  lapply(ind_var, \(x) {
    fit <- lm(reformulate(x, y), data = Ech_final_nom_SD, weights = pond)
    data.frame(
      dependent = y,
      independent = x,
      pval_intercept = summary(fit)$coefficients[1, 4],
      pval_indep = summary(fit)$coefficients[2, 4]
    )
  })
})
res <- do.call(rbind, unlist(res, recursive = FALSE))
res
#>    dependent  independent pval_intercept pval_indep
#> 1         Q1          Man   9.160936e-04 0.80277173
#> 2         Q1 Below_Master   2.026789e-03 0.18690481
#> 3         Q1       Worker   2.231686e-04 0.34553732
#> 4         Q1        Rural   2.331972e-04 0.54473730
#> 5         Q1 Less_than_26   1.768920e-02 0.05845715
#> 6         Q2          Man   3.287724e-05 0.27313922
#> 7         Q2 Below_Master   8.982180e-05 0.72446583
#> 8         Q2       Worker   3.896158e-05 0.56210550
#> 9         Q2        Rural   1.712317e-05 0.38768132
#> 10        Q2 Less_than_26   7.359017e-04 0.93892908
#> 11        Q3          Man   6.976609e-04 0.07359019
#> 12        Q3 Below_Master   3.373039e-02 0.26225475
#> 13        Q3       Worker   1.177118e-02 0.74104616
#> 14        Q3        Rural   3.748440e-03 0.63631632
#> 15        Q3 Less_than_26   3.502726e-02 0.90655491
#> 16        Q4          Man   1.337132e-03 0.47136168
#> 17        Q4 Below_Master   8.923744e-05 0.12647490
#> 18        Q4       Worker   2.199469e-04 0.46193445
#> 19        Q4        Rural   2.862687e-04 1.00000000
#> 20        Q4 Less_than_26   8.392753e-03 0.43036607
#> 21        Q5          Man   3.584230e-05 0.11143429
#> 22        Q5 Below_Master   3.701074e-04 1.00000000
#> 23        Q5       Worker   2.603420e-04 0.35588376
#> 24        Q5        Rural   1.346479e-04 0.74828807
#> 25        Q5 Less_than_26   2.975062e-03 0.62876991
#> 26        Q6          Man   2.911220e-03 0.06558768
#> 27        Q6 Below_Master   1.270703e-03 0.74043945
#> 28        Q6       Worker   1.257419e-03 0.78704595
#> 29        Q6        Rural   7.220649e-04 0.82866767
#> 30        Q6 Less_than_26   5.691830e-03 0.82866767

Created on 2022-08-22 by the reprex package (v2.0.1)

But the simplest solution is to remember that lm can fit multiple regression models and fit all responses at once.

  resp <- as.matrix(Ech_final_nom_SD[depending_var])
  out <- lapply(ind_var, \(x) {
    fit <- lm(resp ~ get(x), data = Ech_final_nom_SD)
    cfs <- sapply(summary(fit), \(s) s$coefficients[, 4])  
    rownames(cfs)[2] <- x
    cfs
  })
  do.call(rbind.data.frame, out)
#>               Response Q1  Response Q2  Response Q3  Response Q4  Response Q5  Response Q6
#> (Intercept)  0.0009160936 3.287724e-05 0.0006976609 1.337132e-03 0.0000358423 0.0029112201
#> Man          0.8027717312 2.731392e-01 0.0735901860 4.713617e-01 0.1114342894 0.0655876769
#> (Intercept)1 0.0020267889 8.982180e-05 0.0337303906 8.923744e-05 0.0003701074 0.0012707033
#> Below_Master 0.1869048103 7.244658e-01 0.2622547504 1.264749e-01 1.0000000000 0.7404394537
#> (Intercept)2 0.0002231686 3.896158e-05 0.0117711783 2.199469e-04 0.0002603420 0.0012574186
#> Worker       0.3455373232 5.621055e-01 0.7410461602 4.619345e-01 0.3558837647 0.7870459482
#> (Intercept)3 0.0002331972 1.712317e-05 0.0037484405 2.862687e-04 0.0001346479 0.0007220649
#> Rural        0.5447373008 3.876813e-01 0.6363163232 1.000000e+00 0.7482880706 0.8286676662
#> (Intercept)4 0.0176892034 7.359017e-04 0.0350272586 8.392753e-03 0.0029750616 0.0056918299
#> Less_than_26 0.0584571528 9.389291e-01 0.9065549125 4.303661e-01 0.6287699080 0.8286676662

Created on 2022-08-22 by the reprex package (v2.0.1)

Data

x<-'ID  Man     Below_Master    Worker  Rural   Less_than_26    Q1  Q2  Q3  Q4  Q5  Q6
1   1   1   0   0   1   4   3   1   2   2   1
2   0   0   0   1   1   2   2   2   2   3   1
3   0   1   0   0   1   4   3   1   1   3   1
4   0   0   1   0   0   1   4   1   2   3   1
5   0   1   1   0   1   4   3   4   2   3   2
6   1   0   1   1   1   3   2   1   2   2   2
7   0   1   0   0   0   3   3   3   2   1   2
8   1   0   0   0   1   3   3   1   4   1   3
9   1   0   0   1   1   2   3   1   2   2   3
10  1   1   1   0   0   1   1   1   1   2   3'
Ech_final_nom_SD <- read.table(textConnection(x), header = TRUE)
pond <- rep(1, nrow(Ech_final_nom_SD))

Created on 2022-08-22 by the reprex package (v2.0.1)

user_82454060

Do you expect something like this?

Ech_final_nom_SD <- read.table(header = T, text = "
ID  Man     Below_Master    Worker  Rural   Less_than_26    Q1  Q2  Q3  Q4  Q5  Q6
1   1   1   0   0   1   4   3   1   2   2   1
2   0   0   0   1   1   2   2   2   2   3   1
3   0   1   0   0   1   4   3   1   1   3   1
4   0   0   1   0   0   1   4   1   2   3   1
5   0   1   1   0   1   4   3   4   2   3   2
6   1   0   1   1   1   3   2   1   2   2   2
7   0   1   0   0   0   3   3   3   2   1   2
8   1   0   0   0   1   3   3   1   4   1   3
9   1   0   0   1   1   2   3   1   2   2   3
10  1   1   1   0   0   1   1   1   1   2   3
                   ")
ind_var <- list(c("Man","Below_Master","Worker","Rural", "Less_than_26"))
dep_var <- list(c("Q1","Q2","Q3","Q4", "Q5", "Q6"))
pond <- rep(1, nrow(Ech_final_nom_SD))

x_list <- NULL
y_vec <- NULL

for (x in ind_var[[1]]) {
  for (y in dep_var[[1]]) {
    fit <- lm(paste(y, "~", x), data = Ech_final_nom_SD, weights = pond)
    p_value <- summary(fit)$coef[2,4]
    y_vec[y] <- p_value
  }
  x_list[[x]] <- y_vec
}

With output:

$Man
        Q1         Q2         Q3         Q4         Q5         Q6 
0.80277173 0.27313922 0.07359019 0.47136168 0.11143429 0.06558768 

$Below_Master
       Q1        Q2        Q3        Q4        Q5        Q6 
0.1869048 0.7244658 0.2622548 0.1264749 1.0000000 0.7404395 

$Worker
       Q1        Q2        Q3        Q4        Q5        Q6 
0.3455373 0.5621055 0.7410462 0.4619345 0.3558838 0.7870459 

$Rural
       Q1        Q2        Q3        Q4        Q5        Q6 
0.5447373 0.3876813 0.6363163 1.0000000 0.7482881 0.8286677 

$Less_than_26
        Q1         Q2         Q3         Q4         Q5         Q6 
0.05845715 0.93892908 0.90655491 0.43036607 0.62876991 0.82866767

I kept the dep_var and ind_var as lists as you intended. I renamed some of the variables though for clarity.

user_198273990

Here's a tidyverse oriented variation :

data <- read.table(header = T, text = "ID   Man Below_Master    Worker  Rural   Less_than_26    Q1  Q2  Q3  Q4  Q5  Q6
1   1   1   0   0   1   4   3   1   2   2   1
2   0   0   0   1   1   2   2   2   2   3   1
3   0   1   0   0   1   4   3   1   1   3   1
4   0   0   1   0   0   1   4   1   2   3   1
5   0   1   1   0   1   4   3   4   2   3   2
6   1   0   1   1   1   3   2   1   2   2   2
7   0   1   0   0   0   3   3   3   2   1   2
8   1   0   0   0   1   3   3   1   4   1   3
9   1   0   0   1   1   2   3   1   2   2   3
10  1   1   1   0   0   1   1   1   1   2   3")

ind_var <- c("Man","Below_Master","Worker","Rural", "Less_than_26")
depending_var <- c("Q1","Q2","Q3","Q4", "Q5", "Q6")

library(tidyverse)
pval_row <- function(data, depending_var, ind_var) {
  coefs <- summary(lm(reformulate(depending_var, ind_var), data))$coefficients
  set_names(tibble(!!!coefs[,"Pr(>|t|)"]), c("pval_intercept", "pval_indep"))
}

expand_grid(depending_var,ind_var) %>% 
  mutate(pmap_dfr(., pval_row, data = data))
#> # A tibble: 30 × 4
#>    depending_var ind_var      pval_intercept pval_indep
#>    <chr>         <chr>                 <dbl>      <dbl>
#>  1 Q1            Man                  0.226      0.803 
#>  2 Q1            Below_Master         0.893      0.187 
#>  3 Q1            Worker               0.101      0.346 
#>  4 Q1            Rural                0.229      0.545 
#>  5 Q1            Less_than_26         0.981      0.0585
#>  6 Q2            Man                  0.0828     0.273 
#>  7 Q2            Below_Master         0.286      0.724 
#>  8 Q2            Worker               0.251      0.562 
#>  9 Q2            Rural                0.194      0.388 
#> 10 Q2            Less_than_26         0.293      0.939 
#> # … with 20 more rows

Created on 2022-09-15 by the reprex package (v2.0.1)

user_22704750
Related