Получить значения p для конкретной переменной во многих моделях со всеми возможными комбинациями других независимых переменных

Question

Я пытаюсь запустить много моделей регрессии со всеми возможными комбинациями набора независимых переменных.

В этом примере меня интересуют коэффициенты cyl со всеми возможными комбинациями других переменных, перечисленных в xlist.

df <- mtcars
md <- "mpg ~ cyl" 
xlist <- c("disp", "hp", "am")
n <- length(xlist)
# get a list of all possible combinations of xlist 
comb_lst <- unlist(lapply(1:n, function(x) combn(xlist, x, simplify=F)), recursive = F)
# get a list of all models
md_lst <- lapply(comb_lst, function(x) paste(md, "+", paste(x, collapse = "+")))
# run those models and obtain coefficients for cyl
coefs <- unlist(lapply(md_lst, function(x) lm(as.formula(x),data=df)$coe[2]))

Хорошо работает, чтобы получить все коэффициенты для cyl. Однако я не знаю, как получить значение р, соответствующее каждому из этих коэффициентов.

pvalues <- lapply(md, function(x) lm(as.formula(x),data=df)$?[2]))

Любые предложения будут оценены.

Answer 1

Простой подход для каждого коэффициента cyl с без пакета:

pvalues <- lapply(md_lst, function(x) summary(lm(as.formula(x),data=df))$coefficients[,4])

[[1]]
 (Intercept)          cyl         disp 
4.022869e-14 3.366495e-02 5.418572e-02 

[[2]]
 (Intercept)          cyl           hp 
1.620660e-16 4.803752e-04 2.125285e-01 

[[3]]
 (Intercept)          cyl           am 
7.694408e-14 1.284560e-07 5.635445e-02 

[[4]]
 (Intercept)          cyl         disp           hp 
1.537198e-13 1.349044e-01 8.092901e-02 3.249519e-01 

[[5]]
 (Intercept)          cyl         disp           am 
2.026114e-12 2.823412e-02 1.544849e-01 1.610559e-01 

[[6]]
 (Intercept)          cyl           hp           am 
9.270924e-12 8.635578e-02 1.692706e-02 5.464020e-03 

[[7]]
 (Intercept)          cyl         disp           hp           am 
3.724625e-11 2.800850e-01 4.760672e-01 4.416647e-02 2.520516e-02 

Использование broom::glance для всех:

pvalues <- lapply(md_lst, function(x) glance(summary(lm(as.formula(x),data=df)))$p.value)

[[1]]
[1] 1.057904e-09

[[2]]
[1] 3.161781e-09

[[3]]
[1] 1.093687e-09

[[4]]
[1] 5.053802e-09

[[5]]
[1] 3.060153e-09

[[6]]
[1] 4.790959e-10

[[7]]
[1] 2.540038e-09

Answer 2

Отказ от ответственности :: В этом ответе используется версия для разработчиков manymodelr , написанная мной. Затем вы можете выбрать только те переменные, которые вас интересуют.

Map(function(x) manymodelr::extract_model_info(x,"p_value"),
  lapply(md_lst, function(x) do.call(lm, list(formula = x, data=mtcars))))

# Just cyl
Map(function(x) manymodelr::extract_model_info(x,"p_value")["cyl"],
  lapply(md_lst, function(x) do.call(lm, list(formula = x, data=mtcars))))

Если вы не хотите использовать пакеты:

Map(function(x) coef(summary(x))[,4]["cyl"],
  lapply(md_lst, function(x) do.call(lm, list(formula = x, data=mtcars))))

Результаты: :( первая часть)

[[1]]
 (Intercept)          cyl         disp 
4.022869e-14 3.366495e-02 5.418572e-02 

[[2]]
 (Intercept)          cyl           hp 
1.620660e-16 4.803752e-04 2.125285e-01 

[[3]]
 (Intercept)          cyl           am 
7.694408e-14 1.284560e-07 5.635445e-02 

[[4]]
 (Intercept)          cyl         disp           hp 
1.537198e-13 1.349044e-01 8.092901e-02 3.249519e-01 

[[5]]
 (Intercept)          cyl         disp           am 
2.026114e-12 2.823412e-02 1.544849e-01 1.610559e-01 

[[6]]
 (Intercept)          cyl           hp           am 
9.270924e-12 8.635578e-02 1.692706e-02 5.464020e-03 

[[7]]
 (Intercept)          cyl         disp           hp           am 
3.724625e-11 2.800850e-01 4.760672e-01 4.416647e-02 2.520516e-02