цикл с итерациями по двум спискам переменных для множественной регрессии в R

Question

Я хочу написать цикл в R для запуска нескольких регрессий с одной зависимой переменной и двумя списками независимых переменных (все непрерывные переменные).Модель является аддитивной, и цикл должен выполняться путем итерации по двум спискам переменных, чтобы она брала первый столбец из первого списка + первый столбец из второго списка, затем то же самое для второго столбца в двух списках и т. Д.Проблема в том, что я не могу заставить его правильно перебирать списки, вместо этого мой цикл запускает больше моделей, чем должен.

Фрейм данных, который я здесь описываю, является лишь подмножеством, которое мне фактически нужно будет выполнить 3772 раза (я работаю над выражением транскрипта RNA-seq).

Мой фрейм данных называется сухим и содержит 22 переменных (столбцы) и 87 наблюдений (строки).Столбец 1 содержит genotypeID, столбец 2:11 содержит один набор независимых переменных для итерации, столбец 12:21 содержит второй набор независимых переменных для итерации, а столбец 23 содержит мою зависимую переменную под названием FITNESS_DRY.Вот как выглядит структура:

str(dry)
'data.frame':   87 obs. of  22 variables:
$ geneID     : Factor w/ 87 levels "e10","e101","e102",..: 12 15 17 24 25    30 35 36 38 39 ...
$ RDPI_T1    : num  1.671 -0.983 -0.776 -0.345 0.313 ...
$ RDPI_T2    : num  -0.976 -0.774 -0.532 -1.137 1.602 ...
$ RDPI_T3    : num  -0.197 -0.324 0.805 -0.701 -0.566 ...
$ RDPI_T4    : num  0.289 -0.92 1.117 -1.214 -0.447 ...
$ RDPI_T5    : num  -0.671 1.963 NA -1.024 -0.295 ...
$ RDPI_T6    : num  2.606 -1.116 -0.383 -0.893 0.119 ...
$ RDPI_T7    : num  -0.843 -0.229 -0.297 0.504 -0.712 ...
$ RDPI_T8    : num  -0.227 NA NA -0.816 -0.761 ...
$ RDPI_T9    : num  0.754 -1.304 1.867 -0.514 -1.377 ...
$ RDPI_T10   : num  1.1352 -0.1028 -0.69 2.0242 -0.0925 ...
$ DRY_T1     : num  0.6636 -0.64508 -0.24643 -1.43231 -0.00855 ...
$ DRY_T2     : num  1.008 0.823 -0.658 -0.148 0.272 ...
$ DRY_T3     : num  -0.518 -0.357 1.294 0.408 0.771 ...
$ DRY_T4     : num  0.0723 0.2834 0.5198 1.6527 0.4259 ...
$ DRY_T5     : num  0.1831 1.9984 NA 0.0923 0.1232 ...
$ DRY_T6     : num  -1.55 0.366 0.692 0.902 -0.993 ...
$ DRY_T7     : num  -2.483 -0.334 -1.077 -1.537 0.393 ...
$ DRY_T8     : num  0.396 NA NA -0.146 -0.468 ...
$ DRY_T9     : num  -0.694 0.353 2.384 0.665 0.937 ...
$ DRY_T10    : num  -1.24 -1.57 -1.36 -3.88 -1.4 ...
$ FITNESS_DRY: num  1.301 3.365 0.458 0.346 1.983 ...

Цель состоит в том, чтобы запустить 10 множественных регрессий, выглядящих следующим образом:

lm1<-lm(FITNESS_DRY~DRY_T1+RDPI_T1)
lm2<-lm(FITNESS_DRY~DRY_T2+RDPI_T2)

и так далее, перебирая все десять столбцов для обоих списков.эквивалентно следующему с точки зрения индексации

lm1<-lm(FITNESS_DRY~dry[,12]+dry[,2])
lm1<-lm(FITNESS_DRY~dry[,12]+dry[,2])

и т. д.

Мой цикл должен затем вычислять итоги для каждой модели и объединять все значения (4-й столбец сводки lm) ввыходной объект.

Сначала я определил свои списки переменных

var_list<-list(
var1=dry[,12:21],
var2=dry[,2:11]
)

Это цикл, который я пробовал, который не работает должным образом:

