#r #time-series #data-modeling #smoothing
#r #временные ряды #моделирование данных #сглаживание
Вопрос:
У меня есть данные о процентном показателе по дате и времени (дата и часы: минуты: секунды). Я хочу графически «исправить» / выделить точку данных, которая не является репрезентативной.
Предыстория
У меня есть данные о том, как люди ежедневно оценивают свой уровень счастья по непрерывной шкале от 0 до> 1, где 0 означает «крайне несчастный», а 1 означает «чрезвычайно счастливый». Я спрашиваю многих людей и хочу со временем почувствовать «счастье в группе».
Данные
library(tidyverse)
library(lubridate)
set.seed(1234)
original_df <-
seq(as.POSIXct('2020-09-01', tz = "UTC"), as.POSIXct('2020-09-15', tz = "UTC"), by="1 mins") %>%
sample(15000, replace = T) %>%
as_tibble %>%
rename(date_time = value) %>%
mutate(date = date(date_time)) %>%
add_column(score = runif(15000))
original_df
## # A tibble: 15,000 x 3
## date_time date score
## <dttm> <date> <dbl>
## 1 2020-09-06 04:11:00 2020-09-06 0.683
## 2 2020-09-06 13:35:00 2020-09-06 0.931
## 3 2020-09-05 23:21:00 2020-09-05 0.121
## 4 2020-09-06 14:45:00 2020-09-06 0.144
## 5 2020-09-07 09:15:00 2020-09-07 0.412
## 6 2020-09-01 10:22:00 2020-09-01 0.564
## 7 2020-09-11 14:00:00 2020-09-11 0.960
## 8 2020-09-08 13:24:00 2020-09-08 0.845
## 9 2020-09-01 15:33:00 2020-09-01 0.225
## 10 2020-09-09 19:27:00 2020-09-09 0.815
## # ... with 14,990 more rows
Однако оказывается, что в один из дней оказывается значительно меньше точек данных. Таким образом, фактический набор данных выглядит следующим образом:
actual_df <-
original_df %>%
filter(date %in% as_date("2020-09-10")) %>%
group_by(date) %>%
slice_sample(n = 15) %>%
ungroup %>%
bind_rows(original_df %>% filter(!date %in% as_date("2020-09-10")))
> actual_df %>% count(date)
## # A tibble: 14 x 2
## date n
## <date> <int>
## 1 2020-09-01 1073
## 2 2020-09-02 1079
## 3 2020-09-03 1118
## 4 2020-09-04 1036
## 5 2020-09-05 1025
## 6 2020-09-06 1089
## 7 2020-09-07 1040
## 8 2020-09-08 1186
## 9 2020-09-09 1098
## 10 2020-09-10 15 ## <- this day has less data
## 11 2020-09-11 1095
## 12 2020-09-12 1051
## 13 2020-09-13 1037
## 14 2020-09-14 1034
Построение этих данных с течением времени
То, что я делал, основывалось на работе со средствами
Я рассматриваю каждый день как фактор и получаю среднее значение за день. По статистике, это решение может быть далеко от идеального, как прокомментировал @BrianLang ниже. Однако прямо сейчас я выбрал именно этот метод.
library(emmeans)
model_fit <-
actual_df %>%
mutate(across(date, factor)) %>%
lm(score ~ date, data = .)
emmeans_fit_data <- emmeans(model_fit, ~ date, CIs = TRUE)
emmeans_fit_data %>%
as_tibble %>%
ggplot(data = ., aes(x = date, y = emmean))
geom_line(color = "#1a476f", group = 1, lwd = 1)
geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), alpha = 0.5, color = "#90353b", width = 0.2)
geom_text(aes(label = paste0(round(100*emmean, 1), "%") , color = "90353b"), vjust = -4, hjust = 0.5, size = 3.5)
geom_point(color = "1a476f")
scale_y_continuous(labels = function(x) paste0(100*x, "%"))
ylab("Level of Happiness")
xlab("Date")
ggtitle("Mood Over Time")
theme(plot.title = element_text(hjust = 0.5, size = 14),
axis.text.x=element_text(angle = -60, hjust = 0),
axis.title.x = element_blank(),
legend.title = element_blank(),
plot.caption = element_text(hjust = 0, size = 8),
legend.position = "none")
Но затем я получаю этот всплеск в 2020-09-10, что связано только с малым количеством точек данных.
Одним из графических решений было бы сделать что-то вроде удаления проблемной линии и «завершения» того, как это выглядело бы с достаточным количеством точек данных. Возможно, на основе усреднения за день до и на следующий день? Я не хочу избавляться от реальных данных, но хочу графически подчеркнуть, что это нерепрезентативно, и что реальное значение должно было быть намного ближе к дню до и после. Я думал, что использование пунктирных линий является разумным графическим решением.
В противном случае я надеялся, что может быть другой метод моделирования / построения таких «временных» данных с использованием ggplot сглаживания, что даст мне более плавную линию тренда и доверительную ленту, которая будет учитывать проблемный день. Но я понимаю, что это может выходить за рамки этого вопроса, поэтому я просто добавляю его в качестве примечания; на случай, если кто-то захочет предложить решение, основанное на другом моделировании, вместо простой графической коррекции. Но я буду благодарен и за то, и за другое.
Комментарии:
1. Похоже, что вы подходите для этой модели с данными как категорией. Я думаю, вам, вероятно, следует рассматривать date как непрерывную переменную. Ваша проблема в том, что вы не понимаете имеющиеся у вас данные или модели, которые были бы уместны для данных. Для этого вам необходимо узнать о моделях, потенциально моделях arima.
2. Спасибо. Да, мой вопрос, безусловно, связан с недостатком знаний, которые я надеюсь получить с помощью вопроса, который я опубликовал. Не могли бы вы уточнить или направить меня дальше? Модели ARIMA — это один из типов, о котором я не знал. Что-нибудь еще, что вы могли бы придумать? Я могу свободно использовать Google, но это то, что я делал до публикации здесь, поэтому мне нужно больше внимания. Спасибо!
Ответ №1:
Не желая углубляться в модели временных рядов, вы можете представить себе преобразование вашей временной переменной с помощью ограниченных кубических сплайнов.
Мне нужно было немного изменить ваш код, чтобы я мог избежать установки новейших версий некоторых пакетов ;-).
Обратите внимание, что я изменил некоторые имена переменных, потому date что это имя функции, и его не следует использовать как имя переменной.
library(chron)
## added a numeric version of your date variable.
actual_df <- original_df %>%
filter(datez %in% lubridate::date("2020-09-10")) %>%
sample_n(size = 15) %>%
group_by(datez) %>%
ungroup %>%
bind_rows(original_df %>% filter(!datez %in% lubridate::date("2020-09-10"))) %>%
mutate(num_date = as.numeric(datez))
## How many knots across the dates do you want?
number_of_knots = 15
## This is to make sure that visreg is passed the actual knot locations! RMS::RCS does not store them in the model fits.
knots <- paste0("c(", paste0(attr(rms::rcs(actual_df$num_date, number_of_knots), "parms"), collapse = ", "), ")")
## We can construct the formula early.
formula <- as.formula(paste("score ~ rms::rcs(num_date,", knots,")"))
## fit the model as a gaussian glm and pass it to visreg for it's prediction function. This will give you predicted means and 95% CI for that mean. Then I convert the numeric dates back to real dates.
glm_rcs <- glm(data = actual_df, formula = formula, family = "gaussian") %>% visreg::visreg(plot = F) %>% .$fit %>%
mutate(date_date = chron::as.chron(num_date) %>% as.POSIXct())
## plot it!
ggplot(data = glm_rcs, aes(date_date,
y = visregFit))
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = .5)
geom_line()
РЕДАКТИРОВАТЬ: вы собираете данные по дням, но вы можете добавить дрожание к дате таким образом, чтобы они распределялись в течение дня.
actual_df <- original_df %>%
filter(datez %in% lubridate::date("2020-09-10")) %>%
sample_n(size = 15) %>%
group_by(datez) %>%
ungroup %>%
bind_rows(original_df %>% filter(!datez %in% lubridate::date("2020-09-10"))) %>%
mutate(num_date = as.numeric(datez)) %>%
## Here we add random noise (uniform -.5 to .5) to each numeric date.
mutate(jittered_date = num_date runif(n(), -.5, .5))
## You can lower this number to increase smoothing.
number_of_knots = 15
knots <- paste0("c(", paste0(attr(rms::rcs(actual_df$jittered_date, number_of_knots), "parms"), collapse = ", "), ")")
formula <- as.formula(paste("score ~ rms::rcs(jittered_date,", knots,")"))
glm_rcs <- glm(data = actual_df, formula = formula, family = "gaussian") %>% visreg::visreg(plot = F) %>% .$fit %>%
mutate(date_date = chron::as.chron(jittered_date) %>% as.POSIXct())
ggplot(data = glm_rcs, aes(date_date,
y = visregFit))
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = .5)
geom_line()
Правка 2:
Дрожание точек не так необходимо, если у вас есть вектор даты и времени, а не простой день. В вашем исходном коде для создания используемых вами поддельных данных lubridate::date() , которые берут ваш вектор posix datetime и сокращают до простой даты! Вы можете избежать этого с помощью чего-то вроде этого:
original_df <- tibble(datez = seq(as.POSIXct('2020-09-01', tz = "UTC"), as.POSIXct('2020-09-15', tz = "UTC"), by="1 mins") %>%
sample(15000, replace = T)) %>%
mutate(datez_day = lubridate::date(datez)) %>%
add_column(score = runif(15000))
actual_df <- original_df %>%
filter(datez_day %in% lubridate::date("2020-09-10")) %>%
sample_n(size = 15) %>%
bind_rows(original_df %>% filter(!datez_day %in% lubridate::date("2020-09-10"))) %>%
mutate(num_date = as.numeric(datez))
теперь у вас есть datez_day , что является значением дня, datez что является датой-временем, а num_date что является числовым представлением даты-времени.
оттуда вы можете напрямую моделировать num_date , не добавляя никакого дрожания.
number_of_knots = 20
knots <- paste0("c(", paste0(attr(rms::rcs(actual_df$num_date, number_of_knots), "parms"), collapse = ", "), ")")
formula <- as.formula(paste("score ~ rms::rcs(num_date,", knots,")"))
glm_rcs <- glm(data = actual_df, formula = formula, family = "gaussian") %>%
visreg::visreg(plot = F) %>%
.$fit %>%
as_tibble() %>%
## Translate the num_date back into a datetime object so it is correct in the figures!
mutate(date_date = as.POSIXct.numeric(round(num_date), origin = "1970/01/01"))
ggplot(data = glm_rcs, aes(date_date,
y = visregFit))
geom_ribbon(aes(ymin = visregLwr, ymax = visregUpr), alpha = .5)
geom_line()
Комментарии:
1. Это очень полезно. Однако в моих реальных данных (а не в тех, которые были смоделированы для вопроса) я получаю более сильное падение, чем показано в вашем решении. (смотрите здесь: i.stack.imgur.com/oNQ3B.png ) Хотя это, очевидно, связано с различиями между наборами данных, и не видя реальных данных, не могли бы вы выдвинуть гипотезу, почему это дает такой доминирующий спад?
2. На самом деле возникает вопрос, который я поднял в сообщении. Если у меня такое доминирующее падение, которое явно связано с малым количеством точек данных за заданное время, то как я могу решить эту проблему, если сглаживание не сглаживает ее? Мотивация для решения этой проблемы исходит из необходимости донести сюжет до других. Я хочу, чтобы люди, впервые увидевшие сюжет, не были отвлечены или введены в заблуждение таким падением. Итак, я либо исправляю это косметически / графически, либо с помощью методологии сглаживания.
3. Вы можете изменить количество узлов, если хотите. Уменьшение этого числа приведет к оценке взаимосвязи с меньшим количеством степеней свободы, увеличивая сглаживание.
4. Я подумал об этом немного больше. У вас есть весь вес данных конкретно для значений x в виде дней. Если вы добавите дрожание к датам таким образом, чтобы они охватывали весь день, мы можем получить более плавные подгонки.
5. Я понимаю. Но если мы учтем, что данные имели время как непрерывную дату и время и не были жестко запрограммированы по дням? Взвешивание по дням было выбором, но можно ли было сделать иначе? Можно ли моделировать данные с помощью X как непрерывный временной поток, а затем отмечать сегменты графика как дни?