用于线性回归 lm(y~x) 的 R 循环并将模型输出保存为数据集

R loop for linear regression lm(y~x) and save model output as a dataset

我想用一个 y 和几个 x 的数据集做一个回归循环 lm(y~x),并且 运行 每个 [=14] 的回归=],然后还将结果(估计值、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 值)。

我找到了一个代码,允许我 select 以“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 包“扫帚”找到了以下方法:

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

但我总是收到 group_by()do()

的错误

非常感谢您的帮助!

一种选择是使用 lapply 对每个自变量执行回归。使用 broom 库中的 tidy 将结果存储为整齐的格式。

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 

第一步

您的数据很乱,让我们整理一下。

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)))

reprex package (v0.3.0)

于 2020-07-13 创建

第二步

迭代每种类型,并使用这个漂亮的函数获得 P 值

# From 
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

reprex package (v0.3.0)

于 2020-07-13 创建