Pytanie szybki / elegancki sposób konstruowania tabeli podsumowań średniej / wariancji


Mogę osiągnąć to zadanie, ale czuję, że musi istnieć "najlepszy" (najlepszy, najbardziej kompaktowy, najczystszy kod, najszybszy?) Sposób na zrobienie tego i nie doszedłem do tej pory ...

Dla określonego zestawu czynników jakościowych chcę skonstruować tabelę średnich i wariancji według grupy.

generować dane:

set.seed(1001)
d <- expand.grid(f1=LETTERS[1:3],f2=letters[1:3],
                 f3=factor(as.character(as.roman(1:3))),rep=1:4)
d$y <- runif(nrow(d))
d$z <- rnorm(nrow(d))

Pożądane wyjście:

  f1 f2  f3    y.mean      y.var
1  A  a   I 0.6502307 0.09537958
2  A  a  II 0.4876630 0.11079670
3  A  a III 0.3102926 0.20280568
4  A  b   I 0.3914084 0.05869310
5  A  b  II 0.5257355 0.21863126
6  A  b III 0.3356860 0.07943314
... etc. ...

za pomocą aggregate/merge:

library(reshape)
m1 <- aggregate(y~f1*f2*f3,data=d,FUN=mean)
m2 <- aggregate(y~f1*f2*f3,data=d,FUN=var)
mvtab <- merge(rename(m1,c(y="y.mean")),
      rename(m2,c(y="y.var")))

za pomocą ddply/summarise (prawdopodobnie najlepiej, ale nie udało się, aby to działało):

mvtab2 <- ddply(subset(d,select=-c(z,rep)),
                .(f1,f2,f3),
                summarise,numcolwise(mean),numcolwise(var))

prowadzi do

Error in output[[var]][rng] <- df[[var]] : 
  incompatible types (from closure to logical) in subassignment type fix

za pomocą melt/cast (może najlepiej?)

mvtab3 <- cast(melt(subset(d,select=-c(z,rep)),
          id.vars=1:3),
     ...~.,fun.aggregate=c(mean,var))
## now have to drop "variable"
mvtab3 <- subset(mvtab3,select=-variable)
## also should rename response variables

Czy (?) Nie będzie działać reshape2. Wyjaśniając ...~. dla kogoś może być trudne!


21
2017-09-16 18:58


pochodzenie




Odpowiedzi:


Jestem trochę zdziwiony. Czy to nie działa:

mvtab2 <- ddply(d,.(f1,f2,f3),
            summarise,y.mean = mean(y),y.var = var(y))

Daje mi to coś takiego:

   f1 f2  f3    y.mean       y.var
1   A  a   I 0.6502307 0.095379578
2   A  a  II 0.4876630 0.110796695
3   A  a III 0.3102926 0.202805677
4   A  b   I 0.3914084 0.058693103
5   A  b  II 0.5257355 0.218631264

Która jest w odpowiedniej formie, ale wygląda na to, że wartości są różne od podanych.

Edytować

Oto, jak utworzyć wersję numcolwise praca:

mvtab2 <- ddply(subset(d,select=-c(z,rep)),.(f1,f2,f3),summarise,
                y.mean = numcolwise(mean)(piece),
                y.var = numcolwise(var)(piece)) 

Zapomniałeś przekazać rzeczywiste dane numcolwise. A potem jest trochę ddply trik, że każdy kawałek jest nazywany piece wewnętrznie. (Na co Hadley wskazuje w komentarzach nie należy polegać, ponieważ może się zmienić w przyszłych wersjach plyr.)


11
2017-09-16 19:09



problem z wartościami jest taki, że tak naprawdę nie zrobiłem set.seed gdzie twierdziłem - dodałem go później (Ups, poprawiono teraz, więc twój komentarz o "innym" nie ma już sensu - przepraszam). Wygląda na to, że to działa. Muszę wymyślić, co jest grane numcolwise - Myślałem, że automatycznie obsłuży to fakt, że w ramce danych pozostała tylko jedna zmienna numeryczna, ale może jestem w tym wszystkim zdezorientowana. - Ben Bolker
trudne wezwanie do zaakceptowania, ale to była pierwsza odpowiedź, a ja lubię elegancję (żadnych anonimowych funkcji do zdefiniowania), nawet jeśli jest wolniejsza. - Ben Bolker
Nie powinieneś polegać na użyciu jakichkolwiek zmiennych wewnętrznych - które zepsują się, gdy w końcu znajdę sposób na ustalenie zakresu w plyr - hadley
@hadley Bummer! Pomyślałem, że ta mała sztuczka była bardzo przydatna w wielu przypadkach. - joran
Jeśli robisz to, to prawdopodobnie znak, że nie powinieneś używać podsumowania. Muszę jednak myśleć o zaimkach. - hadley


Oto rozwiązanie z wykorzystaniem data.table

