R - вычислить разницу между всеми комбинациями наблюдений внутри групп

Я хочу сравнить влияние различных доз удобрений на различные сорта сельскохозяйственных культур в разных местах. Мой набор данных похож на сгенерированный ниже:

locs <- rep(c("loc1","loc2","loc3"), length.out=180)
cults <- rep(c("cult1","cult2","cult3","cult4","cult5"), length.out=180)
doses <- rep(c("no_fert","40kg","50kg","60kg"), length.out=180)
set.seed(123); yld <- runif (3*length(unique(locs))*length(unique(cults))*length(unique(doses)), min=3, max=15)

dat <- data.frame(location=locs,
                  cultivar=cults,
                  fert_dose=doses,
                  yield=yld)

Обратите внимание, что есть три повторения каждой дозы удобрения (но в моем фактическом наборе данных их больше).

Первое, что мне нужно сделать, это рассчитать среднее значение для трех повторов каждой комбинации место-сорт-удобрение.

Я могу сделать это - возможно, не очень эффективным способом - вот так:

d1 <- d2 <- d3 <- list()
for (i in 1:length(unique(locs))){
  for (j in 1:length(unique(cults))){
    for (k in 1:length(unique(doses))){
      d1[[k]] <- data.frame(location=locs[i],
                            cultivar=cults[j],
                            fert_dose=doses[k],
                            mean_yield=mean(dat[dat$location==locs[i]&dat$cultivar==cults[j]&dat$fert_dose==doses[k],]$yield))
    }
    d2[[j]] <- do.call(rbind,d1)
  }
  d3[[i]] <- do.call(rbind,d2)
}

(mean_dat <- do.call(rbind, d3))

Далее, что мне нужно сделать, это: для каждого места найти разницу в урожайности среди всех комбинаций сортов и доз удобрений.

Например, учитывая только loc1 и cult1, ожидаемый результат будет таким:

res <- "
location cultivar dose dose_mean other_cultivar other_dose other_mean diff
loc1 cult1 no_fert 9.402345 cult1 40kg 9.251377 0.150968
loc1 cult1 no_fert 9.402345 cult1 50kg 10.764692 -1.362347
loc1 cult1 no_fert 9.402345 cult1 60kg 10.119129 -0.716784

loc1 cult1 40kg 9.251377 cult1 no_fert 9.402345 -0.150968
loc1 cult1 40kg 9.251377 cult1 50kg 10.764692 -1.513315
loc1 cult1 40kg 9.251377 cult1 60kg 10.119129 -0.867752

loc1 cult1 50kg 10.764692 cult1 no_fert 9.402345 1.362347
loc1 cult1 50kg 10.764692 cult1 40kg 9.251377 1.513315
loc1 cult1 50kg 10.764692 cult1 60kg 10.119129 0.645563

loc1 cult1 60kg 10.119129 cult1 no_fert 9.402345 0.716784
loc1 cult1 60kg 10.119129 cult1 40kg 9.251377 0.867752
loc1 cult1 60kg 10.119129 cult1 50kg 10.764692 -0.645563
"
(res <- read.table(textConnection(res), sep = " ", header=T, stringsAsFactors=F))

В этой таблице я повторяю значения выхода для каждой дозы, полученные на предыдущем этапе (таблица mean_dat), и вычисляю простую разницу между ними. Полученная таблица продолжит этот анализ, включая другие сорта в столбце other_cultivar.

Я считаю, что ожидаемая таблица выглядит не очень хорошо, но она будет использоваться для подачи интерактивной панели инструментов, и это именно тот формат, который ей нужен, так что я не думаю, что у меня есть большой выбор.

Есть ли программный способ достичь этих двух результатов всего за один шаг?

r dplyr data.table

26.10.2022 19:27

Стоит ли изучать PHP в 2023-2024 годах?

Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...

Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией

В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.

Приемы CSS-макетирования - floats и Flexbox

Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...

Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest

В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...

Концепция локализации и ее применение в приложениях React ⚡️

Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...

Пользовательский скаляр GraphQL

Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...

Перейти к ответу Данный вопрос помечен как решенный

Ответы 4

С помощью data.table вы можете выполнить следующее соединение:

library(data.table)

locs <- rep(c("loc1","loc2","loc3"), length.out=180)
cults <- rep(c("cult1","cult2","cult3","cult4","cult5"), length.out=180)
doses <- rep(c("no_fert","40kg","50kg","60kg"), length.out=180)
set.seed(123); yld <- runif (3*length(unique(locs))*length(unique(cults))*length(unique(doses)), min=3, max=15)

dat <- data.frame(location=locs,
                  cultivar=cults,
                  fert_dose=doses,
                  yield=yld)


setDT(dat)

