Линейная модель и dplyr - лучшее решение?

Я получил много хороших отзывов в вопросе который я недавно спросил, и руководствовался использованием 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

И я хочу это в форме этого (значения для Slope и Pearson, смоделированные для иллюстрации):

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

Наклон, являющийся линейно-наименьшим квадратным наклоном. Теоретически script будет выглядеть так:

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 вызывает ошибку, так как не может найти Time или Temp. Он может запускать lm(data$Temp ~ data$Time)$coeff[2], но возвращает наклон для всего набора данных, а не подмножествую форму, которую я ищу. cor() похоже, отлично работает в разделе group_by, так что есть специальный синтаксис, который мне нужно передать в lm(), чтобы он выполнялся аналогичным образом или полностью использовал другую функцию, чтобы получить наклон, пройденный из подмножество?

Ответы

Ответ 1

У вас есть несколько вопросов.

  • Если вы группируете свои данные по трем переменным (или даже по 2), у вас недостаточно отдельных значений, чтобы запустить модель линейной регрессии
  • Pearson требует двух числовых значений, тогда как Time является фактором, который преобразует в числовое значение, не будет иметь большого смысла
  • Третий вопрос: вам нужно использовать 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