r - 将多个模型的回归系数打印到共享数据框

Question

这有点初级，我知道。基本上，我想将 coef 函数中的保存数据用于模型的共享数据框，这些模型都从更大的共享数据集中提取有限的可能变量。

我有 3 套 14 个模型。每个集合使用 100 个变量数据集中的 15-25 个变量，每个模型使用大约 12 个变量的混合，这些变量因模型而异。我想做的是将 14 个模型中每一个的系数保存到一个数据框中。

Coefs=data.frame(col.names = names(EST))

系数看起来像这样:

Coefficients:
                    Estimate   Std. Error  t value           Pr(>|t|)    
RT_SCORE_USER       0.2427506  0.0310486   7.818 0.0000000000000836 ***
VOD.Window..weeks.  0.0092641  0.0009985   9.278            < 2e-16 ***
PX_WK3              0.0300395  0.0098943   3.036           0.002600 ** 

换好10-15个变量。例如，PX 有 14 周(WK1、2 等)。我想将估计值保存到这个网格中，其中每一行有 100 列列出所有可能的变量。其中大部分将是 0。这个表将被导入到 excel 中，我可以简单地对每周的模型进行交叉乘法。

我的努力是弄清楚如何将不同周的所有不同系数记录到一个 data.frame 中，其中每个模型都有一个单独的行:

       PX_WK1     PXWK_2   RT_SCORE_USER  IMAVARIABLE etc.
ESTWK1   .030     0         .24            0
ESTWK2   0        .023      .44            etc
ESTWK3   0        0         etc etc etc

我了解如何使用 coef(ESTWK1)，但是当我尝试将其粘贴到一行中时，我自然会得到一个错误，混淆了两个向量的长度，例如在这个模型中的潜在 100 中的 15。

我想自动化这个过程，以便在处理新数据并运行回归时，我可以运行我的代码保存新系数的数据，然后我可以将其输出到 CSV(我有的那部分)。想法？

Answer 1

第一步是将您的系数组合到一个数据框中，每个模型和术语组合都有一行。然后，您就可以将它展开到一个表中，每个模型一行，每个术语一列。

我的 broom package有一个有用的功能，tidy将线性拟合转换为系数数据框:

fit <- lm(mpg ~ wt + disp + qsec, mtcars)
library(broom)
tidy(fit)
#          term  estimate std.error statistic p.value
# 1 (Intercept) 19.777558    5.9383    3.3305 0.00244
# 2          wt -5.034410    1.2241   -4.1127 0.00031
# 3        disp -0.000128    0.0106   -0.0121 0.99042
# 4        qsec  0.926649    0.3421    2.7087 0.01139

(请注意，与 coef 不同，它返回数据框而不是矩阵，并将术语合并为列而不是行名)。您可以将此功能应用于每个模型，然后重新组合，例如使用 plyr 的 ldply。我们使用 20 个与您的“模型”相同的模型生成一个示例:

models <- replicate(20, lm(mpg ~ wt + disp + qsec, mtcars), simplify = FALSE)
names(models) <- paste0("MODEL", 1:20)

那么我们的“整理和重组”代码将是:

all_coefs <- plyr::ldply(models, tidy, .id = "model")
head(all_coefs)
#    model        term  estimate std.error statistic p.value
# 1 MODEL1 (Intercept) 19.777558    5.9383    3.3305 0.00244
# 2 MODEL1          wt -5.034410    1.2241   -4.1127 0.00031
# 3 MODEL1        disp -0.000128    0.0106   -0.0121 0.99042
# 4 MODEL1        qsec  0.926649    0.3421    2.7087 0.01139
# 5 MODEL2 (Intercept) 19.777558    5.9383    3.3305 0.00244
# 6 MODEL2          wt -5.034410    1.2241   -4.1127 0.00031

然后，您需要删除 std.error、statistic 和 p.value 列并展开 estimate期满。这可以通过 dplyr 和 tidyr 包来完成:

library(dplyr)
library(tidyr)
results <- all_coefs %>% select(-(std.error:p.value)) %>%
    spread(term, estimate)

这产生:

     model (Intercept)      disp  qsec    wt
1   MODEL1        19.8 -0.000128 0.927 -5.03
2   MODEL2        19.8 -0.000128 0.927 -5.03
3   MODEL3        19.8 -0.000128 0.927 -5.03
4   MODEL4        19.8 -0.000128 0.927 -5.03
5   MODEL5        19.8 -0.000128 0.927 -5.03

这是您想要的输出。 (这个输出很无聊，因为所有的模型都是一样的，但你的模型可能不同)。如果某些模型具有其他模型没有的系数，则缺失值将用 NA 填充。