R中的中断时间序列回归；创建函数以进行多次迭代答案

【问题标题】：Interrupted time series regression in R; creating function to conduct multiple iterationsR中的中断时间序列回归；创建函数以进行多次迭代
【发布时间】：2021-04-19 13:38:59
【问题描述】：

我正在进行间断时间序列回归分析，以帮助查看在特定事件之后值是否存在显着的非零变化。这是两个模拟数据框；一个具有日期和值，另一个具有事件名称和该事件的相应日期：-


#dataset 1
eventDate<-structure(c(18262, 18263, 18264, 18265, 18266, 18267, 18268, 
                       18269, 18270, 18271, 18272, 18273, 18274, 18275, 18276, 18277, 
                       18278, 18279, 18280, 18281, 18282, 18283, 18284, 18285, 18286, 
                       18287, 18288, 18289, 18290, 18291, 18292, 18293, 18294, 18295, 
                       18296, 18297, 18298, 18299, 18300, 18301, 18302, 18303, 18304, 
                       18305, 18306, 18307, 18308, 18309, 18310, 18311, 18312, 18313, 
                       18314, 18315, 18316, 18317, 18318, 18319, 18320, 18321, 18322, 
                       18323, 18324, 18325, 18326, 18327, 18328, 18329, 18330, 18331, 
                       18332, 18333, 18334, 18335, 18336, 18337, 18338, 18339, 18340, 
                       18341, 18342, 18343, 18344, 18345, 18346, 18347, 18348, 18349, 
                       18350, 18351, 18352, 18353, 18354, 18355, 18356, 18357, 18358, 
                       18359, 18360, 18361, 18362, 18363, 18364, 18365, 18366, 18367, 
                       18368, 18369, 18370, 18371, 18372, 18373, 18374, 18375, 18376, 
                       18377, 18378, 18379, 18380, 18381, 18382, 18383, 18384, 18385, 
                       18386, 18387, 18388, 18389, 18390, 18391, 18392, 18393, 18394, 
                       18395, 18396, 18397, 18398, 18399, 18400, 18401, 18402, 18403, 
                       18404, 18405, 18406, 18407, 18408, 18409, 18410, 18411, 18412, 
                       18413, 18414, 18415, 18416, 18417, 18418, 18419, 18420, 18421, 
                       18422, 18423, 18424, 18425, 18426, 18427, 18428, 18429, 18430, 
                       18431, 18432, 18433, 18434, 18435, 18436, 18437, 18438, 18439, 
                       18440, 18441, 18442, 18443, 18444, 18445, 18446, 18447, 18448, 
                       18449, 18450, 18451, 18452, 18453, 18454, 18455, 18456, 18457, 
                       18458, 18459, 18460, 18461, 18462, 18463, 18464, 18465, 18466, 
                       18467, 18468, 18469, 18470, 18471, 18472, 18473, 18474, 18475, 
                       18476, 18477, 18478, 18479, 18480, 18481, 18482, 18483, 18484, 
                       18485, 18486, 18487, 18488, 18489, 18490, 18491, 18492, 18493, 
                       18494, 18495, 18496, 18497, 18498, 18499, 18500, 18501, 18502, 
                       18503, 18504, 18505, 18506, 18507, 18508, 18509, 18510, 18511, 
                       18512, 18513, 18514, 18515, 18516, 18517, 18518, 18519, 18520, 
                       18521, 18522, 18523, 18524, 18525, 18526, 18527, 18528, 18529, 
                       18530, 18531, 18532, 18533, 18534, 18535, 18536, 18537, 18538, 
                       18539, 18540, 18541, 18542, 18543, 18544, 18545, 18546, 18547, 
                       18548, 18549, 18550, 18551, 18552, 18553, 18554, 18555, 18556, 
                       18557, 18558, 18559, 18560, 18561, 18562, 18563, 18564, 18565, 
                       18566, 18567, 18568, 18569, 18570, 18571, 18572, 18573, 18574, 
                       18575, 18576, 18577, 18578, 18579, 18580, 18581, 18582, 18583, 
                       18584, 18585, 18586, 18587, 18588, 18589, 18590, 18591, 18592, 
                       18593, 18594, 18595, 18596, 18597, 18598, 18599, 18600, 18601, 
                       18602, 18603, 18604, 18605, 18606, 18607, 18608, 18609, 18610, 
                       18611, 18612, 18613, 18614, 18615, 18616, 18617, 18618, 18619, 
                       18620, 18621, 18622, 18623, 18624, 18625, 18626, 18627), class = "Date")

