R按组将回归系数存储在数据框中答案

【问题标题】：R Storing regression coefficients in data frame column by groupR按组将回归系数存储在数据框中
【发布时间】：2019-09-28 15:32:35
【问题描述】：

我有一个包含调查结果的数据框。结果以垂直格式存储。数据框是这样的 -

set.seed(1000)

df = data.frame(RESP_ID=c(rep(1,6),rep(2,8),rep(3,9),rep(4,5),rep(5,4),rep(6,10),rep(7,4),rep(8,8),rep(9,8),rep(10,10)),
                CLIENT=c(rep("A",6),rep("A",8),rep("A",9),rep("A",5),rep("A",4),rep("B",10),rep("B",4),rep("B",8),rep("B",8),rep("B",10)),
                QST=c(paste0("Q",c(1:6)),paste0("Q",c(1:8)),paste0("Q",c(1:9)),paste0("Q",c(1:5)),paste0("Q",c(1:4)),paste0("Q",c(1:10)),paste0("Q",c(1:4)),paste0("Q",c(1:8)),paste0("Q",c(1:8)),paste0("Q",c(1:10))),
                VALUE=round(runif(72,1,4),0))

数据框说明

RESP_ID = 受访者 ID。每个 ID 对应一个受访者。在这个示例数据框中，我们有 10 位受访者。

CLIENT = 接受调查的客户姓名的通讯员。在这个示例数据框中，我们有两个客户端（A 和 B）。

QST = 对应于调查中的问题编号。

VALUE = 对应于问题的答案选项。所有问题都有 4 个答案选项（1 到 4）。

目标

对于每个客户和问题组合，我想创建一个单独的列，用于存储该问题的回归系数，该回归系数在 QST 列中回归到 Q2。

所以在回归模型中，Q2 是因变量，其他所有问题都是自变量。

我的尝试

我的尝试没有给我想要的结果。

slopesdf = df %>%
  spread(QST, VALUE, fill = 0) %>%
  group_by(CLIENT) %>%
  mutate(COEFFICIENT=lm(Q2 ~ .))

我正在尝试首先按CLIENT 和QST 分组，然后找到每个问题的斜率与 Q2 回归。我确信有更好的方法来做到这一点。

目前，我的尝试给了我以下错误消息 -

mutate_impl(.data, dots) 中的错误：评估错误：'.'丹斯拉公式和参数“数据”

期望的输出

我希望最终数据框包含三列：一列用于 CLIENT，一列用于 QST，第三列称为 COEFFICIENT，其中 CLIENT 和 QST 的每个组合的系数以 Q2 作为响应变量进行回归.

对此的任何帮助将不胜感激。

【问题讨论】：

您正在使用Q2 这不是一个列
谢谢。您有什么方法可以推荐来获得所需的输出吗？
我不清楚你在问什么。 Q2 既不是列也不是数值（线性模型的响应变量必须是数值）。 QST 是一列，但它是分类的，所以它不能是回归中的依赖（响应）变量：我可以想象你想要每个 CLIENT 的 VALUE~QST （这实际上是一个方差分析，因为QST 是分类的），但 VALUE~Q2 没有意义，因为预测变量只有一个值......你能告诉我们一个你想要的特定回归的输出吗？手工完成？
你能运行这个模型并验证这是你真正想做的吗？这是该数据的正确模型吗？
试试df %>% spread(QST, VALUE, fill = 0) %>% group_by(CLIENT) %>% nest %>% mutate(data = map(data, ~ .x %>% summarise(out = list(lm(Q2 ~ ., data = .x) %>% tidy)))) %>% unnest %>% unnest

标签： r dplyr regression

【解决方案1】：

对于这样的任务，我喜欢R for Data Science 的“多模型”方法。它符合 tidyverse 风格，使用嵌套数据框和 purrr::map 创建模型列表列。然后broom::tidy 提供实用程序来提取您需要的有关模型的信息。

我删除了ID列只是为了在数据传播后让它不碍事，并按CLIENT分组和嵌套：

library(tidyverse)

df %>%
  spread(key = QST, value = VALUE, fill = 0) %>%
  select(-RESP_ID) %>%
  group_by(CLIENT) %>%
  nest()
#> # A tibble: 2 x 2
#>   CLIENT data             
#>   <fct>  <list>           
#> 1 A      <tibble [5 × 10]>
#> 2 B      <tibble [5 × 10]>

之后，创建一列线性模型。将quick = T 传递给broom::tidy 返回模型诊断表的简化版本；如果不进行设置，您还将获得模型中每个变量的标准误差、检验统计量和 p 值。

df %>%
  spread(key = QST, value = VALUE, fill = 0) %>%
  select(-RESP_ID) %>%
  group_by(CLIENT) %>%
  nest() %>%
  mutate(lm_mod = map(data, function(d) lm(Q2 ~ ., data = d))) %>%
  mutate(mod_tidy = map(lm_mod, broom::tidy, quick = T)) %>%
  unnest(mod_tidy) %>%
  head()
#> # A tibble: 6 x 3
#>   CLIENT term        estimate
#>   <fct>  <chr>          <dbl>
#> 1 A      (Intercept)    2.67 
#> 2 A      Q1             0.333
#> 3 A      Q10           NA    
#> 4 A      Q3            -0.333
#> 5 A      Q4            -1.   
#> 6 A      Q5             1.

【讨论】：

【解决方案2】：

我不能 100% 确定这个输出是你所追求的，但是，这是在正确的轨道上吗？

df2 <- df %>%
  spread(QST, VALUE, fill = 0) %>%
  split(.$CLIENT) %>%
  lapply(., function(x) { lm(Q2 ~ ., x[, -c(1,2)])$coefficients }) %>%
  do.call(rbind, .) %>%
  data.frame(.) %>%
  mutate(CLIENT = rownames(.)) %>%
  gather(QST, COEFFICIENT, -CLIENT) %>%
  arrange(CLIENT)


