R l oop для линейной регрессии lm (y ~ x) и сохранения вывода модели в виде набора данных - PullRequest
Новая
       106

R l oop для линейной регрессии lm (y ~ x) и сохранения вывода модели в виде набора данных

0 голосов
/ 13 июля 2020

Я хотел бы сделать регрессию l oop lm (y ~ x) с набором данных с одним y и несколькими x, и запустить регрессию для каждого x, а затем также сохранить результаты (оценка, p-значения) в data.frame(), поэтому мне не нужно копировать их вручную (тем более, что мои реальные данные намного больше). Я думаю, что это не должно быть слишком сложно, но я изо всех сил стараюсь заставить его работать и ценю вашу помощь: Вот мой пример набора данных: 

sample_data <- data.frame(
  fit = c(0.8971963, 1.4205607, 1.4953271, 0.8971963, 1.1588785, 0.1869159, 1.1588785, 1.142857143, 0.523809524),
  Xbeta = c(2.8907744,  -0.7680777,  -0.7278847, -0.06293916, -0.04047017, 2.3755812, 1.3043990,  -0.5698354, -0.5698354),
  Xgamma = c( 0.1180758, -0.6275700, 0.3731964,  -0.2353454,-0.5761923,  -0.5186803, 0.43041835, 3.9111749, -0.5030638),
  Xalpha = c(0.2643091, 1.6663923,  0.4041057, -0.2100472, -0.2100472, 7.4874195, -0.2385278,  0.3183102, -0.2385278),
  Xdelta = c(0.1498646, -0.6325119, -0.5947564, -0.2530748, 3.8413339, 0.6839322, 0.7401834,  3.8966404,  1.2028175)
)
#yname <- ("fit")
#xnames <- c("Xbeta ","Xgamma", "Xalpha", "Xdelta")

Простая регрессия с первой независимой переменной Xbeta будет выглядеть так lm(fit~Xbeta, data= sample_data), и я хотел бы запустить регрессию для каждой переменной, начинающейся с «X», а затем сохранить результат (оценку, значение p).

Я нашел код, который позволяет я выбираю переменные, начинающиеся с "X", а затем использую их для модели, но код выдает мне ошибку начиная с mutate() (обозначено #). 

library(tidyverse)
library(tsibble)

sample_data %>% 
  gather(stock, return, starts_with("X")) %>%  
  group_nest(stock) 
#  %>% 
#  mutate(model = map(data,
#                     ~lm(formula = "fit~ return",
#                         data = .x))
# ),
#           resid = map(model, residuals)
#           ) %>%
#           unnest(c(data,resid)) %>%  
#           summarise(sd_residual = sd(resid))

Для последующего сохранения результаты регрессии я также нашел следующую оценку с использованием пакета R "broom": r для l oop для регрессии lm (y ~ x) 

sample_data%>% 
  group_by(y,x)%>%                            # get combinations of y and x to regress
  do(tidy(lm(fRS_relative~xvalue, data=.)))

Но я всегда получить ошибку для group_by() и do()

Я очень благодарен за вашу помощь!

Ответы [ 2 ]

0 голосов
/ 13 июля 2020

Шаг первый

Ваши данные беспорядочные, позвольте нам привести их в порядок. 

sample_data <- data.frame(
  fit = c(0.8971963, 1.4205607, 1.4953271, 0.8971963, 1.1588785, 0.1869159, 1.1588785, 1.142857143, 0.523809524),
  Xbeta = c(2.8907744,  -0.7680777,  -0.7278847, -0.06293916, -0.04047017, 2.3755812, 1.3043990,  -0.5698354, -0.5698354),
  Xgamma = c( 0.1180758, -0.6275700, 0.3731964,  -0.2353454,-0.5761923,  -0.5186803, 0.43041835, 3.9111749, -0.5030638),
  Xalpha = c(0.2643091, 1.6663923,  0.4041057, -0.2100472, -0.2100472, 7.4874195, -0.2385278,  0.3183102, -0.2385278),
  Xdelta = c(0.1498646, -0.6325119, -0.5947564, -0.2530748, 3.8413339, 0.6839322, 0.7401834,  3.8966404,  1.2028175)
)

tidyframe = data.frame(fit = sample_data$fit,
           X = c(sample_data$Xbeta,sample_data$Xgamma,sample_data$Xalpha,sample_data$Xdelta),
           type = c(rep("beta",9),rep("gamma",9),rep("alpha",9),rep("delta",9)))

Создано 13.07.2020 пакет репекс (v0.3.0)

Шаг второй

Перебираем каждый тип и получаем P-значение, используя этот изящный function 

# From https://stackoverflow.com/a/5587781/3212698
lmp <- function (modelobject) {
  if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
  f <- summary(modelobject)$fstatistic
  p <- pf(f[1],f[2],f[3],lower.tail=F)
  attributes(p) <- NULL
  return(p)
}

Затем сделайте несколько хитроумных конвейеров 


tidyframe %>% group_by(type) %>%
  summarise(type = type, p = lmp(lm(formula = fit ~ X))) %>%
  unique()
#> `summarise()` regrouping output by 'type' (override with `.groups` argument)
#> # A tibble: 4 x 2
#> # Groups:   type [4]
#>   type       p
#>   <fct>  <dbl>
#> 1 alpha 0.0902
#> 2 beta  0.148 
#> 3 delta 0.904 
#> 4 gamma 0.529

Создано 13 июля 2020 года пакетом репекс (v0.3.0)

0 голосов
/ 13 июля 2020

Один из вариантов - использовать lapply для выполнения регрессии с каждой из независимых переменных. Используйте tidy из библиотеки broom, чтобы сохранить результаты в удобном формате. 

lapply(1:length(xnames), 
       function(i) broom::tidy(lm(as.formula(paste0('fit ~ ', xnames[i])), data = sample_data))) -> test

, а затем объедините все результаты в один фрейм данных: 

do.call('rbind', test)


# term        estimate std.error statistic   p.value
# <chr>          <dbl>     <dbl>     <dbl>     <dbl>
#   1 (Intercept)   1.05      0.133      7.89  0.0000995
# 2 Xbeta        -0.156     0.0958    -1.62  0.148    
# 3 (Intercept)   0.968     0.147      6.57  0.000313 
# 4 Xgamma        0.0712    0.107      0.662 0.529    
# 5 (Intercept)   1.09      0.131      8.34  0.0000697
# 6 Xalpha       -0.0999    0.0508    -1.96  0.0902   
# 7 (Intercept)   0.998     0.175      5.72  0.000723 
# 8 Xdelta       -0.0114    0.0909    -0.125 0.904
...