应用多年的分段线性模型

Question

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

library(tidyverse); library(segmented); library(seas); library(SiZer)

## get mscdata from "seas" packages
data(mscdata)
dat <- (mksub(mscdata, id=1108447))
dat$julian.date <- as.numeric(format(dat$date, "%j"))

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

然后我想将每年分成两部分（210 天前和 210 天后），然后应用 SiZer 的分段线性模型来识别年度断点。我可以像

那样做一年

data <- subset(df, year == 1975)

sub1 <- filter(data, julian.date < 210)
sub2 <- filter(data, julian.date > 210)

sub1.mod <- piecewise.linear(x= sub1$julian.date, y = sub1$rain_cs,
                             middle = 1,
                             CI = T,
                             bootstrap.samples = 1000)
sub1.mod

sub2.mod <- piecewise.linear(x= sub2$julian.date, y = sub2$rain_cs,
                             CI = T,
                             bootstrap.samples = 1000)
sub2.mod

现在如何动态拟合所有年份的分段线性模型？

Answer 1

您可以尝试使用一个函数 base R，创建一个列表，然后保存模型。我在最后一行包含了一种导出列表外所有模型的方法：

library(tidyverse); library(segmented); library(seas); library(SiZer)

## get mscdata from "seas" packages
data(mscdata)
dat <- (mksub(mscdata, id=1108447))
dat$julian.date <- as.numeric(format(dat$date, "%j"))

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

#Create list

Listyear <- split(df,df$year)

#Function for year process
model_function<-function(x)
{
  data <- x
  
  sub1 <- filter(data, julian.date < 210)
  sub2 <- filter(data, julian.date > 210)
  
  sub1.mod <- piecewise.linear(x= sub1$julian.date, y = sub1$rain_cs,
                               middle = 1,
                               CI = T,
                               bootstrap.samples = 1000)
  sub1.mod
  
  sub2.mod <- piecewise.linear(x= sub2$julian.date, y = sub2$rain_cs,
                               CI = T,
                               bootstrap.samples = 1000)
  sub2.mod
  
  #Group elements
  list.model <- list(v1=sub1.mod,v2=sub2.mod)
  names(list.model)<-paste0(c("sub.mod1.","sub.mod2."),unique(x$year))
  return(list.model)
}
#Iterate over all models
z1 <- lapply(Listyear,model_function)
#Export elements to envir
lapply(z1,list2env,.GlobalEnv)

你最终会得到 z1:

$`1975`
$`1975`$sub.mod1.1975
[1] "Threshold alpha: 85.0000277968913"
[1] ""
[1] "Model coefficients: Beta[0], Beta[1], Beta[2]"
(Intercept)           x           w 
  26.730070    3.376754   -2.406744 
      Change.Point Initial.Slope Slope.Change Second.Slope
2.5%      82.87297      3.259395    -2.515015    0.9283611
97.5%     87.90540      3.478656    -2.273062    1.0153773

$`1975`$sub.mod2.1975
[1] "Threshold alpha: 274.000071675723"
[1] ""
[1] "Model coefficients: Beta[0], Beta[1], Beta[2]"
(Intercept)           x           w 
 -37.968273    2.150220    5.115431 
      Change.Point Initial.Slope Slope.Change Second.Slope
2.5%      272.0000      1.969573     4.750341     7.057207
97.5%     276.0001      2.371539     5.468130     7.504963

在运行最后一行，您将获得全局环境中的模型：

希望对您有所帮助。

导出为csv的代码。

我添加了一个额外的功能，它从模型中获取一些结果并创建数据框，以便在对列表进行一些调整后可以轻松地将其导出到 .csv。接下来是函数：

model_export<-function(x)
{
  data <- x
  
  sub1 <- filter(data, julian.date < 210)
  sub2 <- filter(data, julian.date > 210)
  
  sub1.mod <- piecewise.linear(x= sub1$julian.date, y = sub1$rain_cs,
                               middle = 1,
                               CI = T,
                               bootstrap.samples = 1000)
  sub1.mod
  
  sub2.mod <- piecewise.linear(x= sub2$julian.date, y = sub2$rain_cs,
                               CI = T,
                               bootstrap.samples = 1000)
  sub2.mod
  
  #Group elements for models
  #Model 1
  modelname <- rep('sub1.mod',2)
  year <- rep(unique(x$year),2)
  changepoint <- rep(sub1.mod$change.point,2)
  coefs <- as.data.frame(t(sub1.mod$model$coefficients))
  intervals <- as.data.frame(sub1.mod$intervals)
  intervals <- cbind(data.frame(confidence=rownames(intervals)),intervals)
  rownames(intervals)<-NULL
  #Build DF
  DF1 <- data.frame(modelname,year,changepoint,coefs,intervals)
  #Model 2
  modelname <- rep('sub2.mod',2)
  changepoint <- rep(sub2.mod$change.point,2)
  coefs <- as.data.frame(t(sub2.mod$model$coefficients))
  intervals <- as.data.frame(sub2.mod$intervals)
  intervals <- cbind(data.frame(confidence=rownames(intervals)),intervals)
  rownames(intervals)<-NULL
  #Build DF
  DF2 <- data.frame(modelname,year,changepoint,coefs,intervals)
  #Bind DFs
  DFG <- rbind(DF1,DF2)
  return(DFG)
}

那么您可以申请：

#Apply new function to list
z2 <- lapply(Listyear,model_export)
#DF to export
MyDF <- do.call(rbind,z2)
#Export
write.csv(MyDF,file='Myfile.csv')

我已经使用它两年了，结果保存在 MyDF 中，然后导出到 .csv 文件。正如考虑到如果 rbind 因任何原因无法工作，您可以尝试 plyr 包中的 rbind.fill()。

应用多年的分段线性模型

Applying piecewise linear model for multiple year

r

tidyverse

piecewise