Penalized generalized linear models
glm_lasso(X, y, λ, loss = ProximalBase.Logistic(), options = CDOptions();
intercept = true, standardize_x = true, weights, offset)
GLMLassoPath(X, y, loss = ProximalBase.Logistic(); nλ = 100, λminratio, kwargs...)minimize ∑ᵢ wᵢ ℓ(β₀ + xᵢ'β + oᵢ, yᵢ) / ∑ᵢ wᵢ + λ ∑ⱼ ωⱼ |βⱼ| for a scalar loss ℓ of ProximalBase, with glmnet's conventions and defaults: an unpenalized intercept, loadings ωⱼ equal to the standard deviations of the columns, observation weights w and an offset o. They run proximal Newton (glmnet's method): each step minimizes the weighted least-squares model of the loss by coordinate descent, then backtracks.
using HighDimensionalModels, ProximalBase
using LinearAlgebra, Printf, Random, Statistics
Random.seed!(10)
logistic(z) = 1 / (1 + exp(-z))
n, p = 400, 60
X = randn(n, p)
β = zeros(p)
β[1:5] = [1.2, -1.0, 0.8, -0.6, 0.5]
y = Float64.(rand(n) .< logistic.(-0.5 .+ X * β))
support(x) = findall(!iszero, x)support (generic function with 1 method)Logistic regression
sol = glm_lasso(X, y, 0.03)
@printf "intercept %.3f, selected %s\n" sol.intercept support(sol.x)
@printf "explained deviance %.3f, proximal Newton took %d coordinate descent passes\n" (1 - sol.deviance / sol.nulldeviance) sol.iterationsintercept -0.477, selected [1, 2, 3, 4, 5, 8, 15, 18, 31, 47, 53, 54]
explained deviance 0.338, proximal Newton took 48 coordinate descent passesthe optimality conditions: |∇ⱼ| ≤ λωⱼ at the zeros, ∇ⱼ = -λωⱼ sign(βⱼ) elsewhere, and the gradient over the intercept is zero
μ = logistic.(sol.linearpredictor)
G = X' * (μ .- y) / n
ω = vec(std(X; dims = 1, corrected = false))
viol = maximum(j -> sol.x[j] == 0 ? max(abs(G[j]) - 0.03ω[j], 0) : abs(G[j] + 0.03ω[j] * sign(sol.x[j])), 1:p)
@printf "largest violation of the optimality conditions %.1e, mean residual %.1e\n" viol mean(μ .- y)largest violation of the optimality conditions 8.2e-09, mean residual 2.7e-15The path, and λ chosen on a validation set
train, valid = 1:300, 301:400
path = GLMLassoPath(X[train, :], y[train]; nλ = 50)GLMLassoPath{Float64, ProximalBase.Logistic}([0.19510695512592358, 0.1616739682987798, 0.13396996538952125, 0.11101324359961898, 0.09199032199998082, 0.07622711549785939, 0.06316504835286048, 0.052341260814614536, 0.04337220749613243, 0.0359400662844239 … 0.00010591723353333637, 8.7767550088161e-5, 7.272794607171619e-5, 6.0265486896892846e-5, 4.993856016967893e-5, 4.138122697137324e-5, 3.4290254661687766e-5, 2.8414371704754242e-5, 2.354536375835124e-5, 1.9510695512592358e-5], SparseArrays.SparseVector{Float64, Int64}[sparsevec(Int64[], Float64[], 60), sparsevec([1, 2], [0.1325199768797053, -0.08843277278210537], 60), sparsevec([1, 2], [0.24381309467004642, -0.1994359862550953], 60), sparsevec([1, 2, 4], [0.3457516504940985, -0.305246704801011, -0.039079263098870276], 60), sparsevec([1, 2, 4], [0.4325264252273989, -0.40205301893802387, -0.11973651377586769], 60), sparsevec([1, 2, 4, 31], [0.5120464327088226, -0.49085104717751005, -0.19112981261625417, 0.0042157300472036585], 60), sparsevec([1, 2, 3, 4, 5, 8, 31], [0.5845039587879217, -0.5668286244950901, 0.017645054814523276, -0.2511811728415721, 0.030832692333852052, -0.04500730304524203, 0.06415576670373546], 60), sparsevec([1, 2, 3, 4, 5, 8, 15, 31, 54], [0.6551464669496736, -0.6343699378144864, 0.06077559496688412, -0.30212975514030005, 0.08383363149783146, -0.09676649871395915, 0.006029880591114478, 0.11805641114528709, -0.014435809439577063], 60), sparsevec([1, 2, 3, 4, 5, 8, 15, 31, 53, 54], [0.7205727462089757, -0.6978924840661913, 0.0987972596828949, -0.3482509919946637, 0.12742789457960454, -0.1413730523904955, 0.029219924789951307, 0.16701352674780445, -0.019221657334421133, -0.05294169592273937], 60), sparsevec([1, 2, 3, 4, 5, 8, 11, 15, 18, 31, 35, 53, 54], [0.7805137961843701, -0.759005125994126, 0.13420341732866667, -0.3925467424343456, 0.16549134934665666, -0.18295357961522074, 0.009897946482867542, 0.04700064028788045, -0.0061604086076087194, 0.21197141995382282, -0.02183709452555388, -0.05480849159723236, -0.08492578598426709], 60) … sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [1.9944269931909981, -1.9289210495952576, 0.6711769954576101, -1.2151730112326429, 0.9075422022860564, -0.24395452767769407, -0.15452070258432316, -0.8299461593522135, 0.0515818967106458, 0.09205145204541539 … -0.3500954756454031, 0.12320009828914676, -0.5505865601800081, -0.4833551511093641, 0.18405032094250273, -0.018876847436186357, 0.2530088068215691, -0.055925713341574354, -0.03737264830913149, 0.14216881225992437], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [1.998206431839425, -1.9325046490957742, 0.6728450982240384, -1.218054414451755, 0.9099780315781067, -0.2450189122185276, -0.15510208178428664, -0.8317991470504197, 0.05113188897133278, 0.09305395245908953 … -0.3514190437671138, 0.12383702337496064, -0.5513781325273299, -0.48432349521305096, 0.18468078840117738, -0.019324142527973508, 0.25347055718733125, -0.05634530027134869, -0.037543562050051445, 0.1431998175807086], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.0014088394629166, -1.9355402031518203, 0.6742496030542848, -1.220490533739135, 0.9120340702300614, -0.2459130826904064, -0.1555913989925556, -0.8333692269366827, 0.05075264761371635, 0.09389306893679338 … -0.3525357570241784, 0.12437250028788911, -0.5520483406419325, -0.48513816867246523, 0.18520652058986037, -0.01969936543985795, 0.2538574550791583, -0.056698886195391554, -0.03768023328635354, 0.14406131086068308], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.0040650683066707, -1.9380555880052377, 0.6754150360624266, -1.2225094618524182, 0.9137416854921954, -0.24665885557303563, -0.1559933621930162, -0.834669320298979, 0.050435481766815456, 0.09458938109552636 … -0.3534612579234538, 0.12481622784234596, -0.5526018359511066, -0.4858169664994135, 0.18564907170729353, -0.020009424717065332, 0.2541778355285959, -0.05699099125583438, -0.03779552100622794, 0.14477699551116013], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.0062490643605244, -1.9401307548699795, 0.6763763977230332, -1.2241734785786333, 0.9151458892468823, -0.24726852679963826, -0.15632427835139845, -0.8357427171514551, 0.05017887883679841, 0.09516379073220951 … -0.35422318659322255, 0.1251829689915382, -0.5530606985028915, -0.4863735130425728, 0.18601219029973576, -0.02026723143162494, 0.25444427938332553, -0.05723144102767988, -0.03789314241154608, 0.14536821254745982], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.0080433059338048, -1.9418349462534288, 0.6771686799356276, -1.225540022135576, 0.9163029822295741, -0.24777392649213126, -0.156601297482797, -0.8366222127709241, 0.04996258543556182, 0.09563890377127302 … -0.3548506048352407, 0.12548899096508986, -0.5534343779316085, -0.4868318707144224, 0.18631109977885663, -0.020481492138928355, 0.2546621663910045, -0.05743054968364382, -0.03797754489216085, 0.14585734935353353], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.0095868928694, -1.9432958750168854, 0.6778445274956214, -1.2267139626748123, 0.9172938037811833, -0.24820369001783196, -0.15683193270654372, -0.8373812784845883, 0.04978009939769791, 0.09604203894296244 … -0.3553872013979839, 0.12574397571997867, -0.5537572187398768, -0.4872254656155466, 0.1865645768600237, -0.020661124555916264, 0.2548488039235045, -0.057598998940505394, -0.03804151258916042, 0.14627059923931549], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.0108438118673764, -1.944487723598133, 0.678396988200026, -1.2276702666941484, 0.918101676983922, -0.24855541821541166, -0.15702265071094998, -0.8379982927126208, 0.049630197707968074, 0.09637208192461189 … -0.35582521507003917, 0.12595483996127801, -0.5540197848317764, -0.4875461995590479, 0.18677199658777457, -0.02080918909205174, 0.25500110406556176, -0.057737140709735645, -0.038097047483187596, 0.14660942505585992], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.011922246039422, -1.9455096897926534, 0.6788675654284079, -1.2284870827988008, 0.9187911663576328, -0.24885418312446922, -0.15718397169520396, -0.838527572128609, 0.049501191288805625, 0.09664969886644854 … -0.35619797700361255, 0.12613374898875748, -0.5542428639870969, -0.48781786987140846, 0.1869488589370988, -0.020933808061983708, 0.2551284061338004, -0.057855698950267265, -0.03814071590944035, 0.14689651910470725], 60), sparsevec([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 51, 52, 53, 54, 55, 56, 57, 58, 59, 60], [2.012753995546527, -1.9462982498948789, 0.6792361333568824, -1.2291232822039562, 0.9193307143715845, -0.24909009090901846, -0.15731165059606664, -0.8389360244113162, 0.049401325500377624, 0.09687261626163603 … -0.356489931143561, 0.12627577047435076, -0.55441795761592, -0.48803342491819035, 0.1870887799576153, -0.021034634331296057, 0.2552320822704048, -0.05794658565548844, -0.038182365673826414, 0.14712541823661768], 60)], [-0.295464212893836, -0.29765342002583134, -0.301662946706352, -0.3141898210597486, -0.3350889272528429, -0.355569932408567, -0.3770707074787938, -0.40175236079292154, -0.42529330010058236, -0.4465871464464334 … -1.1015297356143168, -1.1040291559960926, -1.106141406959424, -1.1078938751991794, -1.109337562089062, -1.110523232819727, -1.1115425133387824, -1.1123724112330582, -1.1130813264645578, -1.1136340797524704], [0.0, 0.05663232259591677, 0.10442321327860304, 0.14833922587823434, 0.1878594902424764, 0.21801042080784838, 0.2541375677957989, 0.28845811336105753, 0.31823689097839747, 0.3438261090762975 … 0.5473493736509671, 0.5473581033614737, 0.5473642011914648, 0.5473683877785255, 0.547371247079, 0.547373202024408, 0.5473745855708438, 0.5473755214267683, 0.5473761804919723, 0.5473766089766634], [0, 5, 5, 8, 9, 10, 12, 11, 11, 12 … 14, 16, 13, 14, 13, 12, 17, 19, 11, 15], Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], true, ProximalBase.Logistic(), 409.41163720647046, 300, ProximalBase.NormL1{Float64, Vector{Float64}}(1.0, [0.9970152780295275, 0.983812924378902, 1.0462316366928042, 1.044938071891675, 1.0066195799053697, 0.9948728499518679, 1.0460126454524794, 0.906438853186948, 0.8873074292119646, 0.9830905313365047 … 1.0089050695180182, 0.9807931956835988, 0.9355297188408502, 0.9659143051869733, 1.0066644305539305, 0.98192419723847, 1.011557965041003, 1.0021427377009162, 1.036736645746164, 0.9684855112374441]), nothing)the validation deviance of each solution
function deviance(η, y)
2 * sum(log1p.(exp.(η)) .- y .* η)
end
dev = [deviance(path.β0path[k] .+ X[valid, :] * path.βpath[k], y[valid]) for k in eachindex(path.λpath)]
k = argmin(dev)
@printf "%d levels, explaining up to %.1f%% of the deviance (a path stops at 99.9%%)\n" length(path.λpath) 100 * last(path.devratio)
@printf "best λ = %.4f (level %d): %d nonzero coefficients, validation deviance %.1f\n" path.λpath[k] k count(!iszero, path.βpath[k]) dev[k]50 levels, explaining up to 54.7% of the deviance (a path stops at 99.9%)
best λ = 0.0298 (level 11): 15 nonzero coefficients, validation deviance 97.7Poisson regression with an exposure offset
Poisson draws by Knuth's method, to avoid a dependency on Distributions
function rpois(μ)
k, t, L = 0, rand(), exp(-μ)
while t > L
k += 1
t *= rand()
end
k
end
exposure = rand(n) .* 4 .+ 1
counts = Float64.(rpois.(exposure .* exp.(0.2 .+ 0.4 .* X[:, 1] .- 0.3 .* X[:, 2])))400-element Vector{Float64}:
4.0
2.0
3.0
5.0
6.0
3.0
2.0
2.0
1.0
5.0
⋮
1.0
7.0
3.0
5.0
0.0
3.0
0.0
4.0
0.0λmax, the first level of a path, is the scale of the penalty for this loss
λmaxp = GLMLassoPath(X, counts, ProximalBase.Poisson(); nλ = 1, offset = log.(exposure)).λpath[1]
solp = glm_lasso(X, counts, 0.2 * λmaxp, ProximalBase.Poisson(); offset = log.(exposure))
@printf "λmax %.3f; at 0.2 λmax: intercept %.3f, selected %s\n" λmaxp solp.intercept support(solp.x)λmax 1.621; at 0.2 λmax: intercept 0.234, selected [1, 2]observation weights: integer weights are replicated rows
w = Float64.(rand(1:3, n))
solw = glm_lasso(X, y, 0.03; weights = w)
rows = reduce(vcat, [fill(i, Int(w[i])) for i in 1:n])
solr = glm_lasso(X[rows, :], y[rows], 0.03)
@printf "weighted fit and replicated rows agree to %.1e\n" norm(Vector(solw.x) - Vector(solr.x), Inf)weighted fit and replicated rows agree to 2.2e-08