Exhaustively testing SimpleRandomSample, software eng improvements

src/SurveyDesign.jl

@@ -15,36 +15,102 @@ abstract type AbstractSurveyDesign end
     SimpleRandomSample <: AbstractSurveyDesign
 
     Survey design sampled by simple random sampling.
+    # Required arguments:
+    data - This is the survey dataset loaded as a DataFrame in memory. 
+           Note: Keeping with Julia conventions, original data object
+           is modified, not copied. Be careful
+    # Optional arguments:
+    sampsize -  Sample size of the survey, given as Symbol name of 
+                column in `data`, an `Unsigned` integer, or a Vector
+    popsize  -  The (expected) population size of survey, given as Symbol 
+                name of column in `data`, an `Unsigned` integer, or a Vector
+    weights  -  Sampling weights, passed as Symbol or Vector
+    probs    -  Sampling probabilities, passed as Symbol or Vector
+    ignorefpc-  Ignore finite population correction and assume all weights equal to 1.0
+
+    Precedence order of using `popsize`, `weights` and `probs` is `popsize` > `weights` > `probs` 
+    Eg. if `popsize` given then assumed ground truth over `weights` or `probs`
+
+    If `popsize` not given, `weights` or `probs` must be given, so that in combination 
+    with `sampsize`, `popsize` can be calculated.
 """
 struct SimpleRandomSample <: AbstractSurveyDesign
     data::AbstractDataFrame
-    sampsize::Union{Nothing,Unsigned}
-    popsize::Union{Nothing,Unsigned}
-    sampfraction::Real
-    fpc::Real
+    sampsize::Union{Unsigned,Nothing}
+    popsize::Union{Unsigned,Nothing}
+    sampfraction::Float64
+    fpc::Float64
     ignorefpc::Bool
     function SimpleRandomSample(data::AbstractDataFrame;
         popsize=nothing,
-        sampsize=nrow(data),
+        sampsize=nrow(data) |> UInt,
         weights=nothing,
         probs=nothing,
         ignorefpc=false
     )
+        # Only valid argument types given to constructor
+        argtypes_weights = Union{Nothing,Symbol,Vector{<:Real}}
+        argtypes_probs = Union{Nothing,Symbol,Vector{<:Real}}
+        argtypes_popsize = Union{Nothing,Symbol,<:Unsigned,Vector{<:Real}}
+        argtypes_sampsize = Union{Nothing,Symbol,<:Unsigned,Vector{<:Real}}
+        # If any invalid type raise error
+        if !(isa(weights, argtypes_weights))
+            error("Invalid type of argument given for `weights` argument")
+        elseif !(isa(probs, argtypes_probs))
+            error("Invalid type of argument given for `probs` argument")
+        elseif !(isa(popsize, argtypes_popsize))
+            error("Invalid type of argument given for `popsize` argument")
+        elseif !(isa(sampsize, argtypes_sampsize))
+            error("Invalid type of argument given for `sampsize` argument")
+        end
+        # If any of weights or probs given as Symbol,
+        # find the corresponding column in `data`
         if isa(weights, Symbol)
             weights = data[!, weights]
         end
         if isa(probs, Symbol)
             probs = data[!, probs]
         end
-
-        if ignorefpc
-            @warn "assuming all weights are equal to 1.0"
-            weights = ones(nrow(data))
+        # If weights/probs vector not numeric/real, ie. string column passed for weights, then raise error
+        if !isa(weights, Union{Nothing,Vector{<:Real}})
+            error("Weights should be Vector{<:Real}. You passed $(typeof(weights))")
+        elseif !isa(probs, Union{Nothing,Vector{<:Real}})
+            error("Sampling probabilities should be Vector{<:Real}. You passed $(typeof(probs))")
         end
-
-        # set population size if it is not given; `weights` and `sampsize` must be given
+        # If popsize given as Symbol or Vector, check all records equal 
+        if isa(popsize, Symbol)
+            if !all(w -> w == first(data[!, popsize]), data[!, popsize])
+                error("popsize must be same for all observations in Simple Random Sample")
+            end
+            popsize = first(data[!, popsize]) |> UInt
+        elseif isa(popsize, Vector{<:Real})
+            if !all(w -> w == first(popsize), popsize)
+                error("popsize must be same for all observations in Simple Random Sample")
+            end
+            popsize = first(popsize) |> UInt
+        end
+        # If sampsize given as Symbol or Vector, check all records equal 
+        if isa(sampsize, Symbol)
+            if !all(w -> w == first(data[!, sampsize]), data[!, sampsize])
+                error("sampsize must be same for all observations in Simple Random Sample")
+            end
+            sampsize = first(data[!, sampsize]) |> UInt
+        elseif isa(sampsize, Vector{<:Real})
+            if !all(w -> w == first(sampsize), sampsize)
+                error("sampsize must be same for all observations in Simple Random Sample")
+            end
+            sampsize = first(sampsize) |> UInt
+        end
+        # If both `weights` and `probs` given, then `weights` is assumed to be ground truth for probs.
+        if !isnothing(weights) && !isnothing(probs)
+            probs = 1 ./ weights
+            data[!, :probs] = probs
+        end
+        # popsize must be nothing or <:Integer by now
         if isnothing(popsize)
-            # check that all weights are equal (SRS is by definition equi-weighted)
+            # If popsize not given, fallback to weights, probs and sampsize to estimate `popsize`
+            @warn "Using weights/probs and sampsize to estimate `popsize`"
+            # Check that all weights (or probs if weights not given) are equal, as SRS is by definition equi-weighted
             if typeof(weights) <: Vector{<:Real}
                 if !all(w -> w == first(weights), weights)
                     error("all frequency weights must be equal for Simple Random Sample")
@@ -53,23 +119,46 @@ struct SimpleRandomSample <: AbstractSurveyDesign
                 if !all(p -> p == first(probs), probs)
                     error("all probability weights must be equal for Simple Random Sample")
                 end
-                weights = 1 / probs
+                weights = 1 ./ probs
+            else
+                error("either weights or probs must be given if `popsize` not given")
             end
-            # estimate population size
+            # Estimate population size
             popsize = round(sampsize * first(weights)) |> UInt
-            # @show popsize, sampsize # Check this Line TODO
             if sampsize > popsize
-                error("sample size cannot be greater than population size")
+                error("population size was estimated to be greater than given sampsize. Please check input arguments.")
             end
-        elseif typeof(popsize) <: Vector{<:Real}
-            if !all(y -> y == first(popsize), popsize) # SRS by definition is equi-weighted
-                error("Simple Random Sample must be equi-weighted. Different sampling weights detected in vectors")
-            end
-            weights = popsize ./ sampsize # ratio estimator for SRS
-            popsize = first(popsize) |> UInt
+        elseif typeof(popsize) <: Unsigned
+            weights = fill(popsize / sampsize, nrow(data)) # If popsize is given, weights vector is made concordant with popsize and sampsize, regardless of given weights argument
         else
-            error("either population size or frequency/probability weights must be specified")
+            error("Something went wrong. Please check validity of inputs.")
+        end
+        # If ignorefpc then set weights to 1 ??
+        # TODO: This works under some cases, but should find better way to process ignoring fpc
+        if ignorefpc
+            @warn "assuming all weights are equal to 1.0"
+            weights = ones(nrow(data))
+            probs = 1 ./ weights
         end
+        # sum of weights must equal to `popsize` for SRS
+        if !isnothing(weights) && !(isapprox(sum(weights), popsize; atol=1e-4))
+            if ignorefpc && !(isapprox(sum(weights), sampsize; atol=1e-4)) # Change if ignorefpc functionality changes
+                error("Sum of sampling weights should be equal to `sampsize` for Simple Random Sample with ignorefpc")
+            elseif !ignorefpc
+                @show sum(weights)
+                error("Sum of sampling weights must be equal to `popsize` for Simple Random Sample")
+            end
+        end
+        # sum of probs must equal popsize for SRS
+        if !isnothing(probs) && !(isapprox(sum(1 ./ probs), popsize; atol=1e-4))
+            if ignorefpc && !(isapprox(sum(1 ./ probs), sampsize; atol=1e-4)) # Change if ignorefpc functionality changes
+                error("Sum of inverse sampling probabilities should be equal to `sampsize` for Simple Random Sample with ignorefpc")
+            elseif !ignorefpc
+                @show sum(1 ./ probs)
+                error("Sum of inverse of sampling probabilities must be equal to `popsize` for Simple Random Sample")
+            end
+        end
+        ## Set remaining parts of data structure
         # set sampling fraction
         sampfraction = sampsize / popsize
         # set fpc
@@ -80,7 +169,7 @@ struct SimpleRandomSample <: AbstractSurveyDesign
             probs = 1 ./ data[!, :weights]
         end
         data[!, :probs] = probs
-
+        # Initialise the structure
         new(data, sampsize, popsize, sampfraction, fpc, ignorefpc)
     end
 end
@@ -105,6 +194,9 @@ struct StratifiedSample <: AbstractSurveyDesign
         probs=nothing,
         ignorefpc=false
     )
+        if isa(popsize, Symbol)
+            popsize = data[!, popsize]
+        end
         if isa(weights, Symbol)
             weights = data[!, weights]
         end
@@ -141,7 +233,7 @@ struct StratifiedSample <: AbstractSurveyDesign
         # set sampling fraction
         sampfraction = sampsize ./ popsize
         # set fpc
-        fpc = ignorefpc ? fill(1,size(data, 1)) : 1 .- (sampsize ./ popsize)
+        fpc = ignorefpc ? fill(1, size(data, 1)) : 1 .- (sampsize ./ popsize)
         # add columns for frequency and probability weights to `data`
         if !isnothing(probs)
             data[!, :probs] = probs
@@ -282,8 +374,8 @@ struct SurveyDesign <: AbstractSurveyDesign
         if ignorefpc # && (isnothing(popsize) || isnothing(weights) || isnothing(probs))
             @warn "Assuming equal weights"
             weights = ones(nrow(data))
-        end    
-        
+        end
+
         # TODO: Do the other case where clusters are given
         if isnothing(clusters)
             # set population size if it is not given; `weights` and `sampsize` must be given
@@ -327,7 +419,7 @@ struct SurveyDesign <: AbstractSurveyDesign
             end
             # @show clusters, strata, sampsize,popsize, sampfraction, fpc, ignorefpc
             new(data, clusters, strata, sampsize, popsize, sampfraction, fpc, ignorefpc)
-        elseif isa(clusters,Symbol)
+        elseif isa(clusters, Symbol)
             # One Cluster sampling - PSU chosen with SRS,
             print("One stage cluster design with PSU SRS")
         elseif typeof(clusters) <: Vector{Symbol}

src/svymean.jl

@@ -54,7 +54,7 @@ function svymean(x::Symbol, design::SimpleRandomSample)
         p.se = sqrt.(p.var)
         return p
     end
-    return DataFrame(mean = mean(design.data[!, x]), sem = sem(x, design::SimpleRandomSample))
+    return DataFrame(mean=mean(design.data[!, x]), sem=sem(x, design::SimpleRandomSample))
 end
 
 function svymean(x::Vector{Symbol}, design::SimpleRandomSample)
@@ -84,52 +84,48 @@ function sem_svyby(x::AbstractVector, design::SimpleRandomSample)
     return sqrt(variance)
 end
 
-function svymean(x::AbstractVector, design::SimpleRandomSample, weights)
-    return DataFrame(mean = mean(x), sem = sem_svyby(x, design))
-end
-
 """
 Inner method for `svyby` for SimpleRandomSample
 """
 function svymean(x::AbstractVector, design::SimpleRandomSample, weights)
-    return DataFrame(mean = mean(x), sem = sem_svyby(x, design))
+    return DataFrame(mean=mean(x), sem=sem_svyby(x, design))
 end
 
 """
 Inner method for `svyby` for StratifiedSample
 Calculates domain mean and its std error, based example 10.3.3 on pg394 Sarndal (1992)
 """
-function svymean(x::AbstractVector, popsize::AbstractVector,sampsize::AbstractVector,sampfraction::AbstractVector,strata::AbstractVector)
-    df = DataFrame(x = x, popsize = popsize, sampsize = sampsize, sampfraction = sampfraction,strata = strata)
+function svymean(x::AbstractVector, popsize::AbstractVector, sampsize::AbstractVector, sampfraction::AbstractVector, strata::AbstractVector)
+    df = DataFrame(x=x, popsize=popsize, sampsize=sampsize, sampfraction=sampfraction, strata=strata)
     nsdh = []
     substrata_domain_totals = []
     Nh = []
     nh = []
     fh = []
     ȳsdh = []
     sigma_ȳsh_squares = []
-    grouped_frame = groupby(df,:strata)
+    grouped_frame = groupby(df, :strata)
     for each_strata in keys(grouped_frame)
         nsh = nrow(grouped_frame[each_strata])#, nrow=>:nsdh).nsdh
-        push!(nsdh,nsh)
+        push!(nsdh, nsh)
         substrata_domain_total = sum(grouped_frame[each_strata].x)
         ȳdh = mean(grouped_frame[each_strata].x)
         push!(ȳsdh, ȳdh)
-        push!(substrata_domain_totals,substrata_domain_total)
+        push!(substrata_domain_totals, substrata_domain_total)
         popsizes = first(grouped_frame[each_strata].popsize)
-        push!(Nh,popsizes)
+        push!(Nh, popsizes)
         sampsizes = first(grouped_frame[each_strata].sampsize)
-        push!(nh,sampsizes)
+        push!(nh, sampsizes)
         sampfrac = first(grouped_frame[each_strata].sampfraction)
-        push!(fh,sampfrac)
-        push!(sigma_ȳsh_squares,sum((grouped_frame[each_strata].x .- ȳdh).^2) )
+        push!(fh, sampfrac)
+        push!(sigma_ȳsh_squares, sum((grouped_frame[each_strata].x .- ȳdh) .^ 2))
     end
-    domain_mean = sum(Nh .* substrata_domain_totals ./ nh)/sum(Nh .* nsdh ./ nh)
-    pdh = nsdh./nh
+    domain_mean = sum(Nh .* substrata_domain_totals ./ nh) / sum(Nh .* nsdh ./ nh)
+    pdh = nsdh ./ nh
     N̂d = sum(Nh .* pdh)
-    domain_var = sum(Nh.^2 .* (1 .- fh) .*  (sigma_ȳsh_squares .+ (nsdh .* (1 .-pdh) .* (ȳsdh .- domain_mean).^2)) ./ (nh .* (nh .- 1) )) ./ N̂d.^2
+    domain_var = sum(Nh .^ 2 .* (1 .- fh) .* (sigma_ȳsh_squares .+ (nsdh .* (1 .- pdh) .* (ȳsdh .- domain_mean) .^ 2)) ./ (nh .* (nh .- 1))) ./ N̂d .^ 2
     domain_mean_se = sqrt(domain_var)
-    return DataFrame(domain_mean = domain_mean, domain_mean_se = domain_mean_se)
+    return DataFrame(domain_mean=domain_mean, domain_mean_se=domain_mean_se)
 end
 
 """
