Fixed bugs identified by R CMD check and R CMD check --as-cran.

DESCRIPTION

@@ -2,13 +2,13 @@ Package: ashr
 Maintainer: Matthew Stephens <mstephens@uchicago.edu>
 Author: Matthew Stephens, Chaoxing Dai, Mengyin Lu, David Gerard, Nan Xiao, Peter Carbonetto
 Version: 2.1-6
-Date: 2017-01-19
+Date: 2017-03-15
 Title: Methods for Adaptive Shrinkage, using Empirical Bayes
 Description: The R package 'ashr' implements an Empirical Bayes approach for large-scale hypothesis testing and false discovery rate (FDR) estimation based on the methods proposed in M. Stephens, 2016, "False discovery rates: a new deal", <DOI:10.1093/biostatistics/kxw041>. These methods can be applied whenever two sets of summary statistics---estimated effects and standard errors---are available, just as 'qvalue' can be applied to previously computed p-values. Two main interfaces are provided: ash(), which is more user-friendly; and ash.workhorse(), which has more options and is geared toward advanced users. The ash() and ash.workhorse() also provides a flexible modeling interface that can accomodate a variety of likelihoods (e.g., normal, Poisson) and mixture priors (e.g., uniform, normal).
 Depends: R (>= 3.1.0)
-Imports: assertthat, truncnorm, SQUAREM, doParallel, pscl, Rcpp (>= 0.10.5), foreach, etrunct
+Imports: stats, assertthat, truncnorm, SQUAREM, doParallel, pscl, Rcpp (>= 0.10.5), foreach, etrunct
 LinkingTo: Rcpp
-Suggests: testthat, roxygen2, covr, knitr, rmarkdown
+Suggests: testthat, roxygen2, covr, knitr, rmarkdown, ggplot2
 Enhances: REBayes, Rmosek
 License: GPL (>=3)
 NeedsCompilation: no

R/igmix.R

@@ -35,7 +35,7 @@ comp_dens.igmix = function(m,y,log=FALSE){
   k=ncomp(m)
   n=length(y)
   d = matrix(rep(y,rep(k,n)),nrow=k)
-  return(matrix(dgamma(1/d, shape=m$alpha, rate=outer(m$beta,1/y^2),log),nrow=k))
+  return(matrix(stats::dgamma(1/d, shape=m$alpha, rate=outer(m$beta,1/y^2),log),nrow=k))
 }
 
 #density of product of each component of a inverse-gamma mixture with Gamma(v/2,v/2) at s

R/lik.R

