Introduction

Here we try susie on some example change point problems, and compare with other methods for change point detection in the changepoint package (penalized methods), bcp package (Bayesian MCMC method), genlasso (L1 penalty method), and L0learn (L0 penalty).

First we define some useful functions to run susie and other methods on changepoint problems and plot the CSs.

library("susieR")
library("genlasso")

Loading required package: Matrix

Loading required package: igraph

Attaching package: 'igraph'

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

    decompose, spectrum

The following object is masked from 'package:base':

    union

library("L0Learn")
library("bcp")

Loading required package: grid

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.

library("ggplot2")
library("DNAcopy")

# we set standardize to false because then the coefficients
# have a consistent interpretation as step size 
susie_cp = function(y,auto=FALSE,standardize=FALSE,...){
  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
    }
  }
  if(auto){
    s = susie_auto(X,y,standardize=standardize,...)
  } else {
    s = susie(X,y,standardize=standardize,...)
  }
  return(s)
}

#plot a time series y with confidence sets from susie fit s overlaid
# does - 0.5 so that singletons show up
# this is a ggplot version
susie_plot_cp = function(s,y){

  df<-data.frame(x = 1:length(y),y = y)
  CS = s$sets$cs

  p= ggplot(df) + geom_point(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
}

# this is just a function to add the changepoints to a base grapics plot
plot_cs = function(s){
  CS = s$sets$cs
  for(i in 1:length(CS)){
    rect(min(CS[[i]]),-5,max(CS[[i]])+1,5,col = rgb(1,0,0,alpha=0.5),border=NA)
  }
}

#ggplot function for changepoint results
plot_cp = function(df){
  #unwritten!
}

get_obj = function(s){
  return(s$elbo[length(s$elbo)])
}

This is a wrapper for L0learn

# wrapper to apply L0Learn to changepoint analysis
#coordinate ascent may not work so well so I use the slower/better algorithm, CDPSI
l0_cp = function(y,algorithm="CDPSI",maxSuppSize=20,...){ 
  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
    }
  }
  y.l0.cv = L0Learn.cvfit(X,y,nFolds=5,seed=1,penalty="L0",maxSuppSize = maxSuppSize,algorithm=algorithm,...) 
  opt = which.min(y.l0.cv$cvMeans[[1]])
  yhat = predict(y.l0.cv, newx = X,lambda=y.l0.cv$fit$lambda[[1]][opt])
  return(list(fit = y.l0.cv$fit,yhat=yhat))
}

Here is a wrapper for trendfiltering:

tf_cp = function(x){
  x.tf = trendfilter(x,ord=0)
  x.tf.cv = cv.trendfilter(x.tf)
  opt = which(x.tf$lambda==x.tf.cv$lambda.min) #optimal value of lambda
  yhat= x.tf$fit[,opt]
  return(list(fit=x.tf,yhat =yhat))
}

Here is a wrapper for segment from DNAcopy:

segment_cp = function(x){
  res = segment(CNA(x,rep(1,length(x)),1:length(x)))
  yhat = rep(res$output$seg.mean,diff(c(0,res$output$loc.end)))
  return(list(fit=res,yhat=yhat))
}

Simple simulated example

This example comes from Killick and Eckley

set.seed(10)
eg1=list()
eg1$x=c(rnorm(100,0,1),rnorm(100,1,1),rnorm(100,0,1),rnorm(100,0.2,1)) 
eg1$true_mean = c(rep(0,100),rep(1,100),rep(0,100),rep(0.2,100))

apply_methods = function(data){
  # Susie
  fit.s = susie_cp(data$x)
  yhat.s = predict(fit.s)
  
  # bcp
  fit.bcp = bcp(data$x)
  yhat.bcp = fit.bcp$posterior.mean
  
  # L0Learn
  res.l0= l0_cp(data$x)
  fit.l0 = res.l0$fit
  yhat.l0 = res.l0$yhat
  
  # trendfilter
  res.tf = tf_cp(data$x)
  fit.tf = res.tf$fit
  yhat.tf = res.tf$yhat
  
  # changepoint
  fit.cp = cpt.mean(data$x,method="PELT")
  d = diff(c(0,cpts(fit.cp),length(data$x))) 
  yhat.cp = rep(coef(fit.cp)$mean,d)
  
  #segment
  res.segment = segment_cp(data$x)
  fit.segment = res.segment$fit
  yhat.segment = res.segment$yhat
  
  return(list(fit = list(susie=fit.s,bcp=fit.bcp,l0=fit.l0,tf=fit.tf, cp = fit.cp, segment = fit.segment), yhat = list(susie=yhat.s,bcp = yhat.bcp,l0=yhat.l0,tf=yhat.tf,cp = yhat.cp, segment = yhat.segment)))
}

