R Расчет тренда за каждый день в течение нескольких лет - PullRequest
Новая
       35

R Расчет тренда за каждый день в течение нескольких лет

2 голосов
/ 13 июня 2019

Я хочу рассчитать тренд каждого дня за несколько лет.Например, тренд 1 мая с 2000 по 2010 год. Вот мой тест-фрейм данных: 

library(lubridate)
date_list = seq(ymd('2000-01-15'),ymd('2010-09-18'),by='day')
testframe = data.frame(Date = date_list)
testframe$Day = substr(testframe$Date, start = 6, stop = 10)
testframe$V1 = rnorm(3900)
testframe$V2 = rnorm(3900)
testframe$V3 = seq(from = 10, to = 25, length.out = 3900)
testframe$V4 = seq(from = 5, to = 45, length.out = 3900)

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

В конце я хотел бы иметь фрейм данных, в котором есть столбец, который содержит каждый день (без года курса), и столбцы, которые содержат тренд / наклон значений от V1 до V4.

Есть идеи?

ОБНОВЛЕНИЕ:

Чтобы сделать это более наглядно.Я хочу и вывод, который выглядит следующим образом (тренды являются случайными) 

Day       V1 Trend   V2 Trend    V3 Trend   V4 Trend
01-01     +0.3          +0.4      +0.9        +0.5
01-02     +0.5          +0.3      +0.8        +0.4
01-03     -0.1          -0.2      +1.0        -0.3
01-04     +0.7          -0.7      +0.9        +0.9
......
......
12-30    -0.3           -0.4      +0.5        +0.8
12-31    -0.7           -0.3      +0.6        +0.9

p-значения, Intercept и все было бы также здорово иметь.

Я нашел этот пример, но он все ещене в выводе, который я хочу иметь: 

#Add year for lm    
testframe$Year = as.numeric(format(testframe$Date,'%Y'))
library(plyr)
# Break up d by state, then fit the specified model to each piece and
# return a list
models <- dlply(testframe, "Day", function(df) 
  lm(Year ~ V4, data = df))

# Apply coef to each model and return a data frame
ldply(models, coef)

# Print the summary of each model
l_ply(models, summary, .print = TRUE)

Ответы [ 2 ]

1 голос
/ 13 июня 2019

По вашему выводу это выглядит так: для каждого Day вы хотите построить линейную модель вида V ~ Year для каждого V1, V2, V3, V4.

Вот подход dplyr: 

library(lubridate)
library(dplyr)

set.seed(23)  # for reproducibility

# data (using your code)
date_list = seq(ymd('2000-01-15'),ymd('2010-09-18'),by='day')
testframe = data.frame(Date = date_list)
testframe$Day = substr(testframe$Date, start = 6, stop = 10)
testframe$V1 = rnorm(3900)
testframe$V2 = rnorm(3900)
testframe$V3 = seq(from = 10, to = 25, length.out = 3900)
testframe$V4 = seq(from = 5, to = 45, length.out = 3900)


testframe %>% 
  mutate(Year = year(Date)) %>%  # extract the year
  select(-Date) %>%              # remove the Date column
  group_by(Day) %>%              # for each day
  summarise_at(vars(matches("V")), ~lm(. ~ Year)$coefficients[2])  # build a model and keep the slope

# # A tibble: 366 x 5
#   Day        V1      V2    V3    V4
#   <chr>   <dbl>   <dbl> <dbl> <dbl>
# 1 01-01  0.108   0.0554  1.41  3.75
# 2 01-02 -0.0543 -0.103   1.41  3.75
# 3 01-03 -0.143  -0.0176  1.41  3.75
# 4 01-04  0.146  -0.0232  1.41  3.75
# 5 01-05 -0.154  -0.0533  1.41  3.75
# 6 01-06 -0.268   0.0470  1.41  3.75
# 7 01-07 -0.164   0.0873  1.41  3.75
# 8 01-08  0.0634  0.266   1.41  3.75
# 9 01-09  0.0115 -0.0320  1.41  3.75
# 10 01-10  0.0576 -0.237   1.41  3.75
# # ... with 356 more rows