Count<-c(46L, 58L, 46L, 60L, 42L, 56L, 44L, 60L, 48L, 43L, 50L, 45L, 
         55L, 57L, 47L, 46L, 42L, 44L, 43L, 58L, 60L, 58L, 50L, 55L, 51L, 
         43L, 60L, 51L, 51L, 59L, 44L, 44L, 42L, 48L, 50L, 60L, 53L, 57L, 
         56L, 60L, 52L, 43L, 50L, 55L, 49L, 53L, 50L, 48L, 45L, 51L, 59L, 
         56L, 53L, 45L, 52L, 122L, 100L, 91L, 82L, 60L, 55L, 58L, 42L, 
         53L, 59L, 47L, 58L, 54L, 56L, 49L, 46L, 41L, 48L, 40L, 59L, 43L, 
         46L, 59L, 47L, 51L, 54L, 46L, 46L, 53L, 50L, 51L, 57L, 48L, 60L, 
         59L, 46L, 53L, 50L, 44L, 42L, 58L, 55L, 59L, 57L, 42L, 52L, 43L, 
         54L, 47L, 53L, 44L, 48L, 42L, 56L, 59L, 46L, 49L, 47L, 52L, 58L, 
         42L, 52L, 41L, 55L, 56L, 58L, 52L, 43L, 40L, 56L, 47L, 46L, 50L, 
         45L, 54L, 53L, 50L, 53L, 48L, 58L, 40L, 43L, 55L, 41L, 46L, 46L, 
         55L, 46L, 46L, 52L, 59L, 59L, 46L, 43L, 59L, 57L, 57L, 41L, 40L, 
         44L, 47L, 55L, 44L, 54L, 58L, 56L, 43L, 58L, 45L, 53L, 42L, 57L, 
         59L, 42L, 40L, 53L, 60L, 58L, 40L, 59L, 54L, 41L, 59L, 48L, 48L, 
         43L, 47L, 50L, 50L, 53L, 50L, 43L, 41L, 43L, 51L, 53L, 40L, 52L, 
         43L, 53L, 51L, 51L, 49L, 53L, 40L, 45L, 59L, 50L, 60L, 60L, 42L, 
         47L, 47L, 45L, 50L, 46L, 60L, 40L, 48L, 43L, 59L, 58L, 55L, 48L, 
         44L, 53L, 60L, 52L, 54L, 42L, 44L, 52L, 51L, 47L, 53L, 45L, 41L, 
         56L, 45L, 56L, 52L, 57L, 48L, 47L, 52L, 58L, 51L, 41L, 53L, 155L, 
         123L, 98L, 90L, 84L, 71L, 58L, 50L, 54L, 45L, 58L, 48L, 49L, 
         60L, 41L, 60L, 46L, 40L, 50L, 49L, 57L, 58L, 56L, 58L, 51L, 53L, 
         41L, 45L, 58L, 51L, 50L, 56L, 42L, 59L, 42L, 53L, 57L, 52L, 50L, 
         43L, 43L, 59L, 41L, 54L, 56L, 45L, 42L, 55L, 50L, 58L, 48L, 54L, 
         41L, 42L, 44L, 46L, 50L, 40L, 51L, 42L, 55L, 44L, 54L, 51L, 45L, 
         58L, 40L, 46L, 46L, 40L, 57L, 53L, 40L, 49L, 52L, 50L, 50L, 59L, 
         42L, 57L, 55L, 52L, 57L, 52L, 43L, 48L, 48L, 40L, 42L, 48L, 41L, 
         47L, 53L, 40L, 50L, 60L, 43L, 60L, 51L, 56L, 53L, 42L, 52L, 56L, 
         48L, 44L, 44L, 53L, 60L, 55L, 56L, 45L, 41L, 54L, 43L, 59L, 48L, 
         46L, 50L, 47L, 56L, 57L)


data<-data.frame(eventDate,Count)


#dataset 2
Event<-c("event_a", "event_b", "event_c", "event_d", "event_e", "event_f", 
  "event_g")

Event_Date<-structure(c(18289, 18317, 18358, 18444, 18506, 18528, 18547), class = "Date")

events_dates<-data.frame(Event, Event_Date)

以下是我正在做的两个示例。一个例子有显着的结果，一个例子有一个不显着的结果。

library(dplyr)
library(ggplot2)



#significant result
event_date<-as.Date("2020-09-01")
before_period<-as.Date(event_date)-21
after_period<-as.Date(event_date)+21

data_filtered<-data%>%
  filter(eventDate>=as.Date(before_period) & eventDate<=as.Date(after_period))

data_filtered<-data_filtered%>%
  mutate(DayNumber=row_number())