plot_results = function(res,data){
  plot(data$x,col="gray")
  lines(data$true_mean)
  for(i in 1:length(res$yhat)){
    lines(res$yhat[[i]],col=(i+1),type="s",lwd=2)
  }
}

compute_error = function(res,data){
  mse=rep(0,length(res$yhat))
  for(i in 1:length(res$yhat)){
    mse[i] = mean((res$yhat[[i]]-data$true_mean)^2)
  }
  names(mse) <- names(res$yhat)
  mse
}

eg1.res = apply_methods(eg1)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 
Analyzing: Sample.1 

plot_results(eg1.res,eg1)
plot_cs(eg1.res$fit$susie)

compute_error(eg1.res,eg1)

     susie        bcp         l0         tf         cp    segment 
0.03127789 0.03669960 0.22505358 0.04564223 0.03673359 0.03093324 

Compare PIPs of bcp and susie:

plot(eg1.res$fit$bcp$posterior.prob[-1], susie_get_PIP(eg1.res$fit$susie))

plot(eg1.res$fit$bcp$posterior.prob[-1],col=3)
points(susie_get_PIP(eg1.res$fit$susie),col=2)

Lai 2005 data

This is a real-data example from the changepoint package. Of course we do not know the truth here so cannot compute errors. But it is an interesting example because susie (run from default 0 initialization) misses a changepoint.

data(Lai2005fig4)
lai = list(x=Lai2005fig4[,5],true_mean = rep(NA,length(Lai2005fig4[,5])))

lai.res=apply_methods(lai)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 
Analyzing: Sample.1 

plot_results(lai.res,lai)
plot_cs(lai.res$fit$susie)

Here we try initializing susie from the results from changepoint and genlasso:

s0.cp = susie_init_coef(cpts(lai.res$fit$cp),diff(unlist(coef(lai.res$fit$cp))),length(lai$x)-1)
lai.s.cp0 = susie_cp(lai$x,s_init=s0.cp,estimate_prior_variance=FALSE)
lai.s.cp = susie_cp(lai$x,s_init=lai.s.cp0,estimate_prior_variance=TRUE)

bhat = lai.res$fit$tf$beta[,10] #result with 10 changepoints
dhat = diff(bhat)
dhat = ifelse(abs(dhat)>1e-8, dhat,0)
s0.tf = susie_init_coef(which(dhat!=0),dhat[dhat!=0],length(lai$x)-1)
lai.s.tf0 = susie_cp(lai$x,s_init=s0.tf,estimate_prior_variance=FALSE)
lai.s.tf = susie_cp(lai$x,s_init=lai.s.tf0,estimate_prior_variance=TRUE)

plot(lai$x)
lines(predict(lai.s.cp),col=2,type="s")
lines(predict(lai.s.tf),col=3,type="s")
get_obj(lai.s.cp)

[1] -216.9681

get_obj(lai.s.tf)

[1] -250.555

plot_cs(lai.s.cp)
plot_cs(lai.s.tf)

Simulation based on Lai et al

These data looked interesting because they seemed to be a bit challenging for some methods. But we do not know the truth of course. So I simulated some data based on the fit.

In the results here we see that trendfilter tends to include too many changepoints (not suprising); other methods produce similar results. The different behavior of bcp here vs the real data suggest to me that the real data may show non-gaussian residuals (eg that one outlier(?) point).

set.seed(1)
lai.mean = rep(unlist(coef(lai.res$fit$cp)),diff(c(0,cpts(lai.res$fit$cp),length(lai$x))))
lai.sd = sd(lai$x-lai.mean)
lai.sim = list(x=rnorm(length(lai.mean),lai.mean,lai.sd),true_mean=lai.mean)

lai.sim.res = apply_methods(lai.sim)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 
Analyzing: Sample.1 

plot_results(lai.sim.res,lai.sim)
plot_cs(lai.sim.res$fit$susie)

compute_error(lai.sim.res,lai.sim)

      susie         bcp          l0          tf          cp     segment 
0.315497051 0.015219952 0.060841408 0.036559321 0.008346073 0.270791108 

s0.cp = susie_init_coef(cpts(lai.sim.res$fit$cp),diff(unlist(coef(lai.sim.res$fit$cp))),length(lai.sim$x)-1)

lai.sim.s.cp = susie_cp(lai.sim$x,s_init=s0.cp,estimate_prior_variance=FALSE)
lai.sim.s.cp2 = susie_cp(lai.sim$x,s_init=lai.sim.s.cp,estimate_prior_variance=TRUE)