Если вы хотите обновить имена столбцов на что-то вроде v_trend, вы можете использовать это вместо: 

summarise_at(vars(matches("V")), list(trend = ~lm(. ~ Year)$coefficients[2]))

Альтернатива (получение дополнительной информации о каждой модели)

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

library(lubridate)
library(tidyverse)
library(broom)

testframe %>% 
  mutate(Year = year(Date)) %>%  
  select(-Date) %>% 
  gather(v, value, -Day, -Year) %>%
  group_by(Day, v) %>%
  nest() %>%
  mutate(dd = map(data, ~tidy(lm(value ~ Year, data = .)))) %>%
  unnest(dd) %>%
  arrange(Day)

# # A tibble: 2,928 x 7
#     Day   v     term          estimate  std.error  statistic  p.value
#   <chr> <chr> <chr>            <dbl>      <dbl>      <dbl>    <dbl>
# 1 01-01 V1    (Intercept)  -217.     162.           -1.34  2.16e- 1
# 2 01-01 V1    Year            0.108    0.0806        1.34  2.16e- 1
# 3 01-01 V2    (Intercept)  -112.     196.           -0.570 5.84e- 1
# 4 01-01 V2    Year            0.0554   0.0976        0.567 5.86e- 1
# 5 01-01 V3    (Intercept) -2800.       0.260    -10756.    6.25e-30
# 6 01-01 V3    Year            1.41     0.000130  10824.    5.94e-30
# 7 01-01 V4    (Intercept) -7489.       0.694    -10787.    6.11e-30
# 8 01-01 V4    Year            3.75     0.000346  10824.    5.94e-30
# 9 01-02 V1    (Intercept)   109.     238.            0.458 6.59e- 1
# 10 01-02 V1    Year           -0.0543   0.119        -0.458 6.59e- 1
# # ... with 2,918 more rows

Затем вы можете запросить этот набор данных и получить все, что захотите. Например, если вы сохраните вышеприведенный вывод как testframe2, вы можете получить тренд / наклон для дня 01-01, столбец V1 следующим образом: 

testframe2 %>% filter(Day == "01-01" & v == "V1" & term == "Year") %>% pull(estimate)

и значение p этого наклона: 

testframe2 %>% filter(Day == "01-01" & v == "V1" & term == "Year") %>% pull(p.value)
1 голос
/ 13 июня 2019

Это обеспечивает отдельный перехват и наклон для каждого дня года для каждого столбца V.(yday - это день года 0, 1, 2, ... и ydayf - это то же самое, что и коэффициент, а yr - это год из четырех цифр.) 

m <- as.matrix(testframe[-(1:2)])
yday <- as.POSIXlt(testframe$Date)$yday
ydayf <- factor(yday)
yr <- as.numeric(format(testframe$Date, "%Y"))

fm2 <- lm(m ~ ydayf + ydayf:yr + 0)
coef(fm2)
dummy.coef(fm2) # expand coefficients
summary(fm2)
broom::tidy(fm2) # data frame

Если вам нужны отдельные уклоны, но только один перехват, используйте для каждого столбца V. 

fm3 <- lm(m ~ ydayf:yr)
coef(fm3) 
dummy.coef(fm3) # expands coefficients
summary(fm3)
broom::tidy(fm3)  # data frame

Если вам нужны отдельные уклоны, но только один уклон для столбца V, тогда: 

fm4 <- lm(m ~ ydayf + yr + 0)
coef(fm4) 
dummy.coef(fm4) # expands coefficients
summary(fm4)
broom::tidy(fm4)  # data frame

Книга Современная прикладная статистика с S Plus является хорошим справочником для lm формул.

...