@@ -162,12 +158,12 @@ function svymean(x::Symbol, design::StratifiedSample)
     s²ₕ = combine(gdf, x => var => :s²h).s²h
     V̂Ȳ̂ = sum((Wₕ .^ 2) .* (1 .- fₕ) .* s²ₕ ./ nₕ)
     SE = sqrt(V̂Ȳ̂)
-    return DataFrame(Ȳ̂ = Ȳ̂, SE = SE)
+    return DataFrame(Ȳ̂=Ȳ̂, SE=SE)
 end
 
 function svymean(::Bool; x::Symbol, design::StratifiedSample)
     gdf = groupby(design.data, design.strata)
     ȳₕ = combine(gdf, x => mean => :mean).mean
     s²ₕ = combine(gdf, x => var => :s²h).s²h
-    return DataFrame(ȳₕ,s²ₕ)
+    return DataFrame(ȳₕ, s²ₕ)
 end

src/svytotal.jl

@@ -43,7 +43,7 @@ function svytotal(x::Symbol, design::SimpleRandomSample)
         return p
     end
     total = design.popsize * mean(design.data[!, x])
-    return DataFrame(total = total, se_total = se_tot(x, design::SimpleRandomSample))
+    return DataFrame(total=total, se_total=se_tot(x, design::SimpleRandomSample))
 end
 
 # Inner methods for `svyby`
@@ -57,7 +57,7 @@ function se_total_svyby(x::AbstractVector, design::SimpleRandomSample, _)
 end
 function svytotal(x::AbstractVector, design::SimpleRandomSample, weights)
     total = wsum(x, weights)
-    return DataFrame(total = total, sem = se_total_svyby(x, design::SimpleRandomSample, weights))
+    return DataFrame(total=total, sem=se_total_svyby(x, design::SimpleRandomSample, weights))
 end
 
 # StratifiedSample
@@ -82,7 +82,7 @@ function svytotal(x::Symbol, design::StratifiedSample)
     # the only difference between total and mean variance is the Nₕ instead of Wₕ
     V̂Ȳ̂ = sum((Nₕ .^ 2) .* (1 .- fₕ) .* s²ₕ ./ nₕ)
     SE = sqrt(V̂Ȳ̂)
-    return DataFrame(grand_total = grand_total, SE = SE)
+    return DataFrame(grand_total=grand_total, SE=SE)
 end
 
 function svytotal(x::Vector{Symbol}, design::SimpleRandomSample)