diff --git a/.travis.yml b/.travis.yml index a63c101..6b3594a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,7 +4,7 @@ os: - osx julia: - 1.0 - - 1.1 + - 1.3 - nightly notifications: email: false @@ -17,8 +17,7 @@ before_install: - if [ $TRAVIS_OS_NAME = linux ]; then sudo apt-get install gfortran -y; else - brew cask uninstall oclint; - brew install gcc; + brew update && brew install gcc; fi after_success: @@ -27,7 +26,7 @@ after_success: jobs: include: - stage: "Documentation" - julia: 1.1 + julia: 1.3 os: linux script: - julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' diff --git a/Project.toml b/Project.toml index 7750b76..f4a3bfd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Lasso" uuid = "b4fcebef-c861-5a0f-a7e2-ba9dc32b180a" -version = "0.5.0" +version = "0.5.1" [deps] DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" @@ -15,7 +15,12 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" [compat] +DSP = "0.6" +Distributions = "0.21" GLM = "1.3" +MLBase = "0.8" +Reexport = "0.2" +StatsBase = "0.32" StatsModels = "0.6" julia = "0.7, 1" diff --git a/docs/src/index.md b/docs/src/index.md index 237ec08..15f8e55 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -41,30 +41,28 @@ julia> data = DataFrame(X=[1,2,3], Y=[2,4,7]) │ 3 │ 3 │ 7 │ julia> m = fit(LassoModel, @formula(Y ~ X), data) -StatsModels.TableRegressionModel{LassoModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredQR{Float64}}},Array{Float64,2}} - -Y ~ X +LassoModel using MinAICc(2) segment of the regulatization path. Coefficients: -────────────────────────────────────────────────────────────────── - Estimate Std. Error t value Pr(>|t|) Lower 95% Upper 95% -────────────────────────────────────────────────────────────────── -x1 3.88915 4.86043 0.800166 0.5704 -57.8684 65.6467 -x2 0.222093 1.83707 0.120895 0.9234 -23.1201 23.5643 -────────────────────────────────────────────────────────────────── - -julia> stderror(m) -2-element Array{Float64,1}: - 4.860427341926979 - 1.8370688588910695 +──────────── + Estimate +──────────── +x1 3.88915 +x2 0.222093 +──────────── julia> predict(m) 3-element Array{Float64,1}: - 4.161154511001072 - 4.433161908091373 - 4.705169305181673 - + 4.111240223622052 + 4.333333333333333 + 4.555426443044614 + +julia> predict(m, data[2:end,:]) +2-element Array{Union{Missing, Float64},1}: + 4.333333333333333 + 4.555426443044614 ``` + To get an entire Lasso regularization path with default parameters: ```julia diff --git a/docs/src/lasso.md b/docs/src/lasso.md index ddab76e..74e20de 100644 --- a/docs/src/lasso.md +++ b/docs/src/lasso.md @@ -73,17 +73,15 @@ julia> data = DataFrame(X=[1,2,3], Y=[2,4,7]) │ 3 │ 3 │ 7 │ julia> m = fit(LassoModel, @formula(Y ~ X), data; select=MinCVmse(Kfold(3,2))) -StatsModels.TableRegressionModel{LassoModel{LinearModel{GLM.LmResp{Array{Float64,1}},GLM.DensePredQR{Float64}}},Array{Float64,2}} - -Y ~ X +LassoModel using MinCVmse(Kfold([3, 1, 2], 2, 1.5)) segment of the regulatization path. Coefficients: -────────────────────────────────────────────────────────────────── - Estimate Std. Error t value Pr(>|t|) Lower 95% Upper 95% -────────────────────────────────────────────────────────────────── -x1 4.33333 5.4365 0.797081 0.5716 -64.744 73.4106 -x2 0.0 2.0548 0.0 1.0000 -26.1088 26.1088 -────────────────────────────────────────────────────────────────── +──────────── + Estimate +──────────── +x1 4.33333 +x2 0.0 +──────────── julia> coef(m) 2-element Array{Float64,1}: diff --git a/src/Lasso.jl b/src/Lasso.jl index 1157bf6..be554e9 100644 --- a/src/Lasso.jl +++ b/src/Lasso.jl @@ -14,7 +14,7 @@ include("TrendFiltering.jl") using Reexport, LinearAlgebra, SparseArrays, Random, .Util, MLBase import Random: Sampler @reexport using GLM, Distributions, .FusedLassoMod, .TrendFiltering -using GLM: FPVector, LinPred, Link, LmResp, GlmResp, DensePredQR, updateμ!, +using GLM: FPVector, LinPred, Link, LmResp, GlmResp, cholpred, updateμ!, linkfun, linkinv, LinPredModel export RegularizationPath, LassoPath, GammaLassoPath, NaiveCoordinateDescent, CovarianceCoordinateDescent, fit, fit!, coef, predict, @@ -597,6 +597,20 @@ function StatsBase.predict(path::RegularizationPath, newX::AbstractMatrix{T}; of end end +"predicted values for data used to estimate `path`" +function StatsBase.predict(path::RegularizationPath; select=AllSeg()) + X = path.m.pp.X + offset = path.m.rr.offset + + # destandardize X if needed + if !isempty(path.Xnorm) + X = X ./ transpose(path.Xnorm) + end + + predict(path, X; offset=offset, select=select) +end + + "distribution of underlying GLM" distfun(path::RegularizationPath{M}) where {M<:LinearModel} = Normal() distfun(path::RegularizationPath) = path.m.rr.d diff --git a/src/segselect.jl b/src/segselect.jl index 1f63a8a..2d13534 100644 --- a/src/segselect.jl +++ b/src/segselect.jl @@ -119,15 +119,17 @@ segselect(path::RegularizationPath, abstract type RegularizedModel <: RegressionModel end "LassoModel represents a selected segment from a LassoPath" -struct LassoModel{M<:LinPredModel} <: RegularizedModel +struct LassoModel{M<:LinPredModel, S<:SegSelect} <: RegularizedModel lpm::M # underlying GLM intercept::Bool # whether path added an intercept + select::S # segment selector end "GammaLassoModel represents a selected segment from a GammaLassoPath" -struct GammaLassoModel{M<:LinPredModel} <: RegularizedModel +struct GammaLassoModel{M<:LinPredModel, S<:SegSelect} <: RegularizedModel lpm::M # underlying GLM intercept::Bool # whether path added an intercept + select::S # segment selector end """ @@ -140,7 +142,7 @@ Predicted values using a selected segment of a regularization path. m = fit(LassoModel, X, y; select=MinBIC()) predict(m, newX) # predict using BIC minimizing segment """ -function StatsBase.predict(m::RegularizedModel, newX::AbstractMatrix{T}; kwargs...) where {T<:AbstractFloat} +function StatsBase.predict(m::RegularizedModel, newX::AbstractMatrix{T}; kwargs...) where T # add an interecept to newX if the model has one if m.intercept newX = [ones(T,size(newX,1),1) newX] @@ -169,8 +171,14 @@ selectmodel(path, MinCVmse(path, 5)) # 5-fold CV mse minimizing model function selectmodel(path::R, select::SegSelect) where R<:RegularizationPath # extract reusable path parts m = path.m + rr = deepcopy(m.rr) pp = m.pp - X = pp.X + X = pp.X + + # destandardize X if needed + if !isempty(path.Xnorm) + X = X ./ transpose(path.Xnorm) + end # add an interecept to X if the model has one if hasintercept(path) @@ -183,17 +191,21 @@ function selectmodel(path::R, select::SegSelect) where R<:RegularizationPath beta0 = Vector{Float64}(coef(path, select)) # create new linear predictor - segpp = DensePredQR(segX, beta0) + pivot = true + p = cholpred(segX, pivot) # rescale weights, which in GLM sum to nobs - m.rr.wts .*= nobs(path) + rr.wts .*= nobs(path) # same things GLM does to init just before fit! - GLM.installbeta!(segpp) - updateμ!(m.rr, linpred(segpp, zero(eltype(m.rr.y)))) + lp = rr.mu + copyto!(p.beta0, beta0) + fill!(p.delbeta, 0) + GLM.linpred!(lp, p, 0) + updateμ!(rr, lp) # create a LinearModel or GeneralizedLinearModel with the new linear predictor - newglm(m, segpp) + newglm(m, rr, p) end """ @@ -275,11 +287,11 @@ function StatsBase.fit(::Type{R}, X::AbstractMatrix{T}, y::V, M = pathtype(R) path = fit(M, X, y, d, l; intercept=intercept, kwargs...) - R(selectmodel(path, select), intercept) + R(selectmodel(path, select), intercept, select) end -newglm(m::LinearModel, pp) = LinearModel(m.rr, pp) -newglm(m::GeneralizedLinearModel, pp) = GeneralizedLinearModel(m.rr, pp, true) +newglm(m::LinearModel, rr, pp) = LinearModel(rr, pp) +newglm(m::GeneralizedLinearModel, rr, pp) = GeneralizedLinearModel(rr, pp, true) # don't add an intercept when using a @formula because we use the intercept keyword arg to add an intercept StatsModels.drop_intercept(::Type{R}) where R<:RegularizedModel = true @@ -291,9 +303,7 @@ for modeltype in (:LassoModel, :GammaLassoModel) StatsBase.deviance, StatsBase.nulldeviance, StatsBase.loglikelihood, StatsBase.nullloglikelihood, StatsBase.dof, StatsBase.dof_residual, StatsBase.nobs, - StatsBase.stderror, StatsBase.vcov, - StatsBase.residuals, StatsBase.response, - StatsBase.coeftable + StatsBase.residuals, StatsBase.response ] end end @@ -306,3 +316,19 @@ StatsBase.fit(::Type{M}, d::UnivariateDistribution=Normal(), l::Link=canonicallink(d); kwargs...) where {M<:Union{RegularizationPath, RegularizedModel}} = fit(M, float(X), float(y), d, l; kwargs...) + +function coeftable(mm::RegularizedModel) + cc = coef(mm) + CoefTable([cc], + ["Estimate"], + ["x$i" for i = 1:size(mm.lpm.pp.X, 2)]) +end + +function Base.show(io::IO, obj::RegularizedModel) + # prefix = isa(obj.m, GeneralizedLinearModel) ? string(typeof(distfun(path)).name.name, " ") : "" + println(io, "$(typeof(obj).name.name) using $(obj.select) segment of the regulatization path.") + + println(io, "\nCoefficients:\n", coeftable(obj)) +end + +StatsBase.vcov(obj::RegularizedModel, args...) = error("variance-covariance matrix for a regularized model is not yet implemented") diff --git a/test/segselect.jl b/test/segselect.jl index 2177662..1cc4eae 100644 --- a/test/segselect.jl +++ b/test/segselect.jl @@ -1,5 +1,18 @@ # tests for segment selection and StatsModels @formula+df interface using Random +using GLM: predict + +"Tests whether the output of `show(expr)` converted to String contains `expected`" +macro test_show(expr, expected) + s = quote + local o = IOBuffer() + show(o, $expr) + String(take!(o)) + end + substr = quote $expected end + :(@test occursin($(esc(substr)), $(esc(s)))) +end + @testset "segment/model selection" begin datapath = joinpath(dirname(@__FILE__), "data") @@ -25,13 +38,18 @@ datapath = joinpath(dirname(@__FILE__), "data") Random.seed!(421) pathpredict = Lasso.predict(path, data; select=select, offset=offset) + @test pathpredict ≈ Lasso.predict(path; select=select) @test pathcoefs == coef(m) if isa(dist, Normal) - @test pathpredict == GLM.predict(m, data) + offset + @test pathpredict ≈ predict(m, data) + offset + @test pathpredict ≈ predict(m) else - @test pathpredict == GLM.predict(m, data; offset=offset) + @test pathpredict ≈ predict(m, data; offset=offset) + @test pathpredict ≈ predict(m) end + + @test_show m string(L) end end end