Скользящее среднее для панельных данных (с несколькими деталями)

ProgramBox

Скользящее среднее для панельных данных (с несколькими деталями)

Post author:admin
Запись опубликована:26 апреля, 2023
Post category:Вопросы по программированию

#r #lag #moving-average #panel-data #tapply

#r #задержка #скользящее среднее #панель-данные #tapply

Вопрос:

Я придумал некоторый код для вычисления скользящего среднего для панельных данных (строка в данных содержит значения одного объекта за один день). Поскольку у меня было несколько более специфичных требований, код стал довольно сложным. Слишком сложно для приложения, на мой взгляд, не слишком редкого.

Вот что мне было нужно:

скользящее среднее (среднее значений за (a) предыдущие 3 дня, исключая «текущий» день, (b) вычисляется только в том случае, если в этом окне есть минимум 2 не пропущенных значения)
соблюдение структуры панели

Не слишком сложно, не так ли?

Для 1. Я решил использовать rollapplyr() и mean( , na.rm = T) , чтобы исключить текущий день (а) я решил использовать самодельную функцию задержки и (б) оператор if. И для 2. Я заключил все в tapply() (с unlist() ), чтобы учесть структуру панели.

Вот пример кода:

 library(zoo)

# example data (with missings)
set.seed(1)
df = data.frame(subject = rep(c("a", "b"), each = 10), day = rep(1:10, 2), value = rnorm(20))
df$value[15:17] = NA

# lag function (sensitive to "single day" subjects)
lag <- function(x, l = 1) { 
  if (length(x) > 1) (c(rep(NA, l), x[1:(length(x)-l)])) else (NA) 
} 

# calculate rolling mean
df$roll_mean3 = unlist(tapply(df$value, df$subject, 
                              FUN = function(x) lag(rollapplyr(x, width = 3, fill = NA, partial = T,
                                                               FUN = function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)))))
df

Как я уже сказал, это решение кажется чрезмерно сложным для ситуации, которая, я думаю, не так уж далека от реальности.

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

Для иллюстрации вывод моего кода выглядит следующим образом:

    subject day      value   roll_mean3
1        a   1 -0.6264538           NA
2        a   2  0.1836433           NA
3        a   3 -0.8356286 -0.221405243
4        a   4  1.5952808 -0.426146366
5        a   5  0.3295078  0.314431838
6        a   6 -0.8204684  0.363053321
7        a   7  0.4874291  0.368106730
8        a   8  0.7383247 -0.001177187
9        a   9  0.5757814  0.135095124
10       a  10 -0.3053884  0.600511703
11       b   1  1.5117812           NA
12       b   2  0.3898432           NA
13       b   3 -0.6212406  0.950812202
14       b   4 -2.2146999  0.426794608
15       b   5         NA -0.815365744
16       b   6         NA -1.417970234
17       b   7         NA           NA
18       b   8  0.9438362           NA
19       b   9  0.8212212           NA
20       b  10  0.5939013  0.882528703

1. Не могли бы вы пересмотреть свой пост, чтобы включить воспроизводимые выборочные данные (используйте set.seed для фиксированного случайного начального значения) и соответствующий ожидаемый результат? Похоже, что rollapply with dplyr::group_by может выполнить эту работу, но ваш ожидаемый результат поможет лучше понять вашу проблему. Кроме того, dplyr уже имеет lag функцию (и data.table имеет эквивалентную shift функцию).

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

3. Хм. Я не уверен, что понимаю ваш ожидаемый результат. Вы хотели бы рассчитать «среднее значение значений (a) за предыдущие 3 дня, исключая «текущий» день» . Не будет ли первое не- NA значение для subject "a" начинаться в строке 4? Почему у вас есть значение в строке 3?

4. Потому что среднее вычисляется, как только в скользящем окне появляются 2 не пропущенных значения.

5. Но, согласно вашему описанию, вы хотите вычислить среднее значение за value предыдущие 3 дня, исключая текущий день . Как это согласуется с вычислением среднего , «как только в скользящем окне появятся 2 не пропущенных значения» ? Похоже, это две разные цели.

Ответ №1:

Используйте ave для выполнения rollapply отдельно по каждому объекту. Затем при использовании rollapply обратите внимание, что width может быть списком, содержащим вектор (или векторы) смещений, что list(-seq(3)) означает предыдущие 3 элемента. Смотрите ?rollapply для получения дополнительной информации об аргументах.

 Mean <- function(x) if (sum(!is.na(x)) >= 2) mean(x, na.rm = TRUE) else NA
roll <- function(x)  rollapply(x, list(-seq(3)), Mean, fill = NA, partial = TRUE)
transform(df, roll = ave(value, subject, FUN = roll))

1. Это приятно. Мне нравится, что вам не нужно использовать функцию задержки. Я не знал, что rollapply обладает подобной функциональностью сам по себе (для меня это имеет большой смысл, который он должен иметь).