lm.test1<-name<-vector()
for (i in 12:length(var_list$var1)){
    for (j in 2:length(var_list$var2)){
lm.tmp<-lm(FITNESS_DRY~dry[,i]+dry[,j], na.action=na.omit, data=dry)
sum.tmp<-summary(lm.tmp)
lm.test1<-rbind(lm.test1,sum.tmp$coefficients[,4]) }
}

Цикл возвращает это сообщение об ошибке:

Warning message:
In rbind(lm.test6, sum.tmp$coefficients[, 4]) :
number of columns of result is not a multiple of vector length (arg 2)

Я могу вызвать объект "lm.test1", но этот объект имеет 27 строк вместо 10, которые я хочу, поэтому итерации здесь не работают должным образом.Может кто-нибудь помочь с этим, пожалуйста?Кроме того, было бы замечательно, если бы я мог добавить имена моих столбцов для каждого списка переменных в резюме.Я пытался использовать это для каждого списка переменных, но без успеха:

name<-append(name, as.character(colnames(var_list$var1))

Есть идеи?Заранее благодарен за любую помощь!

ОБНОВЛЕНИЕ1: Дополнительная информация о полном наборе данных: мои фактические данные будут по-прежнему содержать первый столбец «geneID», затем у меня будет 3772 имен столбцов DRY_T1 .... DRY_T3772, затемеще 3772 столбца с именами RDPI_T1 ... RDPI_T3772 и, наконец, моя зависимая переменная "FITNESS_DRY".Я все еще хочу запустить все аддитивные модели как таковые:

lm1<-lm(FITNESS_DRY~DRY_T1+RDPI_T1)
lm2<-lm(FITNESS_DRY~DRY_T2+RDPI_T2)
lm3772<-lm(FITNESS_DRY~DRY_T3772+RDPI_T3772)

Я имитировал набор данных как таковой:

set.seed(2)
dat3 = as.data.frame(replicate(7544, runif(20)))
names(dat3) = paste0(rep(c("DRY_T","RDPI_T"),each=3772), 1:3772)
dat3 = cbind(dat3, FITNESS_DRY=runif(20))

Затем я запускаю цикл for:

models = list()
for(i in 1:3772) {
vars = names(dat3)[grepl(paste0(i,"$"), names(dat3))]
models2[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" 
+")),
                                 data = dat3)
}

Это прекрасно работает при моделировании данных, но когда я пробую его на моем реальном наборе данных, который настроен точно так же, он не работает.Цикл, вероятно, имеет проблемы с обработкой чисел с двумя или более цифрами.Я получаю это сообщение об ошибке:

Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
0 (non-NA) cases

ОБНОВЛЕНИЕ 2: действительно, у модели были проблемы с обработкой чисел с двумя или более цифрами.Чтобы увидеть, как что-то пошло не так в оригинальной версии, я использовал это: (мой набор данных называется «dry2»):

names(dry2)[grepl("2$", names(dry2))]

Это вернуло все переменные DRY_T и RDPI_T с числами, содержащими «2» вместо одногопара DRY_T и RDPI_T.

Чтобы решить эту проблему, этот новый код работает:

models = list()

for(i in 1:3772) {
vars = names(dry2)[names(dry2) %in% paste0(c("DRY_T", "RDPI_T"), i)]
models[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" +   ")),
data = dry2)
}

Answer 1

Существует несколько способов настроить формулы модели для итерации. Вот один из подходов, который мы демонстрируем, используя цикл for или map из пакета purrr для итерации. Затем мы используем tidy из пакета broom, чтобы получить коэффициенты и значения p.

library(tidyverse)
library(broom)

# Fake data
set.seed(2)
dat = as.data.frame(replicate(20, runif(20)))
names(dat) = paste0(rep(c("DRY_T","RDPI_T"),each=10), 0:9)
dat = cbind(dat, FITNESS_DRY=runif(20))

# Generate list of models

# Using for loop
models = list()

for(i in 0:9) {

  # Get the two column names to use for this iteration of the model
  vars = names(dat)[grepl(paste0(i,"$"), names(dat))]

  # Fit the model and add results to the output list
  models[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" + ")),
                                 data = dat)
}

