Merge pull request #18 from juliangehring/two-step-adjust-method

Add adjustment method selection for TwoStep π0 procedure

Author: juliangehring

README.md

@@ -54,15 +54,18 @@ adjust(pvals, ForwardStop())
 * Censored BUM
 * BUM
 * FlatGrenander
+* Oracle for known π0
 
 ```julia
 estimate_pi0(pvals, Storey())
 estimate_pi0(pvals, StoreyBootstrap())
 estimate_pi0(pvals, LeastSlope())
 estimate_pi0(pvals, TwoStep())
 estimate_pi0(pvals, TwoStep(0.05))
+estimate_pi0(pvals, TwoStep(0.05, BenjaminiHochbergAdaptive(0.9))
 estimate_pi0(pvals, RightBoundary())
 estimate_pi0(pvals, CensoredBUM())
 estimate_pi0(pvals, BUM())
 estimate_pi0(pvals, FlatGrenander())
+estimate_pi0(pvals, Oracle(0.9))
 ```

src/pi0-estimators.jl

@@ -160,18 +160,20 @@ Reference: Benjamini, Krieger and Yekutieli, 2006
 """
 immutable TwoStep <: Pi0Estimator
     α::AbstractFloat
-    ## method::PValueAdjustmentMethod
+    method::PValueAdjustmentMethod
 
-    TwoStep(α) = isin(α, 0., 1.) ? new(α) : throw(DomainError())
+    TwoStep(α, method) = isin(α, 0., 1.) ? new(α, method) : throw(DomainError())
 end
 
 TwoStep() = TwoStep(0.05)
 
+TwoStep(α) = TwoStep(α, BenjaminiHochberg())
+
 function estimate_pi0{T<:AbstractFloat}(pValues::Vector{T}, pi0estimator::TwoStep)
-    twostep_pi0(pValues, pi0estimator.α)
+    twostep_pi0(pValues, pi0estimator.α, pi0estimator.method)
 end
 
-function twostep_pi0{T<:AbstractFloat}(pValues::Vector{T}, alpha::T, method::PValueAdjustmentMethod = BenjaminiHochberg())
+function twostep_pi0{T<:AbstractFloat}(pValues::Vector{T}, alpha::T, method::PValueAdjustmentMethod)
     padj = adjust(pValues, method)
     pi0 = sum(padj .>= (alpha/(1+alpha))) / length(padj)
     return(pi0)

test/test-pi0-estimators.jl

@@ -127,24 +127,27 @@ using StatsBase
         alpha = 0.05
 
         ## checked against mutoss::TSBKY_pi0_est
-        @test twostep_pi0(p, alpha) ≈ 0.665
-        @test twostep_pi0(p0, alpha) ≈ 1.0
-        @test twostep_pi0(p1, alpha) ≈ 0.29
+        @test twostep_pi0(p, alpha, BenjaminiHochberg()) ≈ 0.665
+        @test twostep_pi0(p0, alpha, BenjaminiHochberg()) ≈ 1.0
+        @test twostep_pi0(p1, alpha, BenjaminiHochberg()) ≈ 0.29
 
         @test issubtype(typeof(TwoStep(alpha)), Pi0Estimator)
         @test estimate_pi0(p, TwoStep()) ≈ 0.665
+        @test estimate_pi0(p, TwoStep()) ≈ 0.665
         @test estimate_pi0(p, TwoStep(alpha)) ≈ 0.665
         @test estimate_pi0(p0, TwoStep(alpha)) ≈ 1.0
         @test estimate_pi0(p1, TwoStep(alpha)) ≈ 0.29
+        @test estimate_pi0(p, TwoStep(alpha, BenjaminiHochberg())) ≈ 0.665
 
         @test estimate_pi0(p, TwoStep(0.1)) ≈ 0.63
         @test estimate_pi0(p, TwoStep(0.0)) ≈ 1.0
         @test estimate_pi0(p, TwoStep(1.0)) ≈ 0.415
+        @test estimate_pi0(p, TwoStep(0.1, BenjaminiHochberg())) ≈ 0.63
 
         ## unsorted p-values
         p_unsort = unsort(p)
         @test !issorted(p_unsort)
-        @test twostep_pi0(p_unsort, alpha) ≈ 0.665
+        @test estimate_pi0(p_unsort, TwoStep(alpha)) ≈ 0.665
         @test !issorted(p_unsort)
 
     end