CODEWITHSUNDEEP
rtidyversepiecewise
UseR10085
UseR10085

Reputation: 8198

Applying piecewise linear model for multiple year

I have daily rainfall data which I have converted to yearwise cumulative value using following code

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

Then I want to divide every year into 2 parts (before 210 days and after 210 days) then apply piecewise linear model from SiZer to identify yearwise breakpoints. I could able to do it for single year like

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

Now how to dynamically fit piecewise linear model for all the years?

Upvotes: 1

Views: 122

Answers (1)

Duck
Duck

Reputation: 39613

You can try using a function and base R, creating a list and then saving the models. I include in last line a way to export all models outside the list:

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)

You will end up with 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

And by running last line you will get the models in the global environment:

enter image description here

I hope this can help.

Code for exporting to csv.

I include an additional function that takes some results from the models and creates dataframes so that it can be easily exported to .csv after doing some adjusts to lists. The function is next:

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)
}

Then you can apply:

#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')

I have used it for two years having the results saved in MyDF and then exported to .csv file. Just as consideration if rbind would not work for any reason you could try rbind.fill() from plyr package.

Upvotes: 1

Related Questions