线性模型和 dplyr - 更好的解决方案？

Question

我在question I recently asked 上得到了很多很好的反馈，并被指导使用 dplyr 来转换一些数据。我遇到了 lm() 的问题，并试图从这些转换后的数据中找到一个斜率，并认为我会提出一个新问题。

首先我有如下数据：

Var1    Var2    Var3    Time           Temp
a       w       j       9/9/2014       20
a       w       j       9/9/2014       15
a       w       k       9/20/2014       10
a       w       j       9/10/2014       0
b       x       L       9/12/2014       30
b       x       L       9/12/2014       10
b       y       k       9/13/2014       20
b       y       k       9/13/2014       15
c       z       j       9/14/2014       20
c       z       j       9/14/2014       10
c       z       k       9/14/2014       11
c       w       l       9/10/2014       45
a       d       j       9/22/2014       20
a       d       k       9/15/2014       4
a       d       l       9/15/2014       23
a       d       k       9/15/2014       11

我希望它采用这种形式（模拟斜率和皮尔逊的值以进行说明）：

V1  V2  V3  Slope   Pearson
a   w   j   -3      -0.9
a   w   k   2       0
a   d   j   1.5     0.6
a   d   k   0       0.5
a   d   l   -0.5    -0.6
b   x   L   12      0.7
b   y   k   4       0.6
c   z   j   -1      -0.5
c   z   k   -3      -0.4
c   w   l   -10     -0.9

斜率是线性最小二乘斜率。理论上，脚本应该是这样的：

library(dplyr)

data <- read.table("clipboard",sep="\t",quote="",header=T)

newdata = summarise(group_by(data
                              ,Var1
                              ,Var2
                              ,Var3                            
                              )
                     ,Slope = lm(Temp ~ Time)$coeff[2]                 
                     ,Pearson = cor(Time, Temp, method="pearson")
                     )

但是 R 会抛出一个错误，比如找不到时间或温度。它可以运行lm(data$Temp ~ data$Time)$coeff[2]，但返回整个数据集的斜率，而不是我正在寻找的子集形式。 cor() 似乎在group_by 部分运行得很好，所以我需要传递给lm() 的特定语法以使其以类似的方式运行或完全使用不同的函数来获得从子集？

Answer 1

这里有几个问题。

1. 如果您将数据按 3 个（甚至 2 个）变量分组，则您没有足够的不同值来运行线性回归模型
2. Pearson 需要两个数值，而 Time 是一个转换为数值没有多大意义的因素
3. 这里的第三个问题是您需要使用do 才能运行您的线性模型

这是仅在 V1 上分组的示例

data %>%
  group_by(Var1) %>% # You can add here additional grouping variables if your real data set enables it
  do(mod = lm(Temp ~ Time, data = .)) %>%
  mutate(Slope = summary(mod)$coeff[2]) %>%
  select(-mod)
# Source: local data frame [3 x 2]
# Groups: <by row>
#   
#   Var1     Slope
# 1    a  12.66667
# 2    b  -2.50000
# 3    c -31.33333 

---

如果你确实有两个数值变量，你也可以使用do 来计算相关性，例如（我将创建一些虚拟数值变量来说明）

data %>%
  mutate(test1 = sample(1:3, n(), replace = TRUE), # Creating some numeric variables
         test2 = sample(1:3, n(), replace = TRUE)) %>%
  group_by(Var1) %>%
  do(mod = lm(Temp ~ Time, data = .),
     mod2 = cor(.$test1, .$test2, method = "pearson")) %>%
  mutate(Slope = summary(mod)$coeff[2],
         Pearson = mod2[1]) %>%
  select(-mod, -mod2)


# Source: local data frame [3 x 3]
# Groups: <by row>
#   
#   Var1     Slope     Pearson
# 1    a  12.66667  0.25264558
# 2    b  -2.50000 -0.09090909
# 3    c -31.33333  0.30151134

---

奖励解决方案：您也可以使用 data.table 包非常有效/轻松地做到这一点

library(data.table)
setDT(data)[, list(Slope = summary(lm(Temp ~ Time))$coeff[2]), Var1]
#    Var1     Slope
# 1:    a  12.66667
# 2:    b  -2.50000
# 3:    c -31.33333

或者如果我们也想创建一些虚拟变量

library(data.table)
setDT(data)[, `:=`(test1 = sample(1:3, .N, replace = TRUE), 
                   test2 = sample(1:3, .N, replace = TRUE))][, 
                   list(Slope = summary(lm(Temp ~ Time))$coeff[2],
                        Pearson = cor(test1, test2, method = "pearson")), Var1]
#    Var1     Slope     Pearson
# 1:    a  12.66667 -0.02159168
# 2:    b  -2.50000 -0.81649658
# 3:    c -31.33333 -1.00000000