Skip to content

Commit 4167f0c

Browse files
BenjaminDoranBenjaminDoran
authored
Add methods for calling multi-response gaussian family (#66)
Fixes #65
1 parent 47da7a0 commit 4167f0c

19 files changed

+1111
-10
lines changed

‎Project.toml

+4-2
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "GLMNet"
22
uuid = "8d5ece8b-de18-5317-b113-243142960cc6"
3-
version = "0.7.2"
3+
version = "0.7.3"
44

55
[deps]
66
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -22,6 +22,8 @@ julia = "1.6"
2222
[extras]
2323
CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597"
2424
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
25+
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
26+
2527

2628
[targets]
27-
test = ["Test", "CategoricalArrays"]
29+
test = ["Test", "CategoricalArrays", "CSV"]

‎README.md

+2-1
Original file line numberDiff line numberDiff line change
@@ -145,7 +145,7 @@ julia> vline!([lambdamin(iris_cv)])
145145

146146
## Fitting models
147147

148-
`glmnet` has two required parameters: the n x m predictor matrix `X` and the dependent variable `y`. It additionally accepts an optional third argument, `family`, which can be used to specify a generalized linear model. Currently, `Normal()` (least squares, default), `Binomial()` (logistic), `Poisson()` , `Multinomial()`, `CoxPH()` (Cox model) are supported.
148+
`glmnet` has two required parameters: the n x m predictor matrix `X` and the dependent variable `y`. It additionally accepts an optional third argument, `family`, which can be used to specify a generalized linear model. Currently, `Normal()` (least squares, default), `MvNormal()` Multi-response gaussian, `Binomial()` (logistic), `Poisson()` , `Multinomial()`, `CoxPH()` (Cox model) are supported.
149149

150150
- For linear and Poisson models, `y` is a numerical vector of length n.
151151
- For logistic models, `y` is either a string vector of length n or a n x 2 matrix, where the first column is the count of negative responses for each row in `X` and the second column is the count of positive responses.
@@ -166,6 +166,7 @@ julia> vline!([lambdamin(iris_cv)])
166166
- `lambda`: The λ values to consider. By default, this is determined from `nlambda` and `lambda_min_ratio`.
167167
- `tol`: Convergence criterion, with the default value of `1e-7`.
168168
- `standardize`: Whether to standardize predictors so that they are in the same units, with the default value of `true`. Beta values are always presented on the original scale.
169+
- `standardize_response`: (only for `MvNormal()`), Whether to standardize the response variables so that they are the same units. defaults to `false`.
169170
- `intercept`: Whether to fit an intercept term. The intercept is always unpenalized. The default value is `true`.
170171
- `maxit`: The maximum number of iterations of the cyclic coordinate descent algorithm. If convergence is not achieved, a warning is returned. The default value is `1e6`.
171172

‎src/GLMNet.jl

+5-4
Original file line numberDiff line numberDiff line change
@@ -235,7 +235,7 @@ macro validate_and_init()
235235
if !isempty(lambda)
236236
# user-specified lambda values
237237
nlambda == 100 || error("cannot specify both lambda and nlambda")
238-
lambda_min_ratio == (length(y) < size(X, 2) ? 1e-2 : 1e-4) ||
238+
lambda_min_ratio == (size(y, 1) < size(X, 2) ? 1e-2 : 1e-4) ||
239239
error("cannot specify both lambda and lambda_min_ratio")
240240
nlambda = length(lambda)
241241
lambda_min_ratio = 2.0
@@ -514,6 +514,7 @@ function show(io::IO, cv::GLMNetCrossValidation)
514514
end
515515

516516
include("Multinomial.jl")
517+
include("MultiResponseNet.jl")
517518
include("CoxNet.jl")
518519

519520
function glmnetcv(X::AbstractMatrix, y::Union{AbstractVector{<:Number},AbstractMatrix{<:Number}},
@@ -534,10 +535,10 @@ function glmnetcv(X::AbstractMatrix, y::Union{AbstractVector{<:Number},AbstractM
534535
y = convert(Array{Float64}, y)
535536
end
536537
# Fit full model once to determine parameters
537-
offsets = (offsets != nothing) ? offsets : isa(family, Multinomial) ? y*0.0 : zeros(size(X, 1))
538+
offsets = (offsets != nothing) ? offsets : isa(family, Union{Multinomial, MvNormal}) ? y*0.0 : zeros(size(X, 1))
538539

539540

540-
if isa(family, Normal)
541+
if isa(family, Union{Normal, MvNormal})
541542
path = glmnet(X, y, family; weights = weights, kw...)
542543
else
543544
path = glmnet(X, y, family; weights = weights, offsets = offsets, kw...)
@@ -560,7 +561,7 @@ function glmnetcv(X::AbstractMatrix, y::Union{AbstractVector{<:Number},AbstractM
560561
f = folds .== i
561562
holdoutidx = findall(f)
562563
modelidx = findall(!,f)
563-
if isa(family, Normal)
564+
if isa(family, Union{Normal, MvNormal})
564565
g = glmnet!(X[modelidx, :], isa(y, AbstractVector) ? y[modelidx] : y[modelidx, :], family;
565566
weights=weights[modelidx], lambda=path.lambda, kw...)
566567
else

‎src/MultiResponseNet.jl

+188
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,188 @@
1+
import Distributions.MvNormal
2+
3+
MvNormal() = MvNormal([0, 0], [1 0; 0 1]) # MvNormal(0, ([1;]))
4+
modeltype(::MvNormal) = "MvNormal"
5+
6+
7+
function predict(path::GLMNetPath{<:MvNormal}, X::AbstractMatrix,
8+
model::Union{Int, AbstractVector{Int}}=1:length(path.lambda);
9+
outtype = :link, offsets = zeros(size(X, 1), size(path.betas, 2)))
10+
11+
nresp = size(path.betas, 2);
12+
out = zeros(Float64, size(X, 1), nresp, length(model));
13+
for i = 1:length(model)
14+
out[:, :, i] = repeat(path.a0[:,model[i]]', size(X, 1)) + X * path.betas[:, :, model[i]] + offsets
15+
end
16+
if outtype != :link
17+
for i = 1:size(X, 1), j = 1:length(model)
18+
out = exp.(out)
19+
end
20+
end
21+
if length(model) == 1
22+
return out[:, :, 1]
23+
else
24+
return out
25+
end
26+
end
27+
28+
nactive(g::GLMNetPath{<:MvNormal}, b::AbstractVector{Int}=1:size(g.betas, 3)) =
29+
[nactive(g.betas, j, dims=2) for j in b]
30+
31+
function show(io::IO, g::GLMNetPath{<:MvNormal})
32+
println(io, "$(modeltype(g.family)) GLMNet Solution Path ($(size(g.betas, 3)) solutions for $(size(g.betas, 1)) predictors in $(g.npasses) passes):")
33+
print(io, DataFrame(df=nactive(g), pct_dev=g.dev_ratio, λ=g.lambda))
34+
end
35+
36+
struct MultiMSE <: Loss
37+
y::Matrix{Float64}
38+
end
39+
loss(l::MultiMSE, i, mu) = sum(abs2.(l.y[i,:] .- mu))
40+
41+
devloss(::MvNormal, y) = MultiMSE(y)
42+
43+
function loss(path::GLMNetPath{<:MvNormal}, X::AbstractMatrix{Float64},
44+
y::Union{AbstractVector{Float64}, AbstractMatrix{Float64}},
45+
weights::AbstractVector{Float64}=ones(size(y,1)),
46+
lossfun::Loss=devloss(path.family, y),
47+
model::Union{Int, AbstractVector{Int}}=1:length(path.lambda);
48+
offsets = zeros(size(X, 1), size(path.betas, 2)))
49+
50+
validate_x_y_weights(X, y, weights)
51+
mu = predict(path, X; offsets = offsets)
52+
devs = zeros(size(mu, 3))
53+
for j = 1:size(mu, 3), i = 1:size(mu, 1)
54+
devs[j] += loss(lossfun, i, vec(mu[i, :, j]))*weights[i]
55+
end
56+
devs/sum(weights)
57+
end
58+
59+
loss(path::GLMNetPath{<:MvNormal}, X::AbstractMatrix, y::Union{AbstractVector, AbstractMatrix},
60+
weights::AbstractVector=ones(size(y,1)), va...; kw...) =
61+
loss(path, convert(Matrix{Float64}, X),
62+
convert(Array{Float64}, y),
63+
convert(Vector{Float64}, weights), va...; kw...)
64+
65+
66+
67+
macro check_and_return_mvnormal()
68+
esc(quote
69+
check_jerr(jerr[], maxit,pmax)
70+
lmu = lmu_ref[]
71+
# first lambda is infinity; changed to entry point
72+
if isempty(lambda) && length(alm) > 2
73+
alm[1] = exp(2*log(alm[2])-log(alm[3]))
74+
end
75+
GLMNetPath(family, a0[:, 1:lmu], ca[sortperm(ia), :, 1:lmu],
76+
null_dev[], fdev[1:lmu], alm[1:lmu], Int(nlp[]))
77+
end)
78+
end
79+
80+
# change of parameters from elnet to multelnet
81+
# ka,parm,no,ni,nr, x,y,w,jd, vp,cl,ne,nx, nlam,flmin,ulam,thr, isd, ,intr,maxit,lmu, a0,ca,ia,nin, rsq,alm,nlp,jerr
82+
# parm,no,ni,nr, x,y,w,jd, vp,cl,ne,nx, nlam,flmin,ulam,thr, isd,jsd,intr,maxit,lmu, a0,ca,ia,nin, rsq,alm,nlp,jerr
83+
# multi-response normal
84+
function glmnet!(X::Matrix{Float64}, y::Matrix{Float64},
85+
family::MvNormal=MvNormal();
86+
weights::Vector{Float64}=ones(size(y,1)),
87+
naivealgorithm::Bool=(size(X, 2) >= 500), alpha::Real=1.0,
88+
penalty_factor::Vector{Float64}=ones(size(X, 2)),
89+
constraints::Array{Float64, 2}=[x for x in (-Inf, Inf), y in 1:size(X, 2)],
90+
dfmax::Int=size(X, 2), pmax::Int=min(dfmax*2+20, size(X, 2)), nlambda::Int=100,
91+
lambda_min_ratio::Real=(size(y, 1) < size(X, 2) ? 1e-2 : 1e-4),
92+
lambda::Vector{Float64}=Float64[], tol::Real=1e-7, standardize::Bool=true,
93+
standardize_response::Bool=false,
94+
intercept::Bool=true, maxit::Int=1000000)
95+
96+
@validate_and_init_multi
97+
standardize_response = Int32(standardize_response)
98+
99+
# Compute null deviance
100+
yw = y .* repeat(weights, 1, size(y, 2))
101+
mu = mean(y, dims=1)
102+
if intercept == 0
103+
mu = fill(intercept, 1, size(y,2))
104+
end
105+
# Sum of squared error (weighted by obervation weights)
106+
null_dev[] = sum(weights) .* mean(abs2.(yw .- mu))
107+
108+
109+
ccall((:multelnet_, libglmnet), Cvoid,
110+
111+
(Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Int32},
112+
Ref{Float64}, Ref{Float64}, Ref{Float64}, Ref{Int32},
113+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
114+
Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Float64},
115+
Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32},
116+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
117+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
118+
119+
alpha, nobs, nvars, nresp,
120+
X, y, weights, jd,
121+
penalty_factor, constraints, dfmax, pmax,
122+
nlambda, lambda_min_ratio, lambda, tol,
123+
standardize, standardize_response, intercept, maxit, lmu_ref,
124+
a0, ca, ia, nin,
125+
fdev, alm, nlp, jerr
126+
)
127+
128+
@check_and_return_mvnormal
129+
end
130+
131+
132+
133+
# multi-response sparse normal
134+
function glmnet!(X::AbstractSparseMatrix{Float64}, y::Matrix{Float64},
135+
family::MvNormal=MvNormal();
136+
weights::Vector{Float64}=ones(size(y, 1)),
137+
naivealgorithm::Bool=(size(X, 2) >= 500), alpha::Real=1.0,
138+
penalty_factor::Vector{Float64}=ones(size(X, 2)),
139+
constraints::Array{Float64, 2}=[x for x in (-Inf, Inf), y in 1:size(X, 2)],
140+
dfmax::Int=size(X, 2), pmax::Int=min(dfmax*2+20, size(X, 2)), nlambda::Int=100,
141+
lambda_min_ratio::Real=(size(y, 1) < size(X, 2) ? 1e-2 : 1e-4),
142+
lambda::Vector{Float64}=Float64[], tol::Real=1e-7, standardize::Bool=true,
143+
standardize_response::Bool=false,
144+
intercept::Bool=true, maxit::Int=1000000)
145+
146+
@validate_and_init_multi
147+
standardize_response = Int32(standardize_response)
148+
149+
# Compute null deviance
150+
yw = y .* repeat(weights, 1, size(y, 2))
151+
mu = mean(y, dims=1)
152+
if intercept == 0
153+
mu = fill(intercept, 1, size(y,2))
154+
end
155+
# Sum of squared error (weighted by obervation weights)
156+
null_dev[] = sum(weights) .* mean(abs2.(yw .- mu))
157+
158+
159+
ccall((:multspelnet_, libglmnet), Cvoid,
160+
161+
(Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Int32},
162+
Ref{Float64}, Ref{Int32}, Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Int32},
163+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
164+
Ref{Int32}, Ref{Float64}, Ref{Float64}, Ref{Float64},
165+
Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32},
166+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32},
167+
Ref{Float64}, Ref{Float64}, Ref{Int32}, Ref{Int32}),
168+
169+
alpha, nobs, nvars, nresp,
170+
X.nzval, X.colptr, X.rowval, y, weights, jd,
171+
penalty_factor, constraints, dfmax, pmax,
172+
nlambda, lambda_min_ratio, lambda, tol,
173+
standardize, standardize_response, intercept, maxit, lmu_ref,
174+
a0, ca, ia, nin,
175+
fdev, alm, nlp, jerr
176+
)
177+
178+
@check_and_return_mvnormal
179+
end
180+
181+
glmnet(X::Matrix{Float64}, y::Matrix{Float64}, family::MvNormal; kw...) =
182+
glmnet!(copy(X), copy(y), family; kw...)
183+
glmnet(X::SparseMatrixCSC{Float64,Int32}, y::Matrix{Float64}, family::MvNormal; kw...) =
184+
glmnet!(copy(X), copy(y), family; kw...)
185+
glmnet(X::AbstractMatrix, y::AbstractMatrix{<:Number}, family::MvNormal; kw...) =
186+
glmnet(convert(Matrix{Float64}, X), convert(Matrix{Float64}, y), family; kw...)
187+
glmnet(X::SparseMatrixCSC, y::AbstractMatrix{<:Number}, family::MvNormal; kw...) =
188+
glmnet(convert(SparseMatrixCSC{Float64,Int32}, X), convert(Matrix{Float64}, y), family; kw...)

‎src/Multinomial.jl

+3-3
Original file line numberDiff line numberDiff line change
@@ -62,14 +62,14 @@ loss(path::GLMNetPath{<:Multinomial}, X::AbstractMatrix, y::Union{AbstractVector
6262

6363
# Get number of active predictors for a model in X
6464
# nin can be > non-zero predictors under some circumstances...
65-
nactive(X::Array{Float64, 3}, b::Int) = sum(sum(X[:,:,b] .!= 0., dims=1) .> 0)
65+
nactive(X::Array{Float64, 3}, b::Int; dims=1) = sum(sum(X[:,:,b] .!= 0., dims=dims) .> 0)
6666

6767
nactive(X::Array{Float64, 3}, b::AbstractVector{Int}=1:size(X, 3)) =
6868
[nactive(X, j) for j in b]
6969

7070

7171
function show(io::IO, g::GLMNetPath{<:Multinomial})
72-
println(io, "$(modeltype(g.family)) GLMNet Solution Path ($(size(g.betas, 2)) solutions for $(size(g.betas, 1)) predictors in $(g.npasses) passes):")
72+
println(io, "$(modeltype(g.family)) GLMNet Solution Path ($(size(g.betas, 3)) solutions for $(size(g.betas, 1)) predictors in $(g.npasses) passes):")
7373
print(io, DataFrame(df=nactive(g.betas), pct_dev=g.dev_ratio, λ=g.lambda))
7474
end
7575

@@ -86,7 +86,7 @@ macro validate_and_init_multi()
8686
if !isempty(lambda)
8787
# user-specified lambda values
8888
nlambda == 100 || error("cannot specify both lambda and nlambda")
89-
lambda_min_ratio == (length(y) < size(X, 2) ? 1e-2 : 1e-4) ||
89+
lambda_min_ratio == (size(y, 1) < size(X, 2) ? 1e-2 : 1e-4) ||
9090
error("cannot specify both lambda and lambda_min_ratio")
9191
nlambda = length(lambda)
9292
lambda_min_ratio = 2.0

0 commit comments

Comments
 (0)