pass weights to standardization and test against glmnet

Project.toml

@@ -1,6 +1,6 @@
 name = "Lasso"
 uuid = "b4fcebef-c861-5a0f-a7e2-ba9dc32b180a"
-version = "0.7.5"
+version = "0.7.6"
 
 [deps]
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"

README.md

@@ -78,9 +78,12 @@ More documentation is available at [![][docs-stable-img]][docs-stable-url].
 
 ## TODO
 
- - User-specified weights are untested
  - Maybe integrate LARS.jl
 
+## Changelog
+
+- v0.7.6 fixes a bug where user-specified weights were not used during standardization
+
 ## See also
 
  - [LassoPlot.jl](https://github.com/AsafManela/LassoPlot.jl), a package for

docs/src/index.md

@@ -114,7 +114,6 @@ tend to result in sparser coefficient estimates.
 
 ## TODO
 
- - User-specified weights are untested
  - Maybe integrate LARS.jl
 
 ## See also

src/Lasso.jl

@@ -357,10 +357,10 @@ initpenaltyfactor(penalty_factor::Nothing,p::Int,::Bool) = nothing
 initpenaltyfactor(penalty_factor::Vector,p::Int,standardizeω::Bool) =
     standardizeω ? rescale(penalty_factor, p) : penalty_factor
 
-# Standardize predictors if requested
-function standardizeX(X::AbstractMatrix{T}, standardize::Bool) where T
+# Standardize predictors if requested, accounting for observation weights `wts`
+function standardizeX(X::AbstractMatrix{T}, standardize::Bool, wts::AbstractVector=ones(T, size(X, 1))) where T
     if standardize
-        Xnorm = vec(convert(Matrix{T},std(X; dims=1, corrected=false)))
+        Xnorm = vec(convert(Matrix{T}, std(X, StatsBase.weights(wts), 1; corrected=false)))
         if any(x -> x == zero(T), Xnorm)
             @warn("""One of the predicators (columns of X) is a constant, so it can not be standardized.
                   To include a constant predicator set standardize = false and intercept = false""")
@@ -483,7 +483,7 @@ function StatsBase.fit(::Type{LassoPath},
     n = length(y)
     length(wts) == n || error("length(wts) = $(length(wts)) should be 0 or $n")
 
-    X, Xnorm = standardizeX(X, standardize)
+    X, Xnorm = standardizeX(X, standardize, wts)
 
     # Lasso initialization
     α = convert(T, α)

src/gammalasso.jl

@@ -75,7 +75,7 @@ function StatsBase.fit(::Type{GammaLassoPath},
     n = length(y)
     length(wts) == n || error("length(wts) = $(length(wts)) should be 0 or $n")
 
-    X, Xnorm = standardizeX(X, standardize)
+    X, Xnorm = standardizeX(X, standardize, wts)
 
     # Gamma lasso adaptation
     # Can potentially pass a different γ for each element of X, but if scalar we copy it to all params

test/gammalasso.jl

@@ -90,6 +90,20 @@ Random.seed!(testrng, 6540)
                     @test glp.λ == lp.λ
                     @test glp.b0 ≈ lp.b0
                     @test glp.coefs ≈ lp.coefs
+
+                    # Compare with LassoPath with weights and standardize=true (the
+                    # default), to exercise the same weighted standardizeX code path
+                    # that GammaLassoPath and LassoPath share
+                    Random.seed!(testrng, 5847)
+                    wts = rand(testrng, 10:1000, length(y)) ./ 100
+                    glp_wtd = fit(GammaLassoPath, X, y, dist, link; γ=γ, stopearly=false,
+                        λminratio=0.001, penalty_factor=penalty_factor, wts=wts,
+                        rng=StableRNG(1337))
+                    lp_wtd = fit(LassoPath, X, y, dist, link; stopearly=false,
+                        λminratio=0.001, penalty_factor=penalty_factor, λ=glp_wtd.λ,
+                        wts=wts, rng=StableRNG(1337))
+                    @test glp_wtd.b0 ≈ lp_wtd.b0
+                    @test glp_wtd.coefs ≈ lp_wtd.coefs
                 end
             end
         end

test/lasso.jl

@@ -157,6 +157,66 @@ end
     end
 end
 
+# Test against GLMNet with user-specified observation weights
+@testset "weights" begin
+    @testset "$(typeof(dist).name.name) $(typeof(link).name.name)" for (dist, link) in ((Normal(), IdentityLink()), (Binomial(), LogitLink()), (Poisson(), LogLink()))
+        Random.seed!(testrng, 3921)
+        (X, y) = genrand(Float64, dist, link, 1000, 10, false)
+        wts = rand(testrng, 10:1000, length(y)) ./ 100 # random positive weights, not summing to n
+
+        @testset "$(intercept ? "w/" : "w/o") intercept" for intercept = (false, true)
+            let y = y
+                if isa(dist, Normal)
+                    wymean = sum(wts .* y) / sum(wts)
+                    ypstd = sqrt(sum(wts .* abs2.(y .- wymean)) / sum(wts))
+                    yp = y ./ ypstd
+                    y = yp
+                    g = glmnet(X, yp, dist, intercept=intercept, tol=10*eps(); weights=wts)
+                elseif isa(dist, Binomial)
+                    yp = zeros(size(y, 1), 2)
+                    yp[:, 1] = y .== 0
+                    yp[:, 2] = y .== 1
+                    g = glmnet(X, yp, dist, intercept=intercept, tol=10*eps(); weights=wts)
+                else
+                    g = glmnet(X, y, dist, intercept=intercept, tol=10*eps(); weights=wts)
+                end
+                gbeta = convert(Matrix{Float64}, g.betas)
+
+                l = fit(LassoPath, X, y, dist, link, λ=g.lambda, intercept=intercept,
+                        cd_tol=10*eps(), irls_tol=10*eps(), wts=wts)
+
+                @test l.coefs ≈ gbeta rtol=1e-5
+                @test l.b0 ≈ g.a0 rtol=1e-5
+            end
+        end
+
+        # a constant rescaling of the weights should not change the fit, since
+        # wts is internally rescaled to sum to 1 (src/Lasso.jl)
+        @testset "scale invariance" begin
+            l1 = fit(LassoPath, X, y, dist, link, wts=wts, rng=StableRNG(1337))
+            l2 = fit(LassoPath, X, y, dist, link, wts=3.5 .* wts, rng=StableRNG(1337))
+            @test l1.coefs ≈ l2.coefs
+            @test l1.b0 ≈ l2.b0
+        end
+    end
+
+    # zero-weight observations should be equivalent to excluding them entirely
+    @testset "zero weights exclude observations" begin
+        Random.seed!(testrng, 3921)
+        (X, y) = genrand(Float64, Normal(), IdentityLink(), 200, 10, false)
+        keep = rand(testrng, Bool, length(y))
+        wts = Float64.(keep)
+
+        l_zeroed = fit(LassoPath, X, y, Normal(), IdentityLink(), wts=wts, rng=StableRNG(1337),
+                       cd_tol=1e-12, irls_tol=1e-12)
+        l_subset = fit(LassoPath, X[keep, :], y[keep], Normal(), IdentityLink(),
+                        λ=l_zeroed.λ, rng=StableRNG(1337), cd_tol=1e-12, irls_tol=1e-12)
+
+        @test l_zeroed.coefs ≈ l_subset.coefs rtol=1e-4
+        @test l_zeroed.b0 ≈ l_subset.b0 rtol=1e-4
+    end
+end
+
 # Test for sparse matrices
 
 # @testset "LassoPath Zero in" begin