Merge pull request #183 from ayushpatnaikgit/quantilese

Add SE for quantile
xKDR · Jan 19, 2023 · c7126a1 · c7126a1
2 parents e12c9e5 + b675dd9
commit c7126a1
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 35 deletions.
diff --git a/src/quantile.jl b/src/quantile.jl
@@ -13,32 +13,47 @@ The Julia, R and Python-numpy use the same defaults
 ```jldoctest
 julia> apisrs = load_data("apisrs");
 
-julia> srs = SurveyDesign(apisrs; weights=:pw);
+julia> srs = SurveyDesign(apisrs; weights=:pw) |> bootweights; 
 
 julia> quantile(:api00,srs,0.5)
 1×2 DataFrame
- Row │ probability  quantile 
-     │ Float64      Float64  
-─────┼───────────────────────
-   1 │         0.5     659.0
-
-julia> quantile(:enroll,srs,[0.1,0.2,0.5,0.75,0.95])
-5×2 DataFrame
- Row │ probability  quantile 
-     │ Float64      Float64  
-─────┼───────────────────────
-   1 │        0.1      245.5
-   2 │        0.2      317.6
-   3 │        0.5      453.0
-   4 │        0.75     668.5
-   5 │        0.95    1473.1
+ Row │ 0.5th percentile  SE      
+     │ Float64           Float64 
+─────┼───────────────────────────
+   1 │            659.0  14.9764
 ```
 """
-function quantile(var::Symbol, design::SurveyDesign, p::Union{<:Real,Vector{<:Real}}; 
-    alpha::Float64=0.05, ci::Bool=false, se::Bool=false, qrule="hf7",kwargs...)
+function quantile(var::Symbol, design::ReplicateDesign, p::Real;kwargs...)
     v = design.data[!, var]
     probs = design.data[!, design.allprobs]
-    df = DataFrame(probability=p, quantile=Statistics.quantile(v, ProbabilityWeights(probs), p))
-    # TODO: Add CI and SE of the quantile
+    X = Statistics.quantile(v, ProbabilityWeights(probs), p)
+    Xt = [Statistics.quantile(v, ProbabilityWeights(design.data[! , "replicate_"*string(i)]), p) for i in 1:design.replicates]
+    variance = sum((Xt .- X).^2) / design.replicates
+    df = DataFrame(percentile = X, SE = sqrt(variance))
+    rename!(df, :percentile => string(p) * "th percentile")
     return df
+end
+
+"""
+```jldoctest
+julia> apisrs = load_data("apisrs");
+
+julia> srs = SurveyDesign(apisrs; weights=:pw) |> bootweights; 
+
+julia> quantile(:enroll,srs,[0.1,0.2,0.5,0.75,0.95])
+5×3 DataFrame
+ Row │ percentile  statistic  SE       
+     │ String      Float64    Float64  
+─────┼─────────────────────────────────
+   1 │ 0.1             245.5   20.2964
+   2 │ 0.2             317.6   13.5435
+   3 │ 0.5             453.0   24.9719
+   4 │ 0.75            668.5   34.2487
+   5 │ 0.95           1473.1  142.568
+```
+"""
+function quantile(var::Symbol, design::ReplicateDesign, probs::Vector{<:Real}; kwargs...)
+    df = vcat([rename!(quantile(var, design, prob; kwargs...),[:statistic, :SE]) for prob in probs]...)
+    df.percentile = string.(probs)
+    return df[!, [:percentile, :statistic, :SE]]
 end
diff --git a/test/quantile.jl b/test/quantile.jl
@@ -6,21 +6,8 @@
     apisrs_original[!, :derived_sampsize] = fill(200.0, size(apisrs_original, 1))
     ##############################
     apisrs = copy(apisrs_original)
-    srs_design = SurveyDesign(apisrs; weights=:pw) 
-    @test quantile(:api00, srs_design, 0.5)[!,2][1] ≈ 659.0 atol=1e-4
+    srs_design = SurveyDesign(apisrs; weights=:pw) |> bootweights
+    @test quantile(:api00, srs_design, 0.5)[!,1][1] ≈ 659.0 atol=1e-4
     @test quantile(:api00, srs_design, [0.1753,0.25,0.5,0.75,0.975])[!,2] ≈ [512.8847,544,659,752.5,905] atol = 1e-4
     @test quantile(:enroll,srs_design, [0.1,0.2,0.5,0.75,0.95])[!,2] ≈ [245.5,317.6,453.0,668.5,1473.1] atol = 1e-4 
-end
-
-@testset "quantile_Stratified" begin
-    ##### StratifiedSample tests
-    # Load API datasets
-    apistrat_original = load_data("apistrat")
-    apistrat_original[!, :derived_probs] = 1 ./ apistrat_original.pw
-    apistrat_original[!, :derived_sampsize] = apistrat_original.fpc ./ apistrat_original.pw
-    # base functionality
-    apistrat = copy(apistrat_original)
-    dstrat = SurveyDesign(apistrat; strata = :stype, popsize = :fpc)
-    # Check which definition of quantile for StratifiedSample
-    # @test quantile(:enroll, dstrat, [0.1,0.2,0.5,0.75,0.95])[!,2] ≈ [262,309.3366,446.4103,658.8764,1589.7881] atol = 1e-4 
 end