WIP: Add GLM with design based standard errors.

Project.toml

@@ -9,6 +9,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Missings = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

src/Survey.jl

@@ -12,6 +12,8 @@ using AlgebraOfGraphics
 using CategoricalArrays
 using Random
 using Missings
+using GLM
+import GLM: glm
 
 include("SurveyDesign.jl")
 include("bootstrap.jl")
@@ -26,6 +28,7 @@ include("show.jl")
 include("ratio.jl")
 include("by.jl")
 include("jackknife.jl")
+include("reg.jl")
 
 export load_data
 export AbstractSurveyDesign, SurveyDesign, ReplicateDesign
@@ -37,5 +40,6 @@ export boxplot
 export bootweights
 export ratio
 export jackknifeweights, jackknife_variance
+export @formula, glm
 
 end

src/reg.jl

@@ -0,0 +1,15 @@
+"""
+```julia
+apisrs = load_data("apisrs")
+srs = SurveyDesign(apisrs)
+bsrs = bootweights(srs, replicates = 2000)
+glm(@formula(api00~api99), bsrs, Normal())
+```
+"""
+function glm(formula::FormulaTerm, design::ReplicateDesign, args...; kwargs...)
+    X = coef(glm(formula, design.data, args...; wts = design.data[!, design.weights], kwargs...))
+    Xt = [coef(glm(formula, design.data, args...; wts = design.data[!, "replicate_"*string(i)]), kwargs...) for i in 1:design.replicates]
+    Xt = hcat(Xt...)
+    SE = [std(Xt[i,:]) for i in 1:size(Xt)[1]]
+    DataFrame(Coefficients = X, SE = SE)
+end

test/Project.toml

@@ -3,3 +3,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"

test/reg.jl

@@ -0,0 +1,6 @@
+using GLM
+@testset "GLM on simple random sample" begin
+    repdesign = bootweights(srs; replicates = 500)
+   @test glm(@formula(api00~api99), repdesign, Normal()).Coefficients ≈ coef(glm(@formula(api00~api99), repdesign.data, Normal()))
+   @test glm(@formula(api00~api99), repdesign, Normal()).SE ≈ stderror(glm(@formula(api00~api99), repdesign.data, Normal())) rtol = 1
+end

test/runtests.jl

@@ -63,4 +63,5 @@ include("hist.jl")
 include("boxplot.jl")
 include("ratio.jl")
 include("show.jl")
-include("jackknife.jl")
+include("jackknife.jl")
+include("reg.jl")