dat[dat, .(cultivar_1 = cultivar, 
           cultivar_2 = i.cultivar,
           fert_dose_1 = fert_dose,
           fert_dose_2 = i.fert_dose,
           yield_1 = yield,
           yield_2 = i.yield,
           diff = yield - i.yield), on = "location", by = .EACHI][
           !(cultivar_1 == cultivar_2 & fert_dose_1 == fert_dose_2)][
             order(location, cultivar_1,fert_dose_1, cultivar_2, fert_dose_2)]

#>        location cultivar_1 cultivar_2 fert_dose_1 fert_dose_2   yield_1
#>     1:     loc1      cult1      cult1        40kg        50kg  4.665673
#>     2:     loc1      cult1      cult1        40kg        50kg 13.684203
#>     3:     loc1      cult1      cult1        40kg        50kg  9.404255
#>     4:     loc1      cult1      cult1        40kg        50kg  4.665673
#>     5:     loc1      cult1      cult1        40kg        50kg 13.684203
#>    ---                                                                 
#> 10256:     loc3      cult5      cult5     no_fert        60kg  8.794829
#> 10257:     loc3      cult5      cult5     no_fert        60kg  7.265345
#> 10258:     loc3      cult5      cult5     no_fert        60kg  4.829337
#> 10259:     loc3      cult5      cult5     no_fert        60kg  8.794829
#> 10260:     loc3      cult5      cult5     no_fert        60kg  7.265345
#>          yield_2        diff
#>     1: 14.556291 -9.89061803
#>     2: 14.556291 -0.87208812
#>     3: 14.556291 -5.15203544
#>     4:  4.568348  0.09732446
#>     5:  4.568348  9.11585437
#>    ---                      
#> 10256:  7.854123  0.94070539
#> 10257:  7.854123 -0.58877881
#> 10258:  9.981001 -5.15166422
#> 10259:  9.981001 -1.18617243
#> 10260:  9.981001 -2.71565663

^{Created on 2022-10-26 with reprex v2.0.2}

26.10.2022 19:53

Аккуратное решение dplyr с добавлением всех значений для каждого местоположения в новый столбец, а затем фильтрация для удаления нескольких идентичных комбинаций:

library(tidyverse)

myfunc <- function(df) {
  df %>%
    add_column(other = list(.)) %>%
    unnest(other, names_sep = "_") %>%
    filter(!(cultivar == other_cultivar & fert_dose == other_fert_dose)) %>%
    mutate(diff = yield - other_yield)
}

datmeans <- dat %>%
  group_by(location, cultivar, fert_dose) %>%
  summarise(yield = mean(yield), .groups = "drop") %>%
  group_split(location) %>%
  map(myfunc) %>%
  bind_rows()

26.10.2022 19:54

Ответ принят как подходящий

dplyr

library(dplyr)
dat %>%
  group_by(location, cultivar, dose = fert_dose) %>%
  summarize(dose_mean = mean(yield), .groups = "drop") %>%
  full_join(., ., by = "location", suffix = c("", "_other")) %>%
  filter(cultivar != cultivar_other | dose != dose_other) %>%
  mutate(diff = dose_mean - dose_mean_other)
# # A tibble: 1,140 x 8
#    location cultivar dose  dose_mean cultivar_other dose_other dose_mean_other     diff
#    <chr>    <chr>    <chr>     <dbl> <chr>          <chr>                <dbl>    <dbl>
#  1 loc1     cult1    40kg       9.25 cult1          50kg                 10.8  -1.51   
#  2 loc1     cult1    40kg       9.25 cult1          60kg                 10.1  -0.868  
#  3 loc1     cult1    40kg       9.25 cult1          no_fert               9.40 -0.151  
#  4 loc1     cult1    40kg       9.25 cult2          40kg                 10.1  -0.830  
#  5 loc1     cult1    40kg       9.25 cult2          50kg                  8.97  0.282  
#  6 loc1     cult1    40kg       9.25 cult2          60kg                  6.71  2.54   
#  7 loc1     cult1    40kg       9.25 cult2          no_fert              11.5  -2.20   
#  8 loc1     cult1    40kg       9.25 cult3          40kg                 11.9  -2.70   
#  9 loc1     cult1    40kg       9.25 cult3          50kg                  9.21  0.0421 
# 10 loc1     cult1    40kg       9.25 cult3          60kg                  9.26 -0.00416
# # ... with 1,130 more rows

Обратите внимание, что здесь выполняется внешнее соединение cultivar и dose. Мы начали с 180 рядов и закончили с 1140, это будет расти в геометрической прогрессии.

Таблица данных

library(data.table)
DT <- as.data.table(dat)[, .(dose_mean = mean(yield)), by = .(location, cultivar, dose = fert_dose)]
merge(DT, DT, by = "location", all = TRUE, suffix = c("", "_other"), allow.cartesian = TRUE
  )[(cultivar != cultivar_other | dose != dose_other),
  ][, diff := dose_mean - dose_mean_other][]
