Goal

Power analysis of low-dimensional covariance regression model. Run this script after Low_Dimensional.Rmd.

library(CovRegFC)
## Loading required package: Rcpp
## Warning: replacing previous import 'cowplot::ggsave' by 'ggplot2::ggsave'
## when loading 'CovRegFC'
library(ggplot2)
library(magrittr)
library(reshape2)
library(devtools)
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggplot2':
## 
##     ggsave
params
## $num_regions
## [1] "15"

Power Analysis

num_regions = as.integer(params$num_regions)
num_regions
## [1] 15
channel_names = paste0("R",1:num_regions)
num_samples_seq = seq(20,160,20)
load("fit_list_low_15_NA_167.Rdata")
fit_population = fit_list[[1]]
coeff = CovRegFC::plot_coeff(fit_population,
                             response = channel_names,
                             alpha = 0.05/length(channel_names))
compute_power = function(alpha) {
  df = lapply(num_samples_seq,function(n) {
    load(paste0("fit_list_low_15_",n,"_167.Rdata"))
    stats = lapply(fit_list,function(fit) {
      sc = CovRegFC::plot_coeff(fit,
                                response = channel_names,
                                alpha = alpha)
      # True Positive Rate (TPR)
      P = length(unlist(coeff))
      TP1 = sum(sc$set1 %in% coeff$set1) + sum(sc$set2 %in% coeff$set2)
      TP2 = sum(sc$set1 %in% coeff$set2) + sum(sc$set2 %in% coeff$set1)
      TPR = max(TP1,TP2)/P
      # False Discovery Rate (FDR)
      FP1 = sum(!sc$set1 %in% coeff$set1) + sum(!sc$set2 %in% coeff$set2)
      FP2 = sum(!sc$set1 %in% coeff$set2) + sum(!sc$set2 %in% coeff$set1)
      FDR = min(FP1,FP2)/(max(TP1,TP2)+min(FP1,FP2))
      if(is.na(FDR)) FDR = 0
      c(TPR = TPR,FDR = FDR)
    }) %>% do.call(rbind,.) %>% data.frame
    c(n = n,TPR = mean(stats$TPR),FDR = mean(stats$FDR))
  }) %>% do.call(rbind,.) %>% data.frame
  data.frame(df,alpha)
}
alpha_seq = c(0.001,0.01,0.1)
df_list = lapply(alpha_seq,compute_power)
df = df_list %>% do.call(rbind,.)
df$alpha = as.factor(df$alpha)
df_long = melt(df,id.vars = c("n","alpha"),value.name = "rate",variable.name = "statistic")
ggplot(df_long,aes(x = n,y = rate,linetype = statistic,color = alpha)) + 
  geom_point() + 
  geom_line() +
  scale_y_continuous(limits = c(0,1),breaks = seq(0,1,0.2)) + 
  scale_x_continuous(breaks = df$n) +
  ggtitle("Power Analysis")

ggsave(filename = "power.pdf",width = 5,height = 5)

Session Info

session_info()
## Session info -------------------------------------------------------------
##  setting  value                       
##  version  R version 3.3.2 (2016-10-31)
##  system   x86_64, darwin13.4.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  tz       America/Los_Angeles         
##  date     2017-06-22
## Packages -----------------------------------------------------------------
##  package     * version    date       source                          
##  backports     1.1.0      2017-05-22 CRAN (R 3.3.2)                  
##  base        * 3.3.2      2016-10-31 local                           
##  colorspace    1.3-2      2016-12-14 CRAN (R 3.3.2)                  
##  CovRegFC    * 0.1.0      2017-06-22 local                           
##  cowplot     * 0.7.0      2016-10-28 CRAN (R 3.3.0)                  
##  datasets    * 3.3.2      2016-10-31 local                           
##  devtools    * 1.13.2     2017-06-02 CRAN (R 3.3.2)                  
##  digest        0.6.12     2017-01-27 cran (@0.6.12)                  
##  evaluate      0.10       2016-10-11 CRAN (R 3.2.3)                  
##  ggplot2     * 2.2.1      2016-12-30 CRAN (R 3.3.2)                  
##  graphics    * 3.3.2      2016-10-31 local                           
##  grDevices   * 3.3.2      2016-10-31 local                           
##  grid          3.3.2      2016-10-31 local                           
##  gridExtra     2.2.1      2016-02-29 CRAN (R 3.2.4)                  
##  gtable        0.2.0      2016-02-26 CRAN (R 3.2.3)                  
##  htmltools     0.3.6      2017-04-28 cran (@0.3.6)                   
##  inline        0.3.14     2015-04-13 CRAN (R 3.2.0)                  
##  knitr         1.16       2017-05-18 cran (@1.16)                    
##  lazyeval      0.2.0.9000 2017-05-09 Github (hadley/lazyeval@c155c3d)
##  magrittr    * 1.5        2014-11-22 CRAN (R 3.2.0)                  
##  memoise       1.1.0      2017-04-21 CRAN (R 3.3.2)                  
##  methods     * 3.3.2      2016-10-31 local                           
##  munsell       0.4.3      2016-02-13 CRAN (R 3.2.3)                  
##  plyr          1.8.4      2016-06-08 CRAN (R 3.2.5)                  
##  Rcpp        * 0.12.11    2017-05-22 CRAN (R 3.3.2)                  
##  reshape2    * 1.4.2      2016-10-22 CRAN (R 3.3.0)                  
##  rlang         0.1.1      2017-05-18 CRAN (R 3.3.2)                  
##  rmarkdown     1.6        2017-06-15 CRAN (R 3.3.2)                  
##  rprojroot     1.2        2017-01-16 CRAN (R 3.3.2)                  
##  rstan         2.15.1     2017-04-19 CRAN (R 3.3.2)                  
##  scales        0.4.1      2016-11-09 CRAN (R 3.3.2)                  
##  StanHeaders   2.15.0-1   2017-04-19 CRAN (R 3.3.2)                  
##  stats       * 3.3.2      2016-10-31 local                           
##  stats4        3.3.2      2016-10-31 local                           
##  stringi       1.1.5      2017-04-07 cran (@1.1.5)                   
##  stringr       1.2.0      2017-02-18 cran (@1.2.0)                   
##  tibble        1.3.3      2017-05-28 CRAN (R 3.3.2)                  
##  tools         3.3.2      2016-10-31 local                           
##  utils       * 3.3.2      2016-10-31 local                           
##  withr         1.0.2      2016-06-20 CRAN (R 3.2.5)                  
##  yaml          2.1.14     2016-11-12 cran (@2.1.14)