@@ -76,17 +76,17 @@ lik_pois = function(y, scale=1, link=c("identity","log")){
   if (link=="identity"){
     list(name = "pois",
          const = TRUE,
-         lcdfFUN = function(x){pgamma(abs(x),shape=y+1,rate=scale,log.p=TRUE)+y*log(scale)},
-         lpdfFUN = function(x){dgamma(abs(x),shape=y+1,rate=scale,log=TRUE)+y*log(scale)},
+         lcdfFUN = function(x){stats::pgamma(abs(x),shape=y+1,rate=scale,log.p=TRUE)+y*log(scale)},
+         lpdfFUN = function(x){stats::dgamma(abs(x),shape=y+1,rate=scale,log=TRUE)+y*log(scale)},
          etruncFUN = function(a,b){-my_etruncgamma(-b,-a,y+1,scale)},
          e2truncFUN = function(a,b){my_e2truncgamma(-b,-a,y+1,scale)},
          data=list(y=y,link=link))
   }else if (link=="log"){
     y1 = y+1e-5 # add pseudocount
     list(name = "pois",
          const = TRUE,
-         lcdfFUN = function(x){pgamma(exp(-x),shape=y1,rate=scale,log.p=TRUE)-log(y1)+y*log(scale)},
-         lpdfFUN = function(x){dgamma(exp(-x),shape=y1,rate=scale,log=TRUE)-log(y1)+y*log(scale)},
+         lcdfFUN = function(x){stats::pgamma(exp(-x),shape=y1,rate=scale,log.p=TRUE)-log(y1)+y*log(scale)},
+         lpdfFUN = function(x){stats::dgamma(exp(-x),shape=y1,rate=scale,log=TRUE)-log(y1)+y*log(scale)},
          etruncFUN = function(a,b){-my_etruncgamma(exp(-b),exp(-a),y1,scale)},
          e2truncFUN = function(a,b){my_e2truncgamma(exp(-b),exp(-a),y1,scale)},
          data=list(y=y,link=link))
@@ -95,7 +95,7 @@ lik_pois = function(y, scale=1, link=c("identity","log")){
 
 #' @title Likelihood object for Binomial error distribution
 #' @description Creates a likelihood object for ash for use with Binomial error distribution
-#' @details Suppose we have Binomial observations \code{y} where \eqn{y_i\sim Bin(n_i,\p_i)}. 
+#' @details Suppose we have Binomial observations \code{y} where \eqn{y_i\sim Bin(n_i,p_i)}. 
 #'    We either put an unimodal prior g on the success probabilities \eqn{p_i\sim g} (by specifying 
 #'    \code{link="identity"}) or on the logit success probabilities \eqn{logit(p_i)\sim g} 
 #'    (by specifying \code{link="logit"}). Either way, ASH with this Binomial likelihood function 
@@ -116,8 +116,8 @@ lik_binom = function(y,n,link=c("identity","logit")){
   if (link=="identity"){
     list(name = "binom",
          const = TRUE,
-         lcdfFUN = function(x){pbeta(abs(x),shape1=y+1,shape2=n-y+1,log.p=TRUE)-log(n+1)},
-         lpdfFUN = function(x){dbeta(abs(x),shape1=y+1,shape2=n-y+1,log=TRUE)-log(n+1)},
+         lcdfFUN = function(x){stats::pbeta(abs(x),shape1=y+1,shape2=n-y+1,log.p=TRUE)-log(n+1)},
+         lpdfFUN = function(x){stats::dbeta(abs(x),shape1=y+1,shape2=n-y+1,log=TRUE)-log(n+1)},
          etruncFUN = function(a,b){-my_etruncbeta(-b,-a,y+1,n-y+1)},
          e2truncFUN = function(a,b){my_e2truncbeta(-b,-a,y+1,n-y+1)},
          data=list(y=y,n=n,link=link))
@@ -126,8 +126,8 @@ lik_binom = function(y,n,link=c("identity","logit")){
     n1 = n+2e-5
     list(name = "binom",
          const = TRUE,
-         lcdfFUN = function(x){pbeta(1/(1+exp(-x)),shape1=y1,shape2=n1-y1,log.p=TRUE)+log(n1/(y1*(n1-y1)))},
-         lpdfFUN = function(x){dbeta(1/(1+exp(-x)),shape1=y1,shape2=n1-y1,log=TRUE)+log(n1/(y1*(n1-y1)))},
+         lcdfFUN = function(x){stats::pbeta(1/(1+exp(-x)),shape1=y1,shape2=n1-y1,log.p=TRUE)+log(n1/(y1*(n1-y1)))},
+         lpdfFUN = function(x){stats::dbeta(1/(1+exp(-x)),shape1=y1,shape2=n1-y1,log=TRUE)+log(n1/(y1*(n1-y1)))},
          etruncFUN = function(a,b){-my_etruncbeta(-1/(1+exp(-b)),-1/(1+exp(-a)),y1,n1-y1)},
          e2truncFUN = function(a,b){my_e2truncbeta(-1/(1+exp(-b)),-1/(1+exp(-a)),y1,n1-y1)},
          data=list(y=y,n=n,link=link))

R/truncbeta.R

@@ -6,12 +6,12 @@
 #' @export
 my_etruncbeta = function(a, b, alpha, beta){
   tmp = a
-  tmp[a!=b] = (alpha/(alpha+beta)*(pbeta(b,alpha+1,beta)-pbeta(a,alpha+1,beta))/
-                 (pbeta(b,alpha,beta)-pbeta(a,alpha,beta)))[a!=b]
-  # zero denominator case: pbeta(b,alpha,beta) and pbeta(a,alpha,beta) are both 0 or 1 
-  tmp[(pbeta(b,alpha,beta)-pbeta(a,alpha,beta))==0] = 
-    ifelse(dbeta(a,alpha,beta,log=TRUE)>dbeta(b,alpha,beta,log=TRUE), 
-           a, b)[(pbeta(b,alpha,beta)-pbeta(a,alpha,beta))==0]
+  tmp[a!=b] = (alpha/(alpha+beta)*(stats::pbeta(b,alpha+1,beta)-stats::pbeta(a,alpha+1,beta))/
+                 (stats::pbeta(b,alpha,beta)-stats::pbeta(a,alpha,beta)))[a!=b]
+  # zero denominator case: stats::pbeta(b,alpha,beta) and stats::pbeta(a,alpha,beta) are both 0 or 1 
+  tmp[(stats::pbeta(b,alpha,beta)-stats::pbeta(a,alpha,beta))==0] = 
+    ifelse(stats::dbeta(a,alpha,beta,log=TRUE)>stats::dbeta(b,alpha,beta,log=TRUE), 
+           a, b)[(stats::pbeta(b,alpha,beta)-stats::pbeta(a,alpha,beta))==0]
   return(tmp)
 }
 
@@ -24,11 +24,11 @@ my_etruncbeta = function(a, b, alpha, beta){
 my_e2truncbeta = function(a, b, alpha, beta){
   tmp = a^2
   tmp[a!=b] = (alpha*(alpha+1)/((alpha+beta)*(alpha+beta+1))*
-                 (pbeta(b,alpha+2,beta)-pbeta(a,alpha+2,beta))/
-                 (pbeta(b,alpha,beta)-pbeta(a,alpha,beta)))[a!=b]
-  # zero denominator case: pbeta(b,alpha,beta) and pbeta(a,alpha,beta) are both 0 or 1 
-  tmp[(pbeta(b,alpha,beta)-pbeta(a,alpha,beta))==0] = 
-    ifelse(dbeta(a,alpha,beta,log=TRUE)>dbeta(b,alpha,beta,log=TRUE), 
-           a^2, b^2)[(pbeta(b,alpha,beta)-pbeta(a,alpha,beta))==0]
+                 (stats::pbeta(b,alpha+2,beta)-stats::pbeta(a,alpha+2,beta))/
+                 (stats::pbeta(b,alpha,beta)-stats::pbeta(a,alpha,beta)))[a!=b]
+  # zero denominator case: stats::pbeta(b,alpha,beta) and stats::pbeta(a,alpha,beta) are both 0 or 1 
+  tmp[(stats::pbeta(b,alpha,beta)-stats::pbeta(a,alpha,beta))==0] = 
+    ifelse(stats::dbeta(a,alpha,beta,log=TRUE)>stats::dbeta(b,alpha,beta,log=TRUE), 
+           a^2, b^2)[(stats::pbeta(b,alpha,beta)-stats::pbeta(a,alpha,beta))==0]
   return(tmp)
-}
+}

R/truncgamma.R

@@ -7,12 +7,12 @@
 #' @export
 my_etruncgamma = function(a, b, shape, rate){
   tmp = a
-  tmp[a!=b] = (shape/rate*(pgamma(b,shape=shape+1,rate=rate)-pgamma(a,shape=shape+1,rate=rate))/
-                 (pgamma(b,shape=shape,rate=rate)-pgamma(a,shape=shape,rate=rate)))[a!=b]
-  # zero denominator case: pgamma(b,shape,rate) and pgamma(a,shape,rate) are both 0 or 1 
-  tmp[(pgamma(b,shape=shape,rate=rate)-pgamma(a,shape=shape,rate=rate))==0] = 
-    ifelse(dgamma(a,shape=shape,rate=rate,log=TRUE)>dgamma(b,shape=shape,rate=rate,log=TRUE), 
-           a, b)[(pgamma(b,shape=shape,rate=rate)-pgamma(a,shape=shape,rate=rate))==0]
+  tmp[a!=b] = (shape/rate*(stats::pgamma(b,shape=shape+1,rate=rate)-stats::pgamma(a,shape=shape+1,rate=rate))/
+                 (stats::pgamma(b,shape=shape,rate=rate)-stats::pgamma(a,shape=shape,rate=rate)))[a!=b]
+  # zero denominator case: stats::pgamma(b,shape,rate) and stats::pgamma(a,shape,rate) are both 0 or 1 
+  tmp[(stats::pgamma(b,shape=shape,rate=rate)-stats::pgamma(a,shape=shape,rate=rate))==0] = 
+    ifelse(stats::dgamma(a,shape=shape,rate=rate,log=TRUE)>stats::dgamma(b,shape=shape,rate=rate,log=TRUE), 
+           a, b)[(stats::pgamma(b,shape=shape,rate=rate)-stats::pgamma(a,shape=shape,rate=rate))==0]
   return(tmp)
 }
 
@@ -25,11 +25,11 @@ my_etruncgamma = function(a, b, shape, rate){
 #' @export
 my_e2truncgamma = function(a, b, shape, rate){
   tmp = a^2
-  tmp[a!=b] = (shape*(shape+1)/rate^2*(pgamma(b,shape=shape+2,rate=rate)-pgamma(a,shape=shape+2,rate=rate))/
-                 (pgamma(b,shape=shape,rate=rate)-pgamma(a,shape=shape,rate=rate)))[a!=b]
-  # zero denominator case: pgamma(b,shape,rate) and pgamma(a,shape,rate) are both 0 or 1 
-  tmp[(pgamma(b,shape=shape,rate=rate)-pgamma(a,shape=shape,rate=rate))==0] = 
-    ifelse(dgamma(a,shape=shape,rate=rate,log=TRUE)>dgamma(b,shape=shape,rate=rate,log=TRUE), 
-           a^2, b^2)[(pgamma(b,shape=shape,rate=rate)-pgamma(a,shape=shape,rate=rate))==0]
+  tmp[a!=b] = (shape*(shape+1)/rate^2*(stats::pgamma(b,shape=shape+2,rate=rate)-stats::pgamma(a,shape=shape+2,rate=rate))/
+                 (stats::pgamma(b,shape=shape,rate=rate)-stats::pgamma(a,shape=shape,rate=rate)))[a!=b]
+  # zero denominator case: stats::pgamma(b,shape,rate) and stats::pgamma(a,shape,rate) are both 0 or 1 
+  tmp[(stats::pgamma(b,shape=shape,rate=rate)-stats::pgamma(a,shape=shape,rate=rate))==0] = 
+    ifelse(stats::dgamma(a,shape=shape,rate=rate,log=TRUE)>stats::dgamma(b,shape=shape,rate=rate,log=TRUE), 
+           a^2, b^2)[(stats::pgamma(b,shape=shape,rate=rate)-stats::pgamma(a,shape=shape,rate=rate))==0]
   return(tmp)
-}
+}

tests/testthat/test_lik.R

@@ -2,7 +2,7 @@ test_that("general likelihood with multiple df works", {
   df = c(rep(100,50),rep(2,50))
   s = rgamma(100,1,1)
   betahat = s*rt(n=100,df=df)
-  data =set_data(betahat,s,lik = lik_normal),alpha=0)
+  data = set_data(betahat,s,lik = lik_normal(),alpha=0)
   expect_equal(is_normal(data$lik),TRUE)
   expect_equal(is_const(data$lik),TRUE)
 

tests/testthat/test_optmethod.R

@@ -1,14 +1,15 @@
-test_that("optmethod switches to EM when given negative probabilities", {
-	set.seed(777)
-        betahat <- c(0.7775, 0.6091, 0.6880, 0.6897, 0.6565, 0.7505,
-                     0.7125, 0.7201, 0.7498, 0.7553)
-        sebetahat <- rep(0.01, length = length(betahat))
-        ash_out <- ash(betahat = betahat, sebetahat = sebetahat,
-                             mixcompdist = "uniform", outputlevel=4)
-        expect_true(min(ash_out$fitted_g$pi) > -10 ^ -12) 
-        expect_true(ash_out$fit_details$optmethod == "mixEM" | ash_out$fit_details$optmethod != "mixIP")
-#note: had to add check for mixIP to make this run on travis
-})
+# This test is no longer relevant because we have modified this behaviour.
+# test_that("optmethod switches to EM when given negative probabilities", {
+# 	set.seed(777)
+#         betahat <- c(0.7775, 0.6091, 0.6880, 0.6897, 0.6565, 0.7505,
+#                      0.7125, 0.7201, 0.7498, 0.7553)
+#         sebetahat <- rep(0.01, length = length(betahat))
+#         ash_out <- ash(betahat = betahat, sebetahat = sebetahat,
+#                              mixcompdist = "uniform", outputlevel=4)
+#         expect_true(min(ash_out$fitted_g$pi) > -10 ^ -12) 
+#         expect_true(ash_out$fit_details$optmethod == "mixEM" | ash_out$fit_details$optmethod != "mixIP")
+# #note: had to add check for mixIP to make this run on travis
+# })
 
 test_that("control is passed to optmethod correctly when method is mixIP", {
   testthat::skip_if_not_installed("REBayes")