R 计算多年来每一天的趋势

Question

我想计算几年来每一天的趋势。例如从 2000 年到 2010 年 5 月 1 日的趋势。这是我的测试数据框：

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 是值。在 testframe$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-values、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)

Answer 1

从您的输出看来，您希望为每个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)

Answer 2

这为每个 V 列的一年中的每一天提供了单独的截距和斜率。 （yday 是一年中的第 0、1、2、... 天，ydayf 是相同的，但作为一个因素，yr 是数字 4 位年份。）

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

Modern Applied Statistics with S Plus 一书是lm 公式的很好参考。