Last updated: 2018-05-01

Code version: 4093362

source("../code/gdash_lik.R")
Loading required package: EQL
Loading required package: ttutils
Loading required package: SQUAREM
Loading required package: REBayes
Loading required package: Matrix
Loading required package: CVXR

Attaching package: 'CVXR'
The following object is masked from 'package:stats':

    power
Loading required package: PolynomF
Loading required package: Rmosek
Loading required package: ashr
source("../code/gdfit.R")
library(limma)
library(locfdr)
leukemia <- read.csv("http://web.stanford.edu/~hastie/CASI_files/DATA/leukemia_big.csv")
design <- c(rep(0, 20), rep(1, 14), rep(0, 27), rep(1, 11))
lim = limma::lmFit(leukemia, model.matrix(~design))
r.ebayes = limma::eBayes(lim)
p = r.ebayes$p.value[, 2]
t = r.ebayes$t[, 2]
z = -sign(t) * qnorm(p/2)
betahat = lim$coefficients[,2]
sebetahat = betahat / z
hist(p)

hist(z, prob = TRUE, ylim = c(0, dnorm(0)))
lines(seq(-10, 10, by = 0.01), dnorm(seq(-10, 10, by = 0.01)), col = "red")

fit <- gdash(betahat, sebetahat, gd.ord = 10)
rand.sign <- sample(c(rep(1, nrow(leukemia) / 2), rep(-1, nrow(leukemia) / 2)))
fit.sym <- gdash(rand.sign * betahat, sebetahat, gd.ord = 10)
fit.locfdr <- locfdr(z)

hist(z[abs(z) >= 5])

hist(betahat[abs(z) >= 5])

x.plot <- seq(-2, 2, length = 1000)
y.plot <- sapply(fit.sym$fitted_g$sd, function(x) {sapply(x.plot, dnorm, sd = x)})
y.plot <- y.plot %*% fit.sym$fitted_g$pi
plot(x.plot, y.plot, type = "l", col = "red")
y.plot <- sapply(fit$fitted_g$sd, function(x) {sapply(x.plot, dnorm, sd = x)})
y.plot <- y.plot %*% fit$fitted_g$pi
lines(x.plot, y.plot)

x.plot <- seq(-10, 10, length = 1000)
gd.ord <- 10
hermite = Hermite(gd.ord)
gd0.std = dnorm(x.plot)
matrix_lik_plot = cbind(gd0.std)
for (i in 1 : gd.ord) {
  gd.std = (-1)^i * hermite[[i]](x.plot) * gd0.std / sqrt(factorial(i))
  matrix_lik_plot = cbind(matrix_lik_plot, gd.std)
}
y.plot = matrix_lik_plot %*% fit$w * fit$fitted_g$pi[1]
hist(z, prob = TRUE, breaks = 50)
lines(x.plot, y.plot, col = "red")
lines(x.plot, dnorm(x.plot), col = "green", lty = 2)
lines(x.plot, dnorm(x.plot, .09, 1.68) * .93, col = 4, lty = 2)

library(lattice)
test <- seq(z)[fit$lfdr <= 0.1 & fit.locfdr$fdr > 0.1]
hist(z, prob = TRUE)

hist(z[test], xlim = c(-10, 10))

prostate <- read.csv("http://web.stanford.edu/~hastie/CASI_files/DATA/prostmat.csv")
design <- c(rep(0, 50), rep(1, 52))
lim = limma::lmFit(prostate, model.matrix(~design))
r.ebayes = limma::eBayes(lim)
p = r.ebayes$p.value[, 2]
t = r.ebayes$t[, 2]
z = -sign(t) * qnorm(p/2)
betahat = lim$coefficients[,2]
sebetahat = betahat / z
hist(p)

hist(z, prob = TRUE, ylim = c(0, dnorm(0)))
lines(seq(-10, 10, by = 0.01), dnorm(seq(-10, 10, by = 0.01)), col = "red")

fit <- gdash(betahat, sebetahat, gd.ord = 10)
rand.sign <- sample(c(rep(1, nrow(prostate) / 2), rep(-1, nrow(prostate) / 2)))
fit.sym <- gdash(rand.sign * betahat, sebetahat, gd.ord = 10)
Warning in rand.sign * betahat: longer object length is not a multiple of
shorter object length
fit.locfdr <- locfdr(z)

police <- scan("https://web.stanford.edu/~hastie/CASI_files/DATA/police.txt")
fit <- gdash(police, rep(1, length(police)), gd.ord = 10)
Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol

Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol

Warning in REBayes::KWDual(A, rep(1, k), normalize(w), control = control): estimated mixing distribution has some negative values:
               consider reducing rtol
DTI <- scan("https://web.stanford.edu/~hastie/CASI_files/DATA/DTI.txt")

Session information

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.4

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] lattice_0.20-35   locfdr_1.1-8      limma_3.34.4     
 [4] ashr_2.2-2        Rmosek_8.0.69     PolynomF_1.0-1   
 [7] CVXR_0.95         REBayes_1.2       Matrix_1.2-12    
[10] SQUAREM_2017.10-1 EQL_1.0-0         ttutils_1.0-1    

loaded via a namespace (and not attached):
 [1] gmp_0.5-13.1      Rcpp_0.12.16      compiler_3.4.3   
 [4] git2r_0.21.0      R.methodsS3_1.7.1 R.utils_2.6.0    
 [7] iterators_1.0.9   tools_3.4.3       digest_0.6.15    
[10] bit_1.1-12        evaluate_0.10.1   foreach_1.4.4    
[13] yaml_2.1.18       parallel_3.4.3    Rmpfr_0.6-1      
[16] ECOSolveR_0.4     stringr_1.3.0     knitr_1.20       
[19] rprojroot_1.3-2   bit64_0.9-7       grid_3.4.3       
[22] R6_2.2.2          rmarkdown_1.9     magrittr_1.5     
[25] splines_3.4.3     MASS_7.3-47       backports_1.1.2  
[28] codetools_0.2-15  htmltools_0.3.6   scs_1.1-1        
[31] stringi_1.1.6     doParallel_1.0.11 pscl_1.5.2       
[34] truncnorm_1.0-7   R.oo_1.21.0

This R Markdown site was created with workflowr