segmented.R

library("lubridate")
library("strucchange")
library("segmented")

set.seed(1234567)
N <- 60
df  <- data.frame(id  =  1:N)
df$date  <- seq(as.Date("2013-07-01"), by = "day", along = df$id)
df$date2  <-  difftime(df$date, ymd("2013-07-01"), units = "day")
df$date3  <-  ifelse(df$date > as.Date("2013-08-01"), difftime(df$date, ymd("2013-08-01"), units = "day"), 0)
difftime(ymd("2013-08-01"), ymd("2013-07-01"), units = "day")
df$u  <- 1 + 10 * rnorm(N)
alpha  <- .5
df$y  <- ifelse(df$date > as.Date("2013-08-01"), alpha * difftime(ymd("2013-08-01"), ymd("2013-07-01"), units = "day") + - 2 * as.numeric(df$date3) + df$u, alpha * as.numeric(df$date2) + df$u)

bp  <- breakpoints(y ~ date, data = df, breaks = 1)
df$fit  <- fitted(bp)

df$ideal  <- ifelse(df$date > as.Date("2013-08-01"), alpha * difftime(ymd("2013-08-01"), ymd("2013-07-01"), units = "day") + - 2 * as.numeric(df$date3), alpha * as.numeric(df$date2))

m  <- lm(y ~ date2, data = df)
seg  <- segmented(m, seg.Z= ~ date2, psi = list(date2 = c(10)), control=seg.control(display=FALSE))
df$seg  <- fitted(seg)

ggplot(data = df, aes(x = date, y = y)) + 
  theme_bw() + 
  geom_point() + 
  geom_line(aes(y = fit), colour = "blue") + 
  geom_line(aes(y = ideal), colour = "red") + 
  geom_line(aes(y = seg), colour = "green")