data_filtered$intv_trend <- cumsum(data_filtered$eventDate >= as.Date(event_date))
data_filtered$Post_event<-ifelse(data_filtered$eventDate<event_date,0,1)
data_filtered <- data_filtered %>%
  mutate(lag_count = lag(Count))

data_filtered

fit <- glm(Count ~ DayNumber+ Post_event + intv_trend+log(lag_count),
                       family = "poisson", data = data_filtered)


summary(fit)

data_filtered$fit <- exp(c(NA, predict(fit)))
data_filtered$fit2 = c(NA, predict(fit, type="response"))

data_filtered$Group<-ifelse(data_filtered$Post_event==0, "Pre-event","Post-event")
data_filtered$Group<-factor(data_filtered$Group, levels = c("Pre-event","Post-event"))


ggplot(data_filtered, aes(x=DayNumber, y = Count, colour=Group)) + 
  geom_line()+geom_smooth(method="lm", se=F, aes(colour=Group)) +ggtitle("Count before and after event")+
            geom_vline(xintercept = 22,linetype="dotted")+labs(caption = "Dotted line represents time of event")

ggplot(data_filtered, aes(x = DayNumber, y = fit2)) +
  geom_line() +
  geom_smooth(method="lm", se=F, aes(colour=Group)) +
  theme_bw() +
  labs(colour="")+ggtitle("Count (Fitted); Method=lm")+
  geom_vline(xintercept = 22,linetype="dotted")+labs(caption = "Dotted line represents time of event")



#non-significant result
event_date<-as.Date("2020-01-28")
before_period<-as.Date(event_date)-21
after_period<-as.Date(event_date)+21

data_filtered<-data%>%
  filter(eventDate>=as.Date(before_period) & eventDate<=as.Date(after_period))

data_filtered<-data_filtered%>%
  mutate(DayNumber=row_number())
data_filtered$intv_trend <- cumsum(data_filtered$eventDate >= as.Date(event_date))
data_filtered$Post_event<-ifelse(data_filtered$eventDate<event_date,0,1)
data_filtered <- data_filtered %>%
  mutate(lag_count = lag(Count))

data_filtered

fit <- glm(Count ~ DayNumber+ Post_event + intv_trend+log(lag_count),
           family = "poisson", data = data_filtered)


summary(fit)

data_filtered$fit <- exp(c(NA, predict(fit)))
data_filtered$fit2 = c(NA, predict(fit, type="response"))

data_filtered$Group<-ifelse(data_filtered$Post_event==0, "Pre-event","Post-event")
data_filtered$Group<-factor(data_filtered$Group, levels = c("Pre-event","Post-event"))

ggplot(data_filtered, aes(x=DayNumber, y = Count, colour=Group)) + 
  geom_line()+geom_smooth(method="lm", se=F, aes(colour=Group)) +ggtitle("Count before and after event")+
  geom_vline(xintercept = 22,linetype="dotted")+labs(caption = "Dotted line represents time of event")

ggplot(data_filtered, aes(x = DayNumber, y = fit2)) +
  geom_line() +
  geom_smooth(method="lm", se=F, aes(colour=Group)) +
  theme_bw() +
  labs(colour="")+ggtitle("Count (Fitted); Method=lm")+
  geom_vline(xintercept = 22,linetype="dotted")+labs(caption = "Dotted line represents time of event")

每个event_date 变量都来自event_dates 数据框。

我有一长串想要进行分析的日期，所以我需要一个函数来帮助我有效地完成这项工作。我希望它返回回归拟合摘要和每个事件的两个图。到目前为止，这是我的尝试（如您所知，不是一个好的尝试）：-