# Same idea using purrr::map to iterate
models = map(0:9 %>% set_names(), 
             ~ {
               vars = names(dat)[grepl(paste0(.x,"$"), names(dat))]
               form = paste("FITNESS_DRY ~ ", paste(vars, collapse=" + "))
               lm(form, data = dat)
             })

# Check first two models
models[1:2]
#> $`0`
#> 
#> Call:
#> lm(formula = form, data = dat)
#> 
#> Coefficients:
#> (Intercept)       DRY_T0      RDPI_T0  
#>      0.4543       0.3025      -0.1624  
#> 
#> 
#> $`1`
#> 
#> Call:
#> lm(formula = form, data = dat)
#> 
#> Coefficients:
#> (Intercept)       DRY_T1      RDPI_T1  
#>     0.64511     -0.33293      0.06698

# Get coefficients and p-values for each model in a single data frame
results = map_df(models, tidy, .id="run_number")

results
#> # A tibble: 30 x 6
#>    run_number term        estimate std.error statistic p.value
#>    <chr>      <chr>          <dbl>     <dbl>     <dbl>   <dbl>
#>  1 0          (Intercept)   0.454      0.153     2.96  0.00872
#>  2 0          DRY_T0        0.303      0.197     1.53  0.143  
#>  3 0          RDPI_T0      -0.162      0.186    -0.873 0.395  
#>  4 1          (Intercept)   0.645      0.185     3.49  0.00279
#>  5 1          DRY_T1       -0.333      0.204    -1.63  0.122  
#>  6 1          RDPI_T1       0.0670     0.236     0.284 0.780  
#>  7 2          (Intercept)   0.290      0.147     1.97  0.0650 
#>  8 2          DRY_T2        0.270      0.176     1.53  0.144  
#>  9 2          RDPI_T2       0.180      0.185     0.972 0.345  
#> 10 3          (Intercept)   0.273      0.187     1.46  0.162  
#> # … with 20 more rows

^{Создано в 2019-06-28 пакетом Представление (v0.2.1)}

Если вам не нужно сохранять объекты модели, вы можете просто вернуть фрейм данных с коэффициентами и p-значениями:

results = map_df(0:9 %>% set_names(), 
            ~ {
              vars = names(dat)[grepl(paste0(.x,"$"), names(dat))]
              form = paste("FITNESS_DRY ~ ", paste(vars, collapse=" + "))
              tidy(lm(form, data = dat))
            }, .id="run_number")

ОБНОВЛЕНИЕ: В ответ на ваш комментарий, если вы замените все экземпляры 0:9 на 1:10 (извините, не заметил, что суффиксы вашего столбца пошли с 1:10, а не с 0: 9) и все экземпляры dat (мои поддельные данные) с dry2 (или любым другим именем, которое вы используете для своего фрейма данных), код будет работать с вашими данными, если имена столбцов совпадают как те, которые вы использовали в своем вопросе. Если вы используете разные имена столбцов, вам нужно адаптировать код, либо жестко запрограммировав новые имена, либо создав функцию, которая может принимать любые имена столбцов, которые вы используете для моделей, которые вы используете. порождающее.

Чтобы объяснить, что делает код: во-первых, нам нужно получить имена столбцов, которые мы хотим использовать в каждой итерации модели. Например, в цикле for:

vars = names(dry2)[grepl(paste0(i,"$"), names(dry2))]

Когда, например, i=2, это разрешается в:

vars = names(dry2)[grepl("2$", names(dry2))]
vars

[1] "RDPI_T2" "DRY_T2"

Итак, это два столбца, которые мы хотим использовать для генерации формулы регрессии. "2$" - это регулярное выражение (регулярные выражения - это язык сопоставления строк), которое означает: сопоставить значения в names(dry2), заканчивающиеся числом '2'.

Для создания нашей формулы мы делаем:

paste(vars, collapse=" + ")

[1] "RDPI_T2 + DRY_T2"

form = paste("FITNESS_DRY ~ ", paste(vars, collapse=" + "))
form

[1] "FITNESS_DRY ~  RDPI_T2 + DRY_T2"

И теперь у нас есть формула регрессии, которую мы используем внутри lm.

Каждая итерация (с for или map или, по предложению @ RomanLuštrik, mapply), генерирует последовательные модели.

