在 R 中应用逐年分段回归

Question

我有每日降雨量数据，我已使用以下代码将其转换为年度累积值

library(seas)
library(data.table)
library(ggplot2)

#Loading data
data(mscdata)
dat <- (mksub(mscdata, id=1108447))
dat$julian.date <- as.numeric(format(dat$date, "%j"))
DT <- data.table(dat)
DT[, Cum.Sum := cumsum(rain), by=list(year)]

df <- cbind.data.frame(day=dat$julian.date,cumulative=DT$Cum.Sum)

然后我想逐年应用分段回归以获得逐年断点。我可以像一年一样做到这一点

library("segmented")
x <- subset(dat,year=="1984")$julian.date
y <- subset(DT,year=="1984")$Cum.Sum
fit.lm<-lm(y~x)
segmented(fit.lm, seg.Z = ~ x, npsi=3)

我使用npsi = 3 有3 个断点。现在如何最小化地应用它逐年分段回归并获得估计的断点？

Answer 1

这是一个带有自定义函数的简短脚本，以便您可以运行不同的年度回归。

## using tidyverse processes instead of mixing and matching with other data manipulation packages 
library(tidyverse); library(segmented); library(seas)

## get mscdata from "seas" packages
data(mscdata)
dat <- (mksub(mscdata, id=1108447))

## generate cumulative sum of rain by year
d2 <- dat %>% group_by(year) %>% mutate(rain_cs = cumsum(rain)) %>% ungroup

## write a custom function

segmentedlm <- function(data, year){
  subset.df <- data %>% filter(year == year)
  fit.lm <- lm(rain_cs ~ julian.date, subset.df)
  segmented(fit.lm, seg.Z = ~ julian.date, npsi=3)
}

# run the customised function for 1975 data
segmentedlm(d2, "1975") %>% plot(., main="1975")

segmentedlm(d2, "1984") %>% plot(., main = "1984")

将多年的分段线性模型的摘要输出到文本文件中：

sink("output.txt")
lapply(c("1975", "1984"), function(x) segmentedlm(d2, x))
sink()

您可以将 lapply 的参数更改为输入所有年份。

Answer 2

您可以将lm 对象存储在一个列表中，并为每个year 应用segmented。

library(tidyverse)

data <- DT %>%
         group_by(year) %>%
         summarise(fit.lm = list(lm(Cum.Sum~julian.date)), 
                   julian.date1 = list(julian.date)) %>%
         mutate(out = map2(fit.lm, julian.date1, function(x, julian.date) 
                       data.frame(segmented::segmented(x, 
                                  seg.Z = ~julian.date, npsi=3)$psi))) %>%
         unnest_wider(out) %>%
         unnest(cols = c(Initial, Est., St.Err)) %>%
         dplyr::select(-fit.lm, -julian.date1)

# A tibble: 90 x 4
#    year Initial  Est. St.Err
#   <int>   <dbl> <dbl>  <dbl>
# 1  1975    84.8  68.3  1.44 
# 2  1975   168.  167.   9.31 
# 3  1975   282.  281.   0.917
# 4  1976    84.8  68.3  1.44 
# 5  1976   168.  167.   9.33 
# 6  1976   282.  281.   0.913
# 7  1977    84.8  68.3  1.44 
# 8  1977   168.  167.   9.32 
# 9  1977   282.  281.   0.913
#10  1978    84.8  68.3  1.44 
# … with 80 more rows