plot_results(lai.sim.res,lai.sim)
plot_cs(lai.sim.s.cp2)

mean((predict(lai.sim.s.cp2)-lai.sim$true_mean)^2)

[1] 0.008367509

Example from the BCP package

This one is described as a “hard” example (with one change point) in the bcp examples.

set.seed(5)
x <- rep(c(0,1), each=50)
eg2 = list(x = x +  rnorm(50, sd=1), true_mean = x)

eg2.res = apply_methods(eg2)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 
Analyzing: Sample.1 

plot_results(eg2.res,eg2)
plot_cs(eg2.res$fit$susie)

Warning in min(CS[[i]]): no non-missing arguments to min; returning Inf

Warning in max(CS[[i]]): no non-missing arguments to max; returning -Inf

Warning in min(CS[[i]]): no non-missing arguments to min; returning Inf

Warning in max(CS[[i]]): no non-missing arguments to max; returning -Inf

compute_error(eg2.res,eg2)

     susie        bcp         l0         tf         cp    segment 
0.05450532 0.13518020 0.28787237 0.35619464 0.01434871 0.26689895 

Try estimating prior variance of susie initialized from previous results; also try initializing from results of changepoint.

# estimate prior variance
eg2.s2 = susie_cp(eg2$x,s_init=eg2.res$fit$susie,estimate_prior_variance=TRUE)

# intialize from changepoint result
s0.cp = susie_init_coef(cpts(eg2.res$fit$cp),diff(unlist(coef(eg2.res$fit$cp))),length(eg2$x)-1)

eg2.s.cp = susie_cp(eg2$x,s_init=s0.cp,estimate_prior_variance=FALSE)
eg2.s.cp2 = susie_cp(eg2$x,s_init=eg2.s.cp,estimate_prior_variance=TRUE)


get_obj(eg2.s.cp2)

[1] -152.087

get_obj(eg2.s2)

[1] -152.6889

mean((predict(eg2.s.cp2)-eg2$true_mean)^2)

[1] 0.02936138

mean((predict(eg2.s2)-eg2$true_mean)^2)

[1] 0.0568752

DNA segmentation example from bcp package

This example comes from demo(coriell) in the bcp package.

data(coriell)
chrom11 <- list(x=as.vector(na.omit(coriell$Coriell.05296[coriell$Chromosome==11])))
chrom11$true_mean=rep(NA,length(chrom11$x))

chrom11.res = apply_methods(chrom11)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 
Analyzing: Sample.1 

plot_results(chrom11.res,chrom11)

Try susie. Note that this example illustrates a case where a variable (here 66) occurs in multiple CSs… something we don’t yet fully understand the implications of I think.

  plot(chrom11$x, col="grey", pch=20, xlab="Location",
        ylab="Posterior Mean",
        main="Coriell chromosome 11 (DNAcopy)")
  lines(predict(chrom11.res$fit$susie),col=2,lwd=2)
  chrom11.res$fit$susie$sets

$cs
$cs$L1
[1] 66

$cs$L2
[1] 51

$cs$L3
[1] 63 64 66 67 68 69 70 71


$purity
   min.abs.corr mean.abs.corr median.abs.corr
L1    1.0000000      1.000000       1.0000000
L2    1.0000000      1.000000       1.0000000
L3    0.9105706      0.965819       0.9659186

$cs_index
[1] 1 2 3

$coverage
[1] 0.95

  plot_cs(chrom11.res$fit$susie)

Compare different objectives:

An example from the DNAcopy segment function:

set.seed(51)
true_mean = rep(c(-0.2,0.1,1,-0.5,0.2,-0.5,0.1,-0.2),c(137,87,17,49,29,52,87,42))
genomdat = list(x = rnorm(500, sd=0.2) + true_mean, true_mean=true_mean)

genomdat.res = apply_methods(genomdat)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 
Analyzing: Sample.1 

plot_results(genomdat.res,genomdat)
plot_cs(genomdat.res$fit$susie)

compute_error(genomdat.res,genomdat)

       susie          bcp           l0           tf           cp 
0.0015949911 0.0009007416 0.0031717817 0.0018339170 0.0961235912 
     segment 
0.0014030174 

This example from DNAcopy too. (commented out for now as takes too long.)

# data(coriell)
# 
# #Combine into one CNA object to prepare for analysis on Chromosomes 1-23
# 
# CNA.object <-CNA(cbind(coriell$Coriell.05296,coriell$Coriell.13330),coriell$Chromosome,coriell$Position,data.type="logratio",sampleid=c("c05296","c13330"))
# 
# s = susie_cp(CNA.object$c13330[!is.na(CNA.object$c13330)])
# plot(CNA.object$c13330[!is.na(CNA.object$c13330)])
# plot_cs(s)