library(data.table)
d2 = data.table(d)
ans = d2[,list(avg_y = mean(y), var_y = var(y)), 'f1, f2, f3']

17
2017-09-16 19:14





(Głosowałem za Joshua.) Oto rozwiązanie Hmisc :: summary.formula. Zaletą tego jest dla mnie to, że jest dobrze zintegrowany z "kanałem" wyjścia Hmisc :: latex.

summary(y ~ interaction(f3,f2,f1), data=d, method="response", 
                    fun=function(y) c(mean.y=mean(y) ,var.y=var(y) ))
#-----output----------
y    N=108

+-----------------------+-------+---+---------+-----------+
|                       |       |N  |mean.y   |var.y      |
+-----------------------+-------+---+---------+-----------+
|interaction(f3, f2, f1)|I.a.A  |  4|0.6502307|0.095379578|
|                       |II.a.A |  4|0.4876630|0.110796695|

wycięte wyjście pokazujące lateks -> PDF -> wyjście png:

enter image description here


12
2017-09-16 19:32





@joran jest spot-on z ddply odpowiedź. Oto, jak mógłbym to zrobić aggregate. Zwróć uwagę, że unikam interfejsu formuły (jest wolniejszy).

aggregate(d$y, d[,c("f1","f2","f3")], FUN=function(x) c(mean=mean(x),var=var(x)))

11
2017-09-16 19:11





Jestem lekko uzależniony od porównań prędkości, mimo że w dużej mierze nie mają one znaczenia w tej sytuacji ...

joran_ddply <- function(d) ddply(d,.(f1,f2,f3),
                                 summarise,y.mean = mean(y),y.var = var(y))
joshulrich_aggregate <- function(d) {
  aggregate(d$y, d[,c("f1","f2","f3")],
            FUN=function(x) c(mean=mean(x),var=var(x)))
}

formula_aggregate <- function(d) {
  aggregate(y~f1*f2*f3,data=d,
            FUN=function(x) c(mean=mean(x),var=var(x)))
}
library(data.table)
d2 <- data.table(d)
ramnath_datatable <- function(d) {
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1, f2, f3']
}


library(Hmisc)
dwin_hmisc <- function(d) {summary(y ~ interaction(f3,f2,f1), 
                   data=d, method="response", 
                   fun=function(y) c(mean.y=mean(y) ,var.y=var(y) ))
                         }


library(rbenchmark)
benchmark(joran_ddply(d),
          joshulrich_aggregate(d),
          ramnath_datatable(d2),
          formula_aggregate(d),
          dwin_hmisc(d))

aggregate jest najszybszy (nawet szybciej niż data.table, co jest dla mnie zaskoczeniem, chociaż rzeczy mogą się różnić w przypadku większej tabeli do agregacji), nawet przy użyciu interfejsu formuły ...)

                     test replications elapsed relative user.self sys.self
5           dwin_hmisc(d)          100   1.235 2.125645     1.168    0.044
4    formula_aggregate(d)          100   0.703 1.209983     0.656    0.036
1          joran_ddply(d)          100   3.345 5.757315     3.152    0.144
2 joshulrich_aggregate(d)          100   0.581 1.000000     0.596    0.000
3   ramnath_datatable(d2)          100   0.750 1.290878     0.708    0.000

(Teraz potrzebuję tylko Dirka, żeby przyspieszyć i opublikować Rcpp rozwiązanie, które jest 1000 razy szybsze niż cokolwiek innego ...)


7
2017-09-16 21:00



Sprawdziłem na większych stołach z 2700 rzędami i znalazłem to data.table atutem aggregate rozwiązanie oparte o 1,5x. - Ramnath
@Ramnath Jeśli tabele są większe, ulepszenie twojego rozwiązania jest 21 razy szybsze. - marbel
Nic dziwnego, biorąc pod uwagę ile data.table został zoptymalizowany. Pełny kredyt dla Matta Dowle i Aruna :) - Ramnath


Przekonałem się o tym pytaniu i stwierdziłem, że testy porównawcze są wykonywane przy użyciu małych tabel, więc trudno jest stwierdzić, która metoda jest lepsza w przypadku 100 wierszy.

Zmieniłem również trochę dane, aby uczynić je "nieposortowanymi", to byłby bardziej powszechny przypadek, na przykład, ponieważ dane znajdują się w bazie danych. Dodałem kilka prób danych z tabelą danych, aby sprawdzić, czy wcześniejsze ustawienie klucza jest szybsze. Wygląda na to, że ustawienie klawisza wcześniej nie poprawia znacznie wydajności, więc rozwiązanie ramnath wydaje się być najszybsze.

set.seed(1001)
d <- data.frame(f1 = sample(LETTERS[1:3], 30e5, replace = T), f2 = sample(letters[1:3], 30e5, replace = T),
                f3 = sample(factor(as.character(as.roman(1:3))), 30e5, replace = T), rep = sample(1:4, replace = T))

d$y <- runif(nrow(d))
d$z <- rnorm(nrow(d))