#data1 represents data, data2 represents events_dates, n_days is number of days before and after event_date that I want to filter from data 
TS_Intervention_Func<-function(data1,data2,n_days){
  
  myresultslist<-list()
  
  
  for (i in data2$Event_Date) {
    
    event_date<-as.Date(i)
    before_period<-as.Date(event_date)-n_days
    after_period<-as.Date(event_date)+n_days
    
    
    data_filtered<-data1%>%
      filter(eventDate>=as.Date(before_period) & eventDate<=as.Date(after_period))
    
    
    data_filtered<-data_filtered%>%
      mutate(DayNumber=row_number())
    data_filtered$intv_trend <- cumsum(data_filtered$eventDate >= as.Date(event_date))
    data_filtered$Post_event<-ifelse(data_filtered$eventDate<event_date,0,1)
    data_filtered <- data_filtered %>%
      mutate(lag_count = lag(Count))
    
    
    fit <- glm(Count ~ DayNumber+ Post_event + intv_trend+log(lag_count),
               family = "poisson", data = data_filtered)
    
    
    results<-summary(fit)
    
    data_filtered$fit <- exp(c(NA, predict(fit)))
    data_filtered$fit2 = c(NA, predict(fit, type="response"))
    
    data_filtered$Group<-ifelse(data_filtered$Post_event==0, "Pre-event","Post-event")
    data_filtered$Group<-factor(data_filtered$Group, levels = c("Pre-event","Post-event"))
    
    
    plot_1<-ggplot(data_filtered, aes(x=DayNumber, y = Count, colour=Group)) + 
      geom_line()+geom_smooth(method="lm", se=F, aes(colour=Group)) +ggtitle("Count before and after event",paste0(event_date))+
      geom_vline(xintercept = 22,linetype="dotted")+labs(caption = "Dotted line represents time of event")
    
    plot_2<-ggplot(data_filtered, aes(x = DayNumber, y = fit2)) +
      geom_line() +
      geom_smooth(method="lm", se=F, aes(colour=Group)) +
      theme_bw() +
      labs(colour="")+ggtitle("Count (Fitted); Method=lm", paste0(event_date))+
      geom_vline(xintercept = 22,linetype="dotted")+labs(caption = "Dotted line represents time of event")
    
    
    myresultslist[[i]] <- do.call(results,plot1,plot_2) 
    
  }
return(myresultslist)
}



TS_Intervention_Func(data,events_dates,21)
#Error in as.Date.numeric(i) : 'origin' must be supplied

简而言之，我希望函数做三件事：-

从event_dates 数据框中获取每个日期，并在data 中对该日期运行分析迭代
将fit摘要和两个对应的图表存储到该事件日期并将其保存在列表中（如果可以将事件名称与它一起保存以便于查找，那就更好了）
这是可取的；如果列表可以分为两部分，一个保存重要结果，另一个保存非重要结果。

一个很大的问题，但我们一如既往地感谢任何帮助:)

【问题讨论】：

标签： r list function for-loop regression

【解决方案1】：

这不一定是最有效的方法，但我想让它简单易懂。

我稍微清理了您的代码并将其放入一个函数中。首先，我编写了一个处理一个事件日期的函数：

library(dplyr)
library(ggplot2)


reg_fun <- function(event_date, data) {


  before_period <- as.Date(event_date) - 21
  after_period <- as.Date(event_date) + 21
  
  data_filtered <- data %>%
    filter(eventDate >= as.Date(before_period) &
             eventDate <= as.Date(after_period)) %>% 
    mutate(DayNumber = row_number()) %>% 
    mutate(intv_trend = cumsum(eventDate >= event_date)) %>% 
    mutate(Post_event = as.integer(!eventDate < event_date)) %>% 
    mutate(lag_count = lag(Count))
  
  
  fit <- glm(Count ~ DayNumber+ Post_event + intv_trend+log(lag_count),
             family = "poisson", data = data_filtered)
  
  
  fit_summary <- summary(fit)
  
  
  
  data_filtered$fit <- exp(c(NA, predict(fit)))
  data_filtered$fit2 = c(NA, predict(fit, type = "response"))
  
  data_filtered$Group <-
    ifelse(data_filtered$Post_event == 0, "Pre-event", "Post-event")
  data_filtered$Group <-
    factor(data_filtered$Group, levels = c("Pre-event", "Post-event"))
  
  
  plot1 <- ggplot(data_filtered, aes(x = DayNumber, y = Count, colour = Group)) +
    geom_line() + geom_smooth(method = "lm", se = F, aes(colour = Group)) +
    ggtitle("Count before and after event") +
    geom_vline(xintercept = 22, linetype = "dotted") + labs(caption = "Dotted line represents time of event")
  
  plot2 <- ggplot(data_filtered, aes(x = DayNumber, y = fit2)) +
    geom_line() +
    geom_smooth(method = "lm", se = F, aes(colour = Group)) +
    theme_bw() +
    labs(colour = "") + ggtitle("Count (Fitted); Method=lm") +
    geom_vline(xintercept = 22, linetype = "dotted") + labs(caption = "Dotted line represents time of event")
  
  output <- list(
    event_date = event_date,
    fit = fit,
    summary = fit_summary,
    plot1 = plot1,
    plot2 = plot2,
    signif = sum(coef(fit_summary)[, 4] < 0.05) > 1 # see if any is significant (p-value for intercept is always < 0.05 so we want more than one significant value)
  )
  
  return(output)
}

现在我们可以用一个日期来测试它。如上所示，输出是一个列表，我认为在这种情况下最合适。

