diff --git a/Project.toml b/Project.toml index a0ea707..52b0386 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index 54cc6e8..72ff354 100644 --- a/README.md +++ b/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 diff --git a/docs/src/index.md b/docs/src/index.md index 99af39e..cb73821 100644 --- a/docs/src/index.md +++ b/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 diff --git a/src/Lasso.jl b/src/Lasso.jl index 2297431..3f1b149 100644 --- a/src/Lasso.jl +++ b/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, α) diff --git a/src/gammalasso.jl b/src/gammalasso.jl index 465b98b..a005118 100644 --- a/src/gammalasso.jl +++ b/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 diff --git a/test/gammalasso.jl b/test/gammalasso.jl index c94c712..1b1c9ac 100644 --- a/test/gammalasso.jl +++ b/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 diff --git a/test/lasso.jl b/test/lasso.jl index 51742f6..dbfac43 100644 --- a/test/lasso.jl +++ b/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