str(d)

require(Hmisc)
require(plyr)
require(data.table)
d2 = data.table(d)
d3 = data.table(d)

# Set key of d3 to compare how fast it is if the DT is already keyded
setkey(d3,f1,f2,f3)

joran_ddply <- function(d) ddply(d,.(f1,f2,f3),
                                 summarise,y.mean = mean(y),y.var = var(y))

formula_aggregate <- function(d) {
  aggregate(y~f1*f2*f3,data=d,
            FUN=function(x) c(mean=mean(x),var=var(x)))
}

ramnath_datatable <- function(d) {
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

key_agg_datatable <- function(d) {
  setkey(d2,f1,f2,f3)
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

one_key_datatable <- function(d) {
  setkey(d2,f1)
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

including_3key_datatable <- function(d) {
  d[,list(avg_y = mean(y), var_y = var(y)), 'f1,f2,f3']
}

dwin_hmisc <- function(d) {summary(y ~ interaction(f3,f2,f1), 
                                   data=d, method="response", 
                                   fun=function(y) c(mean.y=mean(y) ,var.y=var(y) ))
}

require(rbenchmark)
benchmark(joran_ddply(d),
          joshulrich_aggregate(d),
          ramnath_datatable(d2),
          including_3key_datatable(d3),
          one_key_datatable(d2),
          key_agg_datatable(d2),
          formula_aggregate(d),
          dwin_hmisc(d)
          )

#                         test replications elapsed relative user.self sys.self
#                dwin_hmisc(d)          100 1757.28  252.121   1590.89   165.65
#         formula_aggregate(d)          100  433.56   62.204    390.83    42.50
# including_3key_datatable(d3)          100    7.00    1.004      6.02     0.98
#               joran_ddply(d)          100  173.39   24.877    119.35    53.95
#      joshulrich_aggregate(d)          100  328.51   47.132    307.14    21.22
#        key_agg_datatable(d2)          100   24.62    3.532     19.13     5.50
#        one_key_datatable(d2)          100   29.66    4.255     22.28     7.34
#        ramnath_datatable(d2)          100    6.97    1.000      5.96     1.01

4
2018-01-16 14:35





Znajduję pakiet doBy ma kilka bardzo przydatnych funkcji do takich rzeczy. Na przykład funkcja ? summaryBy jest bardzo przydatny. Rozważać:

> summaryBy(y~f1+f2+f3, data=d, FUN=c(mean, var))
   f1 f2  f3    y.mean       y.var
1   A  a   I 0.6502307 0.095379578
2   A  a  II 0.4876630 0.110796695
3   A  a III 0.3102926 0.202805677
4   A  b   I 0.3914084 0.058693103
5   A  b  II 0.5257355 0.218631264
6   A  b III 0.3356860 0.079433136
7   A  c   I 0.3367841 0.079487973
8   A  c  II 0.6273320 0.041373836
9   A  c III 0.4532720 0.022779672
10  B  a   I 0.6688221 0.044184575
11  B  a  II 0.5514724 0.020359289
12  B  a III 0.6389354 0.104056229
13  B  b   I 0.5052346 0.138379070
14  B  b  II 0.3933283 0.050261804
15  B  b III 0.5953874 0.161943989
16  B  c   I 0.3490460 0.079286849
17  B  c  II 0.5534569 0.207381592
18  B  c III 0.4652424 0.187463143
19  C  a   I 0.3340988 0.004994589
20  C  a  II 0.3970315 0.126967554
21  C  a III 0.3580250 0.066769484
22  C  b   I 0.7676858 0.124945402
23  C  b  II 0.3613772 0.182689385
24  C  b III 0.4175562 0.095933470
25  C  c   I 0.3592491 0.039832864
26  C  c  II 0.7882591 0.084271963
27  C  c III 0.3936949 0.085758343

Tak więc wywołanie funkcji jest proste, łatwe w użyciu, i powiedziałbym, eleganckie.

Teraz, jeśli twoim głównym problemem jest szybkość, wydaje się, że byłoby to rozsądne - przynajmniej w przypadku mniejszych zadań (pamiętaj, że nie mogłem uzyskać ramnath_datatable funkcja działa z dowolnego powodu):

                     test replications elapsed relative user.self 
4           dwin_hmisc(d)          100    0.50    2.778      0.50 
3    formula_aggregate(d)          100    0.23    1.278      0.24 
5       gung_summaryBy(d)          100    0.34    1.889      0.35 
1          joran_ddply(d)          100    1.34    7.444      1.32 
2 joshulrich_aggregate(d)          100    0.18    1.000      0.19 

3
2018-05-30 04:33





A oto rozwiązanie wykorzystujące nowy Hadley Wickham dplyrbiblioteka.

library(dplyr)
d %>% group_by(f1, f2, f3) %>% 
summarise(y.mean = mean(y), z.mean = mean(z))

2
2018-06-10 21:59