Last updated: 2019-12-20

Checks: 7 0

Knit directory: smash-paper/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


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.

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.

The command set.seed(1) 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.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use 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:    dsc/code/Wavelab850/MEXSource/CPAnalysis.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/DownDyadHi.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/DownDyadLo.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FAIPT.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FCPSynthesis.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FMIPT.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FWPSynthesis.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FWT2_PO.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FWT_PBS.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FWT_PO.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/FWT_TI.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/IAIPT.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/IMIPT.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/IWT2_PO.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/IWT_PBS.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/IWT_PO.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/IWT_TI.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/LMIRefineSeq.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/MedRefineSeq.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/UpDyadHi.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/UpDyadLo.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/WPAnalysis.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/dct_ii.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/dct_iii.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/dct_iv.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/dst_ii.mexmac
    Ignored:    dsc/code/Wavelab850/MEXSource/dst_iii.mexmac

Untracked files:
    Untracked:  analysis/temp.R
    Untracked:  files.txt

Unstaged changes:
    Modified:   analysis/gausfullres.Rmd
    Modified:   code/gaussmeanest.functions.R

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.


These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File Version Author Date Message
Rmd 02f9167 Peter Carbonetto 2019-12-20 wflow_publish(“gaussmeanest.Rmd”)
html c8548d5 Peter Carbonetto 2019-12-20 Re-built gaussmeanest analysis after making some of the plotting code
Rmd 4d7aae3 Peter Carbonetto 2019-12-20 wflow_publish(“gaussmeanest.Rmd”)
html 8998bb8 Peter Carbonetto 2019-12-20 Re-built gaussmeanest analysis after Zhengrong’s recent additions.
Rmd 85aa594 Peter Carbonetto 2019-12-20 wflow_publish(“gaussmeanest.Rmd”)
Rmd f0221c5 Zhengrong Xing 2019-10-27 address some reviewer comments
html f0221c5 Zhengrong Xing 2019-10-27 address some reviewer comments
html 74aff51 Peter Carbonetto 2018-12-20 Re-built gauss_shiny_setup and gaussmeanest pages.
Rmd 86da808 Peter Carbonetto 2018-12-10 Added pointers to dsc/README in relevant R Markdown files.
html 8caff70 Peter Carbonetto 2018-12-06 Re-built the workflowr pages after several minor changes to the text
Rmd c589dbb Peter Carbonetto 2018-12-06 wflow_publish(c(“index.Rmd”, “gaussian_signals.Rmd”,
html ee71f27 Peter Carbonetto 2018-12-04 Made a few small adjustments to the text in the “gaussianmeanest” analysis.
Rmd eb6cc34 Peter Carbonetto 2018-12-04 wflow_publish(“gaussmeanest.Rmd”)
html 05684ba Peter Carbonetto 2018-12-04 Ran wflow_publish(“gaussmeanest.Rmd”) to populate the webpage.
Rmd 9a67b48 Peter Carbonetto 2018-12-02 Moved dsc results file.
Rmd 049dcbb Peter Carbonetto 2018-11-08 Moved around some files and revised TOC in home page.

In this analysis, we assess the ability of different signal denoising methods to recover the true signal after being provided with Gaussian-distributed observations of the signal. We consider scenarios in which the data have homoskedastic errors (constant variance) and heteroskedastic errors (non-constant variance).

Since the simulation experiments are computationally intensive, here we only illustrate the application of the signal denoising methods, and create plots summarizing the results of the full experiments; the full experiments were implemented separately. (For instructions on re-running these simulation experiments, see the README in the “dsc” directory of the git repository).

Set up environment

Load the ggplot2 and cowplot packages, and the functions definining the mean and variances used to simulate the data.

library(plyr)
library(smashr)
library(ggplot2)
library(cowplot)
source("../code/signals.R")
source("../code/gaussmeanest.functions.R")

Load results

Load the results of the simulation experiments.

load("../output/dscr.RData")

Simulated data with constant variances

This plot reproduces Fig. 2 of the manuscript, which compares the accuracy of the mean curves estimated from the data sets that were simulated using the “Spikes” mean function with constant variance and a signal-to-noise ratio of 3.

First, extract the results used to generate this plot, and transform them into a data frame suitable for plotting using ggplot2.

pdat <- get.results.homosked(res,"sp.3.v1")

Create the combined boxplot and violin plot using ggplot2.

pdat <-
  transform(pdat,
            method = factor(method,
                            names(sort(tapply(pdat$mise,pdat$method,mean),
                                       decreasing = TRUE))))
p <- ggplot(pdat,aes(x = method,y = mise,fill = method.type)) +
     geom_violin(fill = "skyblue",color = "skyblue") +
     geom_boxplot(width = 0.15,outlier.shape = NA) +
     scale_y_continuous(breaks = seq(6,16,2)) +
     scale_fill_manual(values = c("darkorange","dodgerblue","gold"),
                       guide = FALSE) +
     coord_flip() +
     labs(x = "",y = "MISE") +
     theme(axis.line = element_blank(),
           axis.ticks.y = element_blank())
print(p)

From this plot, we see that the three variants of SMASH all outperformed EbayesThresh and TI thresholding in this setting.

These plots summarize the results for all 7 simulation scenarios and the two signal-to-noise ratios (1 and 3), including the “Spikes” scenario shown in greater detail in the violin plot above.

create.bar.plots(rbind(
  data.frame(sim = "sp",  snr = 1,get.results.homosked(res,"sp.1.v1")),
  data.frame(sim = "bump",snr = 1,get.results.homosked(res,"bump.1.v1")),
  data.frame(sim = "blkp",snr = 1,get.results.homosked(res,"blk.1.v1")),
  data.frame(sim = "ang", snr = 1,get.results.homosked(res,"ang.1.v1")),
  data.frame(sim = "dop", snr = 1,get.results.homosked(res,"dop.1.v1")),
  data.frame(sim = "blip",snr = 1,get.results.homosked(res,"blip.1.v1")),
  data.frame(sim = "cor", snr = 1,get.results.homosked(res,"cor.1.v1")),
  data.frame(sim = "sp",  snr = 3,get.results.homosked(res,"sp.3.v1")),
  data.frame(sim = "bump",snr = 3,get.results.homosked(res,"bump.3.v1")),
  data.frame(sim = "blkp",snr = 3,get.results.homosked(res,"blk.3.v1")),
  data.frame(sim = "ang", snr = 3,get.results.homosked(res,"ang.3.v1")),
  data.frame(sim = "dop", snr = 3,get.results.homosked(res,"dop.3.v1")),
  data.frame(sim = "blip",snr = 3,get.results.homosked(res,"blip.3.v1")),
  data.frame(sim = "cor", snr = 3,get.results.homosked(res,"cor.3.v1"))))

Next, we compare the same methods in simulated data sets with heteroskedastic errors.

Simulated data with heteroskedastic errors: “Spikes” mean signal and “Clipped Blocks” variance

In this scenario, the data sets were simulated using the “Spikes” mean function and the “Clipped Blocks” variance function. The next two plots reproduce part of Fig. 3 in the manuscript.

This plot shows the mean function as a black line, and the +/- 2 standard deviations as orange lines:

t         <- (1:1024)/1024
mu        <- spike.fn(t,"mean")
sigma.ini <- sqrt(cblocks.fn(t,"var"))
sd.fn     <- sigma.ini/mean(sigma.ini) * sd(mu)/3
par(cex.axis = 1,cex.lab = 1.25)
plot(mu,type = "l", ylim = c(-0.05,1),xlab = "position",ylab = "",
     lwd = 1.75,xaxp = c(0,1024,4),yaxp = c(0,1,4))
lines(mu + 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)
lines(mu - 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)

Version Author Date
8998bb8 Peter Carbonetto 2019-12-20
f0221c5 Zhengrong Xing 2019-10-27
05684ba Peter Carbonetto 2018-12-04

Now, we plot the ground-truth signal (the mean function, drawn as a black line) and the signals recovered by TI thresholding (light blue line) and SMASH (the red line) for one simulated dataset as an illustration

t         <- (1:1024)/1024
mu        <- spike.fn(t,"mean")
sigma.ini <- sqrt(cblocks.fn(t,"var"))
sd.fn     <- sigma.ini/mean(sigma.ini) * sd(mu)/3
x.sim  <- rnorm(1024,mu,sd.fn)
mu.smash <- smash(x.sim,family = "DaubLeAsymm",filter.number = 8)
mu.ti    <- ti.thresh(x.sim,method = "rmad",family = "DaubLeAsymm",
                      filter.number = 8)
par(cex.axis = 1)
plot(mu,type = "l",col = "black",lwd = 3,xlab = "position",ylab = "",
     ylim = c(-0.05,1),xaxp = c(0,1024,4),yaxp = c(0,1,4))
lines(mu.ti,col = "dodgerblue",lwd = 3)
lines(mu.smash,col = "orangered",lwd = 3)

Version Author Date
8998bb8 Peter Carbonetto 2019-12-20

Extract the results from running the simulations.

hetero.data.smash <-
  res[res$.id == "sp.3.v5" & res$method == "smash.s8",]
hetero.data.smash.homo <-
  res[res$.id == "sp.3.v5" & res$method == "smash.homo.s8",]
hetero.data.tithresh.homo <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.homo.s8",]
hetero.data.tithresh.rmad <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.rmad.s8",]
hetero.data.tithresh.smash <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.smash.s8",]
hetero.data.tithresh.true <-
  res[res$.id == "sp.3.v5" & res$method == "tithresh.true.s8",]
hetero.data.ebayes <-
  res[res$.id == "sp.3.v5" & res$method == "ebayesthresh",]
hetero.data.smash.true <-
  res[res$.id == "sp.3.v5" & res$method == "smash.true.s8",]

Transform these results into a data frame suitable for ggplot2.

pdat <-
  rbind(data.frame(method      = "smash",
                   method.type = "est",
                   mise        = hetero.data.smash$mise),
        data.frame(method      = "smash.homo",
                   method.type = "homo",
                   mise        = hetero.data.smash.homo$mise),
        data.frame(method      = "tithresh.rmad",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.rmad$mise),
        data.frame(method      = "tithresh.smash",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.smash$mise),
        data.frame(method      = "tithresh.true",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.true$mise),
        data.frame(method      = "ebayesthresh",
                   method.type = "homo",
                   mise        = hetero.data.ebayes$mise),
        data.frame(method      = "smash.true",
                   method.type = "true",
                   mise        = hetero.data.smash.true$mise))
pdat <-
  transform(pdat,
            method = factor(method,
                            names(sort(tapply(pdat$mise,pdat$method,mean),
                                       decreasing = TRUE))))

Create the combined boxplot and violin plot using ggplot2.

p <- ggplot(pdat,aes(x = method,y = mise,fill = method.type)) +
     geom_violin(fill = "skyblue",color = "skyblue") +
     geom_boxplot(width = 0.15,outlier.shape = NA) +
     scale_fill_manual(values=c("darkorange","dodgerblue","limegreen","gold"),
                       guide = FALSE) +
     coord_flip() +
     scale_y_continuous(breaks = seq(10,70,10)) +
     labs(x = "",y = "MISE") +
     theme(axis.line = element_blank(),
           axis.ticks.y = element_blank())
print(p)

Version Author Date
8998bb8 Peter Carbonetto 2019-12-20
f0221c5 Zhengrong Xing 2019-10-27
05684ba Peter Carbonetto 2018-12-04

In this scenario, we see that SMASH, when allowing for heteroskedastic errors, outperforms EbayesThresh and all variants of TI thresholding (including TI thresholding when provided with the true variance). Further, SMASH performs almost as well when estimating the variance compared to when provided with the true variance.

Simulated data with heteroskedastic errors: “Corner” mean signal and “Doppler” variance

In this next scenario, the data sets were simulated using the “Corner” mean function and the “Doppler” variance function. These plots were also used in Fig. 3 of the manuscript.

This plot shows the mean function as a black line, and the +/- 2 standard deviations as orange lines:

mu        <- cor.fn(t,"mean") 
sigma.ini <- sqrt(doppler.fn(t,"var"))
sd.fn     <- sigma.ini/mean(sigma.ini) * sd(mu)/3
plot(mu,type = "l", ylim = c(-0.05,1),xlab = "position",ylab = "",
     lwd = 1.75,xaxp = c(0,1024,4),yaxp = c(0,1,4))
lines(mu + 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)
lines(mu - 2*sd.fn,col = "darkorange",lty = 5,lwd = 1.75)

Version Author Date
8998bb8 Peter Carbonetto 2019-12-20
f0221c5 Zhengrong Xing 2019-10-27
05684ba Peter Carbonetto 2018-12-04

Now, we plot the ground-truth signal (the mean function, drawn as a black line) and the signals recovered by TI thresholding (light blue line) and SMASH (the red line) for one simulated dataset as an illustration

t         <- (1:1024)/1024
mu        <- cor.fn(t,"mean") 
sigma.ini <- sqrt(doppler.fn(t,"var"))
sd.fn     <- sigma.ini/mean(sigma.ini) * sd(mu)/3
x.sim  <- rnorm(1024,mu,sd.fn)
mu.smash <- smash(x.sim,family = "DaubLeAsymm",filter.number = 8)
mu.ti    <- ti.thresh(x.sim,method = "rmad",family = "DaubLeAsymm",
                      filter.number = 8)
par(cex.axis = 1)
plot(mu,type = "l",col = "black",lwd = 3,xlab = "position",ylab = "",
     ylim = c(-0.05,1),xaxp = c(0,1024,4),yaxp = c(0,1,4))
lines(mu.ti,col = "dodgerblue",lwd = 3)
lines(mu.smash,col = "orangered",lwd = 3)

Version Author Date
8998bb8 Peter Carbonetto 2019-12-20

Extract the results from running these simulations.

hetero.data.smash.2 <-
  res[res$.id == "cor.3.v3" & res$method == "smash.s8",]
hetero.data.smash.homo.2 <-
  res[res$.id == "cor.3.v3" & res$method == "smash.homo.s8",]
hetero.data.tithresh.homo.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.homo.s8",]
hetero.data.tithresh.rmad.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.rmad.s8",]
hetero.data.tithresh.smash.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.smash.s8",]
hetero.data.tithresh.true.2 <-
  res[res$.id == "cor.3.v3" & res$method == "tithresh.true.s8",]
hetero.data.ebayes.2 <-
  res[res$.id == "cor.3.v3" & res$method == "ebayesthresh",]
hetero.data.smash.true.2 <-
  res[res$.id == "cor.3.v3" & res$method == "smash.true.s8",]

Transform these results into a data frame suitable for ggplot2.

pdat <-
  rbind(data.frame(method      = "smash",
                   method.type = "est",
                   mise        = hetero.data.smash.2$mise),
        data.frame(method      = "smash.homo",
                   method.type = "homo",
                   mise        = hetero.data.smash.homo.2$mise),
        data.frame(method      = "tithresh.rmad",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.rmad.2$mise),
        data.frame(method      = "tithresh.smash",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.smash.2$mise),
        data.frame(method      = "tithresh.true",
                   method.type = "tithresh",
                   mise        = hetero.data.tithresh.true.2$mise),
        data.frame(method      = "ebayesthresh",
                   method.type = "homo",
                   mise        = hetero.data.ebayes.2$mise),
        data.frame(method      = "smash.true",
                   method.type = "true",
                   mise        = hetero.data.smash.true.2$mise))
pdat <-
  transform(pdat,
            method = factor(method,
                            names(sort(tapply(pdat$mise,pdat$method,mean),
                                       decreasing = TRUE))))

Create the combined boxplot and violin plot using ggplot2.

p <- ggplot(pdat,aes(x = method,y = mise,fill = method.type)) +
     geom_violin(fill = "skyblue",color = "skyblue") +
     geom_boxplot(width = 0.15,outlier.shape = NA) +
     scale_fill_manual(values=c("darkorange","dodgerblue","limegreen","gold"),
                       guide = FALSE) +
     coord_flip() +
     scale_y_continuous(breaks = seq(1,5)) +
     labs(x = "",y = "MISE") +
     theme(axis.line = element_blank(),
           axis.ticks.y = element_blank())
print(p)

Version Author Date
8998bb8 Peter Carbonetto 2019-12-20
f0221c5 Zhengrong Xing 2019-10-27
05684ba Peter Carbonetto 2018-12-04

Similar to the “Spikes” scenario, we see that the SMASH method, when allowing for heteroskedastic variances, outperforms both the TI thresholding and EbayesThresh approaches.


sessionInfo()
# R version 3.4.3 (2017-11-30)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.6
# 
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.4/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] cowplot_0.9.4 ggplot2_3.2.0 smashr_1.2-5  plyr_1.8.4   
# 
# loaded via a namespace (and not attached):
#  [1] wavethresh_4.6.8  tidyselect_0.2.5  xfun_0.7         
#  [4] reshape2_1.4.3    ashr_2.2-39       purrr_0.2.5      
#  [7] lattice_0.20-35   colorspace_1.4-0  htmltools_0.3.6  
# [10] yaml_2.2.0        rlang_0.4.2       mixsqp_0.3-9     
# [13] later_0.8.0       pillar_1.3.1      glue_1.3.1       
# [16] withr_2.1.2.9000  dscr_0.1-8        foreach_1.4.4    
# [19] stringr_1.4.0     munsell_0.4.3     gtable_0.2.0     
# [22] workflowr_1.6.0   caTools_1.17.1.2  codetools_0.2-15 
# [25] evaluate_0.13     labeling_0.3      knitr_1.23       
# [28] pscl_1.5.2        doParallel_1.0.14 httpuv_1.5.0     
# [31] parallel_3.4.3    Rcpp_1.0.1        xtable_1.8-2     
# [34] promises_1.0.1    backports_1.1.2   scales_0.5.0     
# [37] truncnorm_1.0-8   mime_0.6          fs_1.2.7         
# [40] digest_0.6.18     stringi_1.4.3     dplyr_0.8.0.1    
# [43] shiny_1.2.0       grid_3.4.3        rprojroot_1.3-2  
# [46] tools_3.4.3       bitops_1.0-6      magrittr_1.5     
# [49] lazyeval_0.2.1    tibble_2.1.1      crayon_1.3.4     
# [52] whisker_0.3-2     pkgconfig_2.0.2   MASS_7.3-48      
# [55] Matrix_1.2-12     SQUAREM_2017.10-1 data.table_1.12.0
# [58] assertthat_0.2.1  rmarkdown_2.0     iterators_1.0.10 
# [61] R6_2.4.0          git2r_0.26.1      compiler_3.4.3