ОБНОВЛЕНИЕ 2: Как я заметил в комментарии, я понял, что регулярное выражение paste(i, "$") не будет выполнено (при сопоставлении более одного столбца независимой переменной каждого типа), когда окончательное число больше чем одна цифра. Итак, попробуйте это вместо (и аналогично для map версии):

models = list()

for(i in 1:3772) {

  # Get the two column names to use for this iteration of the model
  vars = names(dry2)[names(dry2) %in% paste0(c("DRY_T", "RDPI_T"), i)]

  # Fit the model and add results to the output list
  models[[as.character(i)]] = lm(paste("FITNESS_DRY ~ ", paste(vars, collapse=" + ")),
                                 data = dry2)
}

Чтобы увидеть, как что-то пошло не так в оригинальной версии, запустите, например, names(dry2)[grepl("2$", names(dry2))]

Answer 2

Рассмотрите возможность преобразования вашего очень широкого фрейма данных в длинный формат с помощью reshape, который обычно является предпочтительным форматом данных практически любого приложения для обработки данных.

Для ваших нужд требуется две формы для каждой метрики _T. После изменения формы создайте индикатор T_NUM (т. Е. Снимая число DRY_T## и RDPI_T##) и используйте его вместе с соответствующими FITNESS_DRY до merge двумя метриками.

Наконец, используйте by, чтобы нарезать большой фрейм данных с помощью группировок T_NUM , чтобы создать список моделей. Ниже используется dat3 , который вы смоделировали выше. В целом, все с базой R: reshape -> TNUM <- ... -> merge -> by -> lm. Другие методы, lapply, within и Reduce, являются помощниками для кода DRY-er.

# TWO DATA FRAMES OF FOUR COLUMNS
df_list <- lapply(c("DRY_T", "RDPI_T"), function(i)
  within(reshape(dat3[c(grep(i, names(dat3)), ncol(dat3))],
                 varying = list(names(dat3)[grep(i, names(dat3))]),
                 v.names = i,
                 times = names(dat3)[grep(i, names(dat3))],
                 timevar = "T_NUM",
                 direction = "long"), {
           T_NUM <- as.integer(gsub(i, "", as.character(T_NUM)))
           id <- NULL
  })
)

# MERGE BOTH DFs
long_df <- Reduce(function(x, y) merge(x, y, by=c("T_NUM", "FITNESS_DRY")), df_list)

head(long_df, 10)
#    T_NUM FITNESS_DRY     DRY_T     RDPI_T
# 1      1   0.1528837 0.9438393 0.87948274
# 2      1   0.1925344 0.7023740 0.65120186
# 3      1   0.2193480 0.2388948 0.29875871
# 4      1   0.2743660 0.1291590 0.60097630
# 5      1   0.2877732 0.9763985 0.66921847
# 6      1   0.3082835 0.7605133 0.22456361
# 7      1   0.5196165 0.1848823 0.79543965
# 8      1   0.5603618 0.1680519 0.08759412
# 9      1   0.5789254 0.8535485 0.37942053
# 10     1   0.6291315 0.5526741 0.43043940

# NAMED LIST OF 3,772 MODELS
model_list <- by(long_df, long_df$T_NUM, function(sub) 
                  lm(FITNESS_DRY ~ DRY_T + RDPI_T, sub))

выход

summary(model_list$`1`)$coefficients
#               Estimate Std. Error    t value     Pr(>|t|)
# (Intercept)  0.7085512  0.1415849  5.0044269 0.0001085681
# DRY_T       -0.1423601  0.1985256 -0.7170867 0.4830577281
# RDPI_T      -0.1273237  0.2179249 -0.5842551 0.5667218157

summary(model_list$`2`)$coefficients
#              Estimate Std. Error   t value   Pr(>|t|)
# (Intercept) 0.3907525  0.1524423 2.5632809 0.02015115
# DRY_T       0.1952963  0.1990449 0.9811672 0.34026853
# RDPI_T      0.1979513  0.1884085 1.0506492 0.30812662

summary(model_list$`3`)$coefficients
#               Estimate Std. Error  t value   Pr(>|t|)
# (Intercept) 0.38836708  0.2076638 1.870172 0.07878049
# DRY_T       0.06995811  0.1965336 0.355960 0.72624947
# RDPI_T      0.27144752  0.2115787 1.282962 0.21672143

...