test_res <- reg_fun(event_date = as.Date("2020-09-01"), data = data)

# print summary
test_res$summary
#> 
#> Call:
#> glm(formula = Count ~ DayNumber + Post_event + intv_trend + log(lag_count), 
#>     family = "poisson", data = data_filtered)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -2.9155  -0.8267  -0.1498   0.8527   4.9107  
#> 
#> Coefficients:
#>                 Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)     3.649360   0.353805  10.315  < 2e-16 ***
#> DayNumber       0.004957   0.005517   0.899    0.369    
#> Post_event      0.714854   0.098185   7.281 3.32e-13 ***
#> intv_trend     -0.051807   0.008034  -6.448 1.13e-10 ***
#> log(lag_count)  0.050323   0.089530   0.562    0.574    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 278.729  on 41  degrees of freedom
#> Residual deviance:  95.044  on 37  degrees of freedom
#>   (1 observation deleted due to missingness)
#> AIC: 350.86
#> 
#> Number of Fisher Scoring iterations: 4

# plot 1
test_res$plot1
#> `geom_smooth()` using formula 'y ~ x'


# plot 2
test_res$plot2
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 row(s) containing missing values (geom_path).

我不显示图，但你明白了。

现在我将它包装在一个可以同时获取多个日期的函数中。理论上我们可以在同一个函数中做到这一点。

reg_fun_mult <- function(event_dates, dat) {
  output <- lapply(event_dates, reg_fun, dat)
  names(output) <- event_dates # give the list elements suitable names
  return(output)
}

test_res_mult <- reg_fun_mult(eventDate, data)

# check if any variables were significant
sign <- sapply(test_res_mult, function(x) x$signif)
# look at the first ten results to see which ones have significant values
sign[1:10]
#> 2020-01-01 2020-01-02 2020-01-03 2020-01-04 2020-01-05 2020-01-06 2020-01-07 
#>      FALSE      FALSE      FALSE      FALSE      FALSE      FALSE      FALSE 
#> 2020-01-08 2020-01-09 2020-01-10 
#>      FALSE      FALSE      FALSE


# keep only significant results
test_res_mult_sign <- test_res_mult[sign]

您可以再次查看单个摘要和图表，如下所示：

# summary
test_res_mult$`2020-01-01`$summary
#> 
#> Call:
#> glm(formula = Count ~ DayNumber + Post_event + intv_trend + log(lag_count), 
#>     family = "poisson", data = data_filtered)
#> 
#> Deviance Residuals: 
#>     Min       1Q   Median       3Q      Max  
#> -1.3085  -0.8371  -0.2015   0.9650   1.4188  
#> 
#> Coefficients: (2 not defined because of singularities)
#>                  Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)     4.3952688  0.9248944   4.752 2.01e-06 ***
#> DayNumber       0.0001341  0.0050776   0.026    0.979    
#> Post_event             NA         NA      NA       NA    
#> intv_trend             NA         NA      NA       NA    
#> log(lag_count) -0.1213342  0.2358350  -0.514    0.607    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 18.972  on 20  degrees of freedom
#> Residual deviance: 18.705  on 18  degrees of freedom
#>   (1 observation deleted due to missingness)
#> AIC: 145.57
#> 
#> Number of Fisher Scoring iterations: 4

# plot 1
test_res_mult$`2020-01-01`$plot1
#> `geom_smooth()` using formula 'y ~ x'


# plot 2
test_res_mult$`2020-01-01`$plot2
#> `geom_smooth()` using formula 'y ~ x'
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).

#> Warning: Removed 1 row(s) containing missing values (geom_path).

如果有什么不清楚的地方请告诉我。

【讨论】：

这是一个很好的解决方案，我认为它非常有效，并且可以为我省去很多麻烦！我会在允许的情况下尽快奖励赏金。我想知道，有没有稍微调整一下，所以它只包括重要的Post_event 结果？这将是我最感兴趣的，但如果它太难修改，请不要强调。
是的，您可以运行sign <- sapply(test_res_mult, function(x) x$signif)，然后运行test_res_mult[sign]，只获得具有显着结果的迭代。因为我想打印sign 来显示正在发生的事情，所以我可能已经把那部分隐藏了一点。
太好了，感谢您抽出宝贵时间编写此解决方案！
不客气 :) 我最近为我自己的作品写了一些非常相似的东西，有点难过它是如此具体以至于没有分享的意义。现在我至少可以在这里分享它:)
至少你有机会分享它！当你为一个小众问题想出一些东西，然后以某种能力再次使用它时，这很好。