2. Мне также нравится, что вы также не используете dplyr (или data.table). Однако, это, вероятно, просто личный вкус…

Ответ №2:

В дополнение к моему комментарию выше, я не совсем уверен, каким должен быть ожидаемый вами результат, но, возможно, хорошей отправной точкой является следующее:

 df %>%
    group_by(subject) %>%
    mutate(roll_mean3 = rollapplyr(
        lag(value),
        width = 3,
        fill = NA,
        FUN = function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)))
## A tibble: 20 x 4
## Groups:   subject [2]
#   subject   day   value roll_mean3
#   <fct>   <int>   <dbl>      <dbl>
# 1 a           1  -0.626   NA
# 2 a           2   0.184   NA
# 3 a           3  -0.836   -0.221
# 4 a           4   1.60    -0.426
# 5 a           5   0.330    0.314
# 6 a           6  -0.820    0.363
# 7 a           7   0.487    0.368
# 8 a           8   0.738   -0.00118
# 9 a           9   0.576    0.135
#10 a          10  -0.305    0.601
#11 b           1   1.51    NA
#12 b           2   0.390   NA
#13 b           3  -0.621    0.951
#14 b           4  -2.21     0.427
#15 b           5  NA       -0.815
#16 b           6  NA       -1.42
#17 b           7  NA       NA
#18 b           8   0.944   NA
#19 b           9   0.821   NA
#20 b          10   0.594    0.883

Или с помощью data.table

 custom_mean <- function(x) ifelse(sum(!is.na(x)) > 1, mean(x, na.rm = T), NA)
setDT(df)[, roll_mean3 := rollapplyr(shift(value), width = 3, fill = NA, FUN = custom_mean), by = subject]
df
#   subject day      value   roll_mean3
#1:       a   1 -0.6264538           NA
#2:       a   2  0.1836433           NA
#3:       a   3 -0.8356286 -0.221405243
#4:       a   4  1.5952808 -0.426146366
#5:       a   5  0.3295078  0.314431838
#6:       a   6 -0.8204684  0.363053321
#7:       a   7  0.4874291  0.368106730
#8:       a   8  0.7383247 -0.001177187
#9:       a   9  0.5757814  0.135095124
#10:       a  10 -0.3053884  0.600511703
#11:       b   1  1.5117812           NA
#12:       b   2  0.3898432           NA
#13:       b   3 -0.6212406  0.950812202
#14:       b   4 -2.2146999  0.426794608
#15:       b   5         NA -0.815365744
#16:       b   6         NA -1.417970234
#17:       b   7         NA           NA
#18:       b   8  0.9438362           NA
#19:       b   9  0.8212212           NA
#20:       b  10  0.5939013  0.882528703

1. Отличная идея. Возможно, вы также захотите включить library(zoo) в свой код

2. хм, да, это решение. Однако он заменяет unlist(tapply()) только с помощью dplyrs group_by и мою собственную функцию задержки с помощью dplyrs lag function. В остальном это выглядит так же.

3. @MrMax Да, я согласен. Не сильно отличается, кроме (возможно) того, что оно немного более читабельно; Я не был ясен в вашей постановке проблемы и до сих пор в замешательстве. Спасибо за добавление ожидаемого результата, который, по-видимому, предполагает, что приведенный выше dplyr код не является правильным решением.

4. @MrMax Я внес обновление на основе выборочных данных с помощью set.seed(1) .

5. Не беспокойтесь @MrMax; я также добавил data.table решение (которое должно быть быстрее, чем dplyr метод).

Ответ №3:

Возможно, это не самое элегантное или масштабируемое решение, но оно обеспечивает желаемый результат:

 df %>%
  group_by(subject) %>%
  mutate(n_values = 3 - is.na(lag(value, 1)) - is.na(lag(value, 2)) - is.na(lag(value, 3)),
         roll_mean = ifelse(
           n_values >= 2,
           (coalesce(lag(value), 0)   coalesce(lag(value, 2), 0)   coalesce(lag(value, 3), 0)) / n_values,
           NA)
  )

Объяснение: это dplyr конвейер, который сначала группируется по темам, чтобы группы соблюдались. Далее, есть два вычисленных значения в mutate :

n_values подсчитывает количество значений, отличных от NA, в предыдущих 3 строках, оно равно 3 минус 1 для каждого значения NA. Доступ к предыдущим строкам осуществляется с помощью lag .
roll_mean является условным, используя ifelse : если n_values по крайней мере равно 2, среднее значение может быть вычислено. Он суммирует предыдущие 3 значения, заменяя NAs на 0 с помощью coalesce . Сумма делится на n_values , чтобы получить среднее значение. Если n_values < 2 , возвращается NA.

1. Это интересное решение. Однако, как пользователь R среднего уровня, мне трудно понять, как это работает. Поскольку я ищу более простые решения, возможно, это не то направление, в котором я предполагал.

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