Add HartleyRao for SurveyDesign

src/Survey.jl

@@ -25,6 +25,7 @@ include("boxplot.jl")
 include("show.jl")
 include("ratio.jl")
 include("by.jl")
+include("hr.jl")
 
 export load_data
 export AbstractSurveyDesign, SurveyDesign, ReplicateDesign
@@ -35,5 +36,6 @@ export hist, sturges, freedman_diaconis
 export boxplot
 export bootweights
 export ratio
+export HartleyRao
 
 end

src/hr.jl

@@ -0,0 +1,10 @@
+"""
+    Hartley Rao Approximation of the variance of the Horvitz-Thompson Estimator.
+    Avoids exhaustively calculating joint (inclusion) probabilities πᵢⱼ of the sampling scheme. 
+"""
+function HartleyRao(y::AbstractVector, design::SurveyDesign, HT_total)
+    πᵢ = design.data[!,design.allprobs]
+    hartley_rao_var = sum((1 .- ((design.data[!,design.sampsize] .- 1) ./ design.data[!,design.sampsize]) .* πᵢ) .*
+                          ((y .* design.data[!,design.weights]) .- (HT_total ./ design.data[!,design.sampsize])) .^ 2)
+    return hartley_rao_var
+end

src/total.jl

@@ -18,7 +18,8 @@ julia> total(:api00, dclus1)
 """
 function total(x::Symbol, design::SurveyDesign)
     X = wsum(design.data[!, x], weights(design.data[!, design.weights]))
-    DataFrame(total = X)
+    variance = HartleyRao(design.data[!, x], design, X)
+    DataFrame(total = X, SE = sqrt(variance))
 end
 
 """

test/hr.jl

@@ -0,0 +1,12 @@
+@testset "hartley_rao" begin
+    ### I havent been able to reproduce below in R survey, which give somewhat different SE estimates
+    # base functionality SRS
+    tot = total(:api00, srs)
+    @test tot.SE[1] ≈ 57154.1 rtol = STAT_TOL
+    # Stratified
+    tot = total(:api00, dstrat)
+    @test tot.SE[1] ≈ 699405 rtol = STAT_TOL
+    # one stage cluster
+    tot = total(:api00, dclus1)
+    @test tot.SE[1] ≈ 4880240 rtol = STAT_TOL
+end

test/runtests.jl

@@ -35,3 +35,4 @@ include("hist.jl")
 include("boxplot.jl")
 include("ratio.jl")
 include("show.jl")
+include("hr.jl")