Compare initializations

Here we systematically compare intializing with trendfilter solution (L=20, and cross-validation optimum) with susie_auto and regular susie (run in two steps, first estimating residual variance, then estimating prior variance - this is because I am concerned estimating prior variance straight away may miss things, partly because at time of writing susie_init_coef may not be quite correct - it sets coefficients based on unstandardized X, but internally X are standardized)

#runs susie from both regular trend-filtering initialization
compare_init = function(x,fit.tf){
  bhat = fit.tf$beta[,20] #result with 20 changepoints
  dhat = diff(bhat)
  dhat = ifelse(abs(dhat)>1e-8, dhat,0)
  s0.tf20 = susie_init_coef(which(dhat!=0),dhat[dhat!=0],length(x)-1)
  s0.tf20.2 = susie_cp(x,s_init=s0.tf20,estimate_prior_variance=FALSE)
  s.tf20 = susie_cp(x,s_init=s0.tf20.2,estimate_prior_variance=TRUE)
  
  fit.tf.cv = cv.trendfilter(fit.tf)
  opt = which(fit.tf$lambda==fit.tf.cv$lambda.min) #optimal value of lambda
  bhat = fit.tf$beta[,opt] #result with 20 changepoints
  dhat = diff(bhat)
  dhat = ifelse(abs(dhat)>1e-8, dhat,0)
  s0.tf = susie_init_coef(which(dhat!=0),dhat[dhat!=0],length(x)-1)
  s0.tf.2 = susie_cp(x,s_init=s0.tf,estimate_prior_variance=FALSE)
  s.tf = susie_cp(x,s_init=s0.tf.2,estimate_prior_variance=TRUE)
  
  
  s0 = susie_cp(x,estimate_prior_variance = FALSE)
  s1 = susie_cp(x,estimate_prior_variance = TRUE)
  
  s_auto = susie_cp(x,auto=TRUE)
  return(list(fit = list(tf20=s.tf20, tf=s.tf,s=s1,s_auto=s_auto), obj = list(tf20 = get_obj(s.tf20), tf=get_obj(s.tf),s = get_obj(s1),s_auto= get_obj(s_auto))))
}

compare_init(eg1$x,eg1.res$fit$tf)$obj

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

$tf20
[1] -570.5557

$tf
[1] -570.5558

$s
[1] -570.5559

$s_auto
[1] -570.5637

compare_init(lai$x,lai.res$fit$tf)$obj

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

$tf20
[1] -226.4172

$tf
[1] -226.6029

$s
[1] -250.5503

$s_auto
[1] -217.0181

compare_init(lai.sim$x,lai.sim.res$fit$tf)$obj

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

$tf20
[1] -200.0482

$tf
[1] -200.0483

$s
[1] -248.9327

$s_auto
[1] -216.9614

compare_init(eg2$x,eg2.res$fit$tf)$obj

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

$tf20
[1] -152.6792

$tf
[1] -152.668

$s
[1] -152.087

$s_auto
[1] -152.087

c11 = compare_init(chrom11$x,chrom11.res$fit$tf)

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

c11$obj

$tf20
[1] 172.5228

$tf
[1] 172.523

$s
[1] 156.4556

$s_auto
[1] 172.5233

compare_init(genomdat$x,genomdat.res$fit$tf)$obj

Fold 1 ... Fold 2 ... Fold 3 ... Fold 4 ... Fold 5 ... 

$tf20
[1] 42.98781

$tf
[1] 42.88485

$s
[1] 45.9114

$s_auto
[1] 43.65838

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] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] DNAcopy_1.55.0    ggplot2_3.0.0     changepoint_2.2.2
 [4] zoo_1.8-4         bcp_4.0.3         L0Learn_1.0.7    
 [7] genlasso_1.4      igraph_1.2.2      Matrix_1.2-14    
[10] susieR_0.5.0.0347

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       yaml_2.2.0       
[19] bindrcpp_0.2.2    withr_2.1.2       stringr_1.3.1    
[22] dplyr_0.7.7       knitr_1.20        tidyselect_0.2.5 
[25] rprojroot_1.3-2   glue_1.3.0        R6_2.3.0         
[28] rmarkdown_1.10    reshape2_1.4.3    purrr_0.2.5      
[31] magrittr_1.5      whisker_0.3-2     backports_1.1.2  
[34] scales_1.0.0      htmltools_0.3.6   assertthat_0.2.0 
[37] colorspace_1.3-2  stringi_1.2.4     lazyeval_0.2.1   
[40] munsell_0.5.0     crayon_1.3.4      R.oo_1.22.0