#       location cultivar    dose dose_mean cultivar_other dose_other dose_mean_other         diff
#         <char>   <char>  <char>     <num>         <char>     <char>           <num>        <num>
#    1:     loc1    cult1 no_fert  9.402345          cult4       60kg        8.508675  0.893670057
#    2:     loc1    cult1 no_fert  9.402345          cult2       50kg        8.969489  0.432856209
#    3:     loc1    cult1 no_fert  9.402345          cult5       40kg        9.345814  0.056530679
#    4:     loc1    cult1 no_fert  9.402345          cult3    no_fert       11.243009 -1.840663741
#    5:     loc1    cult1 no_fert  9.402345          cult1       60kg       10.119129 -0.716784445
#    6:     loc1    cult1 no_fert  9.402345          cult4       50kg        9.638162 -0.235817407
#    7:     loc1    cult1 no_fert  9.402345          cult2       40kg       10.081336 -0.678991009
#    8:     loc1    cult1 no_fert  9.402345          cult5    no_fert        9.405199 -0.002854273
#    9:     loc1    cult1 no_fert  9.402345          cult3       60kg        9.255537  0.146807576
#   10:     loc1    cult1 no_fert  9.402345          cult1       50kg       10.764692 -1.362347580
#   ---                                                                                           
# 1131:     loc3    cult5    60kg  8.442893          cult5       40kg        7.217206  1.225686617
# 1132:     loc3    cult5    60kg  8.442893          cult3    no_fert        8.688523 -0.245630492
# 1133:     loc3    cult5    60kg  8.442893          cult1       60kg        7.221926  1.220966527
# 1134:     loc3    cult5    60kg  8.442893          cult4       50kg        7.918912  0.523980425
# 1135:     loc3    cult5    60kg  8.442893          cult2       40kg        7.405098  1.037794838
# 1136:     loc3    cult5    60kg  8.442893          cult5    no_fert        6.963170  1.479722527
# 1137:     loc3    cult5    60kg  8.442893          cult3       60kg        8.183201  0.259691148
# 1138:     loc3    cult5    60kg  8.442893          cult1       50kg        9.444416 -1.001523464
# 1139:     loc3    cult5    60kg  8.442893          cult4       40kg       10.264777 -1.821884187
# 1140:     loc3    cult5    60kg  8.442893          cult2    no_fert        7.196217  1.246675164

Обратите внимание, что это data.table работает хорошо, но на самом деле не уменьшает объем памяти для вычислений на месте или скорость, обычно приписываемую решениям на основе data.table.

26.10.2022 19:56

Решение data.table, которое позволяет избежать более крупного объединения и последующей фильтрации. Это будет быстро и эффективно с точки зрения памяти.

dat_mean <- setDT(dat)[,.(mean_yield = mean(yield)), location:fert_dose][, doseIdx := match(fert_dose, unique(fert_dose))]

rbindlist(
  lapply(
    parse(
      text = c(
        ".(location, doseIdx > doseIdx)",
        ".(location, doseIdx < doseIdx)"
      )
    ),
    function(e) {
      dat_mean[
        dat_mean,
        .(
          location,
          cultivar1 = cultivar,
          fert_dose1 = fert_dose,
          yield1 = mean_yield,
          cultivar2 = i.cultivar,
          fert_dose2 = i.fert_dose,
          yield2 = i.mean_yield,
          diff = mean_yield - i.mean_yield
        ),
        on = eval(e)
      ]
    }
  )
)
#>      location cultivar1 fert_dose1    yield1 cultivar2 fert_dose2   yield2        diff
#>   1:     loc1     cult4       60kg  8.508675     cult1    no_fert 9.402345 -0.89367006
#>   2:     loc1     cult2       50kg  8.969489     cult1    no_fert 9.402345 -0.43285621
#>   3:     loc1     cult5       40kg  9.345814     cult1    no_fert 9.402345 -0.05653068
#>   4:     loc1     cult1       60kg 10.119129     cult1    no_fert 9.402345  0.71678444
#>   5:     loc1     cult4       50kg  9.638162     cult1    no_fert 9.402345  0.23581741
#>  ---                                                                                  
#> 926:     loc3     cult2       40kg  7.405098     cult5       60kg 8.442893 -1.03779484
#> 927:     loc3     cult5    no_fert  6.963170     cult5       60kg 8.442893 -1.47972253
#> 928:     loc3     cult1       50kg  9.444416     cult5       60kg 8.442893  1.00152346
#> 929:     loc3     cult4       40kg 10.264777     cult5       60kg 8.442893  1.82188419
#> 930:     loc3     cult2    no_fert  7.196217     cult5       60kg 8.442893 -1.24667516

26.10.2022 20:30