Last updated: 2018-10-23
workflowr checks: (Click a bullet for more information) ✔ R Markdown file: up-to-date
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
✔ Environment: empty
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
✔ Seed:
set.seed(20180414)
The command set.seed(20180414)
was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
✔ Repository version: 0c1559e
wflow_publish
or wflow_git_commit
). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.Rhistory
Untracked files:
Untracked: analysis/null.Rmd
Untracked: analysis/test.Rmd
Untracked: data/geneMatrix.tsv
Untracked: data/liter_data_4_summarize_ld_1_lm_less_3.rds
Untracked: data/meta.tsv
Untracked: docs/figure/test.Rmd/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
Here I investigate using trendfilter to initialize the susie fit for changepoint problems.
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")
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)
}
The goal here is to define an automatic version of susie for changepoint problems based on initializing from a trendfilter fit.
# runs susie initialized to fit from trendfilter with L changepoints
susie_cp_tf = function(y,fit.tf,L){
bhat = fit.tf$beta[,L] #result with L changepoints
dhat = diff(bhat)
dhat = ifelse(abs(dhat)>1e-8, dhat,0)
s0.tf = susie_init_coef(which(dhat!=0),dhat[dhat!=0],length(y)-1)
print(which(dhat!=0))
print(dhat[dhat!=0])
s0.tf.2 = susie_cp(y,s_init=s0.tf,estimate_prior_variance=FALSE)
s = susie_cp(y,s_init=s0.tf.2,estimate_prior_variance=TRUE)
return(s)
}
get_obj = function(s){
return(s$elbo[length(s$elbo)])
}
# This function simply runs susie for L=1,2,4,8. If the best is not L=8 then
# it returns the best of these. Otherwise it keeps doubling L until the solution gets worse.
susie_cp_auto_tf = function(y,fit.tf){
s=list()
s[[1]] = susie_cp_tf(y,fit.tf,1+1) #+1 because first solution is constant (no breakpoints)
s[[2]] = susie_cp_tf(y,fit.tf,2+1)
s[[3]] = susie_cp_tf(y,fit.tf,4+1)
s[[4]] = susie_cp_tf(y,fit.tf,8+1)
obj = lapply(s,get_obj)
opt = which.max(obj)
if(opt!=4){
return(s[[opt]])
} else {
L=8
s_old = s[[4]]
repeat{
s_new = susie_cp_tf(y,fit.tf,2*L+1)
if(get_obj(s_new) < get_obj(s_old)){break;}
L = L*2
s_old = s_new
}
}
return(s_old)
}
data(Lai2005fig4)
lai = list(x=Lai2005fig4[,5],true_mean = rep(NA,length(Lai2005fig4[,5])))
lai.tf = trendfilter(lai$x,ord=0)
s= susie_cp_auto_tf(lai$x, lai.tf)
[1] 81
[1] 0.1167276
[1] 81 133
[1] 0.8614129 -0.8080183
[1] 81 96 123 133
[1] 1.2501386 -0.4372270 0.3013434 -1.0723904
[1] 76 81 89 96 97 122 123 133
[1] 0.01506203 2.23227569 0.49293550 -2.35436338 -0.12029902 0.65215478
[7] 1.98491595 -2.89626068
[1] 25 33 49 72 76 81 85 89 96 97 98 103 122 123 125 133
[1] 0.21349567 -0.05480517 -0.29508482 0.10029882 0.10605781
[6] 3.25513106 -2.04571074 2.43203643 -3.27060894 -0.28135461
[11] -0.09663054 -0.00358927 0.91343018 1.98491595 0.80171487
[16] -3.75142321
[1] 11 25 26 28 32 33 49 53 54 55 72 76 81 85 89 95 96
[18] 97 98 103 111 115 122 123 125 128 133 147 154 159 176 187
[1] 0.041095098 0.234687244 0.208869253 0.113288205 -0.064854824
[6] -0.312107804 -0.484548824 -0.268979866 0.211009797 0.218026009
[11] 0.129706880 0.106057814 3.746048802 -3.027546224 3.246206677
[16] -0.299096566 -3.294764884 -0.281354613 -0.096630542 -0.155801964
[21] 0.070944289 0.021140435 0.973558155 1.984915946 0.916096260
[26] 0.209723965 -4.152694944 0.043316043 0.003419216 0.089376319
[31] -0.048950754 -0.114062027
get_obj(s)
[1] -226.2105
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
This reproducible R Markdown analysis was created with workflowr 1.1.1