> df2
   CLIENT          QST   COEFFICIENT
1       A X.Intercept. -1.200000e+01
2       A           Q1  1.000000e+00
3       A          Q10            NA
4       A           Q3  2.000000e+00
5       A           Q4  3.000000e+00
6       A           Q5  5.000000e-01
7       A           Q6            NA
8       A           Q7            NA
9       A           Q8            NA
10      A           Q9            NA
11      B X.Intercept.  5.000000e+00
12      B           Q1 -1.326970e-16
13      B          Q10  1.666667e+00
14      B           Q3  3.726559e-15
15      B           Q4 -2.000000e+00
16      B           Q5            NA
17      B           Q6            NA
18      B           Q7            NA
19      B           Q8            NA
20      B           Q9            NA

编辑：

运行拆分组件只会为每个客户端生成一个宽格式数据帧列表：

df %>%
  spread(QST, VALUE, fill = 0) %>%
  split(.$CLIENT) 

$A
  RESP_ID CLIENT Q1 Q10 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9
1       1      A  4   0  1  4  3  3  2  0  0  0
2       2      A  2   0  2  2  3  2  4  4  3  0
3       3      A  2   0  2  3  3  1  2  4  2  3
4       4      A  3   0  3  4  2  1  0  0  0  0
5       5      A  3   0  4  4  3  0  0  0  0  0

$B
   RESP_ID CLIENT Q1 Q10 Q2 Q3 Q4 Q5 Q6 Q7 Q8 Q9
6        6      B  3   2  3  2  3  2  2  1  3  3
7        7      B  2   0  3  2  2  0  0  0  0  0
8        8      B  3   0  2  4  1  3  3  2  3  0
9        9      B  2   0  1  4  2  1  3  1  2  0
10      10      B  3   2  3  3  3  3  4  2  3  3

请注意，如果未回答问题，则原始数据没有值的问题将全部填零。请参阅 Ben Bolker 关于这一点的回答。

如果您现在包含代码以在每个上面运行 lm，您将直接获得系数值，其中包括上面看到的 NA 值：

> df %>%
+   spread(QST, VALUE, fill = 0) %>%
+   split(.$CLIENT) %>%
+   lapply(., function(x) { lm(Q2 ~ ., x[, -c(1,2)])$coefficients })
$A
(Intercept)          Q1         Q10          Q3          Q4          Q5          Q6          Q7          Q8          Q9 
  6.6666667   2.0000000          NA  -1.6666667  -0.6666667  -1.6666667          NA          NA          NA          NA 

$B
(Intercept)          Q1         Q10          Q3          Q4          Q5          Q6          Q7          Q8          Q9 
       13.0        -3.0        -0.5        -2.0          NA         2.0          NA          NA          NA          NA

编辑 2：

只是为了探索更完整的数据集，如果我们使用这个df：

set.seed(42)
df <-
  expand.grid(RESP_ID = 1:10,
              CLIENT = c("A", "B"),
              QST = paste("Q", 1:10, sep = "")) %>%
  mutate(VALUE = round(runif(200, 1, 4), 0))

并运行相同的代码，我们得到没有 NA 值的系数：

> df %>%
+   spread(QST, VALUE, fill = 0) %>%
+   split(.$CLIENT) %>%
+   lapply(., function(x) { lm(Q2 ~ ., x[, -c(1,2)])$coefficients }) %>%
+   do.call(rbind, .) %>%
+   data.frame(.) %>%
+   mutate(CLIENT = rownames(.)) %>%
+   gather(QST, COEFFICIENT, -CLIENT) %>%
+   arrange(CLIENT)
   CLIENT          QST COEFFICIENT
1       A X.Intercept.  6.50000000
2       A           Q1 -4.14285714
3       A           Q3  2.50000000
4       A           Q4  0.85714286
5       A           Q5  1.00000000
6       A           Q6 -0.64285714
7       A           Q7 -1.21428571
8       A           Q8 -1.85714286
9       A           Q9  2.50000000
10      A          Q10 -0.07142857
11      B X.Intercept. -4.69924812
12      B           Q1 -0.86466165
13      B           Q3  1.56390977
14      B           Q4  1.10150376
15      B           Q5 -0.86842105
16      B           Q6  0.87593985
17      B           Q7  0.57142857
18      B           Q8  0.25187970
19      B           Q9  0.79699248
20      B          Q10 -0.12781955

【讨论】：

这看起来非常接近我的想法。你知道为什么会有 NA 吗？
@Varun - 与我想的数据有关 - 如果您直接在拆分数据帧上运行 lm 函数，这些是每个 Q 的输出值 - 我将添加到我的答案中展示

【解决方案3】：

遵循我大脑中逻辑的解决方案（我们需要将Q2 作为单独的变量提供......一旦我们以这种方式重新排列数据，我们就可以运行。（NA 值肯定是由于这个小数据集的缺陷 - 预测变量没有变化的情况，因此无法估计响应......）

(df
    %>% group_by(RESP_ID,CLIENT)
    ## add a new variable for Q2
    %>% mutate(Q2=VALUE[QST=="Q2"])
    ## drop the old one
    %>% filter(QST!="Q2")
    %>% group_by(CLIENT,QST)
    ## run the regression by group; return a data frame
    %>% do(as.data.frame(rbind(coef(lm(Q2~VALUE,data=.)))))
    ## convert wide coefficients to long
    %>% tidyr::gather(coef,value,-c(CLIENT,QST))
    %>% arrange(CLIENT)
)

【讨论】：