changepoint

Introduction

Compare susier on some example change point problems and the changepoint package.

This example comes from Killick and Eckley

library(changepoint)

Loading required package: zoo

Attaching package: 'zoo'

The following objects are masked from 'package:base':

    as.Date, as.Date.numeric

Successfully loaded changepoint package version 2.2.2
 NOTE: Predefined penalty values changed in version 2.2.  Previous penalty values with a postfix 1 i.e. SIC1 are now without i.e. SIC and previous penalties without a postfix i.e. SIC are now with a postfix 0 i.e. SIC0. See NEWS and help files for further details.

set.seed(10)
m.data=c(rnorm(100,0,1),rnorm(100,1,1),rnorm(100,0,1),rnorm(100,0.2,1)) 
ts.plot(m.data,xlab="Index")

Here we apply susier to this example. It finds 2 (out of the three) changepoints.

library("susieR")
susie_cp = function(y,...){
  n=length(y)
  X = matrix(0,nrow=n,ncol=n-1)
  for(j in 1:(n-1)){
    for(i in (j+1):n){
      X[i,j] = 1
    }
  }
  s = susie(X,y,min_abs_corr=0.9,...)
  return(s)
}


s = susie_cp(m.data)
ts.plot(m.data,xlab="Index")
lines(predict(s),col=2,lwd=2)

Note that susieR can provide credible sets for the changepoints. Here it finds two, shown in this plot:

#plot a time series y with confidence sets from susie fit s overlaid
# does +- 0.5 so that singletons show up
plot_cp = function(s,y){
  library("ggplot2")
  df<-data.frame(x = 1:length(y),y = y)
  CS = s$sets$cs

  p= ggplot(df) + geom_line(mapping=aes_string(x="x", y="y"))
  for(i in 1:length(CS)){
    p = p+  annotate("rect", fill = "red", alpha = 0.5, 
        xmin = min(CS[[i]])-0.5, xmax = max(CS[[i]])+0.5,
        ymin = -Inf, ymax = Inf) 
  }
  p
}
plot_cp(s,m.data)

data(Lai2005fig4)
Lai.default=cpt.mean(Lai2005fig4[,5],method="PELT")
plot(Lai.default,pch=20,col="grey",cpt.col="black",type="p",xlab="Index") 

cpts(Lai.default)

[1]  81  85  89  96 123 133

coef(Lai.default)

$mean
[1] 0.2468910 4.6699210 0.4495538 4.5902489 0.2079891 4.2913844 0.2291286

We see in this case susie seems to “miss” one of the changepoints.

fitted = susie_cp(Lai2005fig4[,5])
plot(Lai2005fig4[,5])
lines(predict(fitted),col=2)

See if this is maybe a convergence issue solved by susie_auto:

susie_cp_auto = function(y,...){
  n=length(y)
  X = matrix(0,nrow=n,ncol=n-1)
  for(j in 1:(n-1)){
    for(i in (j+1):n){
      X[i,j] = 1
    }
  }
  s = susie_auto(X,y,min_abs_corr=0.9,...)
  return(s)
}

fitted_a = susie_cp_auto(Lai2005fig4[,5])
plot(Lai2005fig4[,5])
lines(predict(fitted_a),col=2)

Plot confidence sets:

plot_cp(fitted_a,Lai2005fig4[,5])

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ggplot2_3.0.0     susieR_0.5.0.0347 changepoint_2.2.2 zoo_1.8-4        

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.19       bindr_0.1.1        compiler_3.5.1    
 [4] pillar_1.3.0       git2r_0.23.0       plyr_1.8.4        
 [7] workflowr_1.1.1    R.methodsS3_1.7.1  R.utils_2.7.0     
[10] tools_3.5.1        digest_0.6.18      evaluate_0.12     
[13] tibble_1.4.2       gtable_0.2.0       lattice_0.20-35   
[16] pkgconfig_2.0.2    rlang_0.2.2        Matrix_1.2-14     
[19] yaml_2.2.0         bindrcpp_0.2.2     withr_2.1.2       
[22] stringr_1.3.1      dplyr_0.7.7        knitr_1.20        
[25] tidyselect_0.2.5   rprojroot_1.3-2    grid_3.5.1        
[28] glue_1.3.0         R6_2.3.0           rmarkdown_1.10    
[31] purrr_0.2.5        magrittr_1.5       whisker_0.3-2     
[34] backports_1.1.2    scales_1.0.0       matrixStats_0.54.0
[37] htmltools_0.3.6    assertthat_0.2.0   colorspace_1.3-2  
[40] labeling_0.3       stringi_1.2.4      lazyeval_0.2.1    
[43] munsell_0.5.0      crayon_1.3.4       R.oo_1.22.0