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"
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 -------------------------------------------------------------
## 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)