将回归系数从多个模型打印到共享数据帧

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等).我想将Estimate值保存到此网格中,对于每一行,有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),但是当我尝试将其粘贴到一行时,我自然会得到一个混淆两个向量长度的错误,比如说这个模型中的15个是潜在的100.

我想自动化这个过程,所以当处理新数据并运行回归时,我可以运行我的代码保存新系数的数据,然后我可以将其输出到CSV(我已经得到的那部分).思考？

Answer 1

第一步是将系数组合成一个数据框,每个模型和术语组合一行.然后,您将能够将其扩展到一个表中,每个模型一行,每个术语一列.

我的扫帚包有一个有用的功能,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此不同,它返回数据框而不是矩阵,并将术语合并为列而不是rownames).您可以将此功能应用于每个模型,然后重新组合,例如使用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填充.