diff --git a/src/GLM.jl b/src/GLM.jl index 5401b8d1..accae615 100644 --- a/src/GLM.jl +++ b/src/GLM.jl @@ -93,19 +93,19 @@ abstract type DensePred <: LinPred end # linear predictor with dense abstract type LinPredModel <: RegressionModel end # model based on a linear predictor const COMMON_FIT_KWARGS_DOCS = """ - - `dropcollinear::Bool`: Controls whether or not a model matrix + - `dropcollinear::Bool=false`: Controls whether or not a model matrix less-than-full rank is accepted. If `true` (the default) the coefficient for redundant linearly dependent columns is `0.0` and all associated statistics are set to `NaN`. Typically from a set of linearly-dependent columns the last ones are identified as redundant (however, the exact selection of columns identified as redundant is not guaranteed). - - `method::Symbol`: Controls which decomposition method to use. + - `method::Symbol=:qr`: Controls which decomposition method to use. If `method=:qr` (the default), then the `QR` decomposition method will be used. If `method=:cholesky`, then the `Cholesky` decomposition method will be used. The Cholesky decomposition is faster and more computationally efficient than QR, but is less numerically stable and thus may fail or produce less accurate estimates for some models. - - `wts::AbstractWeights`: Weights of observations. + - `wts::AbstractWeights=uweights(0)`: Weights of observations. The weights can be of type `AnalyticWeights`, `FrequencyWeights`, `ProbabilityWeights`, or `UnitWeights`. `AnalyticWeights` describe a non-random relative importance (usually between 0 and 1) for each observation. These weights may @@ -114,7 +114,7 @@ const COMMON_FIT_KWARGS_DOCS = """ `ProbabilityWeights` represent the inverse of the sampling probability for each observation, providing a correction mechanism for under- or over-sampling certain population groups. `UnitWeights` (default) describes the case in which all weights are equal to 1 (so no weighting takes place) - - `contrasts::AbstractDict{Symbol}`: a `Dict` mapping term names + - `contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}()`: a `Dict` mapping term names (as `Symbol`s) to term types (e.g., `ContinuousTerm`) or contrasts (e.g., `HelmertCoding()`, `SeqDiffCoding(; levels=["a", "b", "c"])`, etc.). If contrasts are not provided for a variable, the appropriate diff --git a/src/glmfit.jl b/src/glmfit.jl index 33193eee..f2e59dfe 100644 --- a/src/glmfit.jl +++ b/src/glmfit.jl @@ -381,7 +381,7 @@ end dof(obj::GeneralizedLinearModel) = linpred_rank(obj) + dispersion_parameter(obj.rr.d) function _fit!(m::AbstractGLM, maxiter::Integer, minstepfac::Real, - atol::Real, rtol::Real, start) + atol::Real, rtol::Real, start::Union{AbstractVector,Nothing}) # Return early if model has the fit flag set m.fit && return m @@ -463,7 +463,7 @@ function StatsBase.fit!(m::AbstractGLM; minstepfac::Real=0.001, atol::Real=1e-6, rtol::Real=1e-6, - start=nothing, + start::Union{AbstractVector,Nothing}=nothing, kwargs...) if haskey(kwargs, :verbose) Base.depwarn("""`verbose` argument is deprecated, use `ENV["JULIA_DEBUG"]=GLM` instead.""", @@ -483,7 +483,7 @@ end const FIT_GLM_DOC = """ In the first method, `formula` must be a - [StatsModels.jl `Formula` object](https://juliastats.org/StatsModels.jl/stable/formula/) + [StatsModels.jl `FormulaTerm` object](https://juliastats.org/StatsModels.jl/stable/formula/) and `data` a table (in the [Tables.jl](https://tables.juliadata.org/stable/) definition, e.g., a data frame). In the second method, `X` must be a matrix holding values of the independent variable(s) in columns (including if appropriate the intercept), and `y` must be a vector holding @@ -494,20 +494,20 @@ const FIT_GLM_DOC = """ # Keyword Arguments $COMMON_FIT_KWARGS_DOCS - - `offset::Vector=similar(y,0)`: offset added to `Xβ` to form `eta`. Can be of - length 0 + - `offset::Union{AbstractVector{<:Real},Nothing}=nothing`: offset added to `Xβ` to form `eta`. + Can be of length 0. - `maxiter::Integer=30`: Maximum number of iterations allowed to achieve convergence - `atol::Real=1e-6`: Convergence is achieved when the relative change in deviance is less than `max(rtol*dev, atol)`. - `rtol::Real=1e-6`: Convergence is achieved when the relative change in deviance is less than `max(rtol*dev, atol)`. - `minstepfac::Real=0.001`: Minimum line step fraction. Must be between 0 and 1. - - `start::AbstractVector=nothing`: Starting values for beta. Should have the + - `start::Union{AbstractVector,Nothing}=nothing`: Starting values for beta. Should have the same length as the number of columns in the model matrix. """ """ - fit(GeneralizedLinearModel, formula, data, + fit(GeneralizedLinearModel, formula::FormulaTerm, data, distr::UnivariateDistribution, link::Link = canonicallink(d); ) fit(GeneralizedLinearModel, X::AbstractMatrix, y::AbstractVector, distr::UnivariateDistribution, link::Link = canonicallink(d); ) @@ -524,15 +524,22 @@ function fit(::Type{M}, dropcollinear::Bool=true, method::Symbol=:qr, wts::Union{AbstractWeights,AbstractVector{<:Real}}=uweights(length(y)), - offset::AbstractVector{<:Real}=similar(y, 0), - fitargs...) where {M<:AbstractGLM} + offset::Union{AbstractVector{<:Real},Nothing}=nothing, + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing) where {M<:AbstractGLM} # Check that X and y have the same number of observations if size(X, 1) != size(y, 1) throw(DimensionMismatch("number of rows in X and y must match")) end _wts = convert_weights(wts, length(y)) - rr = GlmResp(y, d, l, offset, _wts) + + off = offset === nothing ? similar(y, 0) : offset + + rr = GlmResp(y, d, l, off, _wts) if method === :cholesky res = M(rr, cholpred(X, dropcollinear, _wts), nothing, false) @@ -542,15 +549,26 @@ function fit(::Type{M}, throw(ArgumentError("The only supported values for keyword argument `method` are `:cholesky` and `:qr`.")) end - return fit!(res; fitargs...) + return fit!(res; maxiter, atol, rtol, minstepfac, start) end function fit(::Type{M}, X::AbstractMatrix, y::AbstractVector, d::UnivariateDistribution, - l::Link=canonicallink(d); kwargs...) where {M<:AbstractGLM} - return fit(M, float(X), float(y), d, l; kwargs...) + l::Link=canonicallink(d); + wts::AbstractVector{<:Real}=uweights(length(y)), + offset::Union{AbstractVector{<:Real},Nothing}=nothing, + dropcollinear::Bool=true, + method::Symbol=:qr, + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing) where {M<:AbstractGLM} + return fit(M, float(X), float(y), d, l; + offset, wts, dropcollinear, method, + maxiter, atol, rtol, minstepfac, start) end function fit(::Type{M}, @@ -559,13 +577,16 @@ function fit(::Type{M}, d::UnivariateDistribution, l::Link=canonicallink(d); offset::Union{AbstractVector,Nothing}=nothing, - wts::Union{AbstractVector,Nothing}=nothing, + wts::AbstractVector{<:Real}=uweights(0), dropcollinear::Bool=true, method::Symbol=:qr, contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}(), - fitargs...) where {M<:AbstractGLM} + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing) where {M<:AbstractGLM} f, (y, X) = modelframe(f, data, contrasts, M) - wts = wts === nothing ? uweights(0) : wts _wts = convert_weights(wts, length(y)) # Check that X and y have the same number of observations if size(X, 1) != size(y, 1) @@ -584,11 +605,11 @@ function fit(::Type{M}, throw(ArgumentError("The only supported values for keyword argument `method` are `:cholesky` and `:qr`.")) end - return fit!(res; fitargs...) + return fit!(res; maxiter, atol, rtol, minstepfac, start) end """ - glm(formula, data, + glm(formula::FormulaTerm, data, distr::UnivariateDistribution, link::Link = canonicallink(distr); ) glm(X::AbstractMatrix, y::AbstractVector, distr::UnivariateDistribution, link::Link = canonicallink(distr); ) @@ -597,7 +618,35 @@ Fit a generalized linear model to data. Alias for `fit(GeneralizedLinearModel, . $FIT_GLM_DOC """ -glm(X, y, args...; kwargs...) = fit(GeneralizedLinearModel, X, y, args...; kwargs...) +glm(X::AbstractMatrix, y::AbstractVector, + d::UnivariateDistribution, l::Link=canonicallink(d); + offset::Union{AbstractVector{<:Real},Nothing}=nothing, + wts::AbstractVector{<:Real}=uweights(length(y)), + dropcollinear::Bool=true, + method::Symbol=:qr, + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing) = + fit(GeneralizedLinearModel, X, y, d, l; + offset, wts, dropcollinear, method, + maxiter, atol, rtol, minstepfac, start) + +glm(formula::FormulaTerm, data, d::UnivariateDistribution, l::Link=canonicallink(d); + offset::Union{AbstractVector{<:Real},Nothing}=nothing, + wts::AbstractVector{<:Real}=uweights(0), + dropcollinear::Bool=true, + method::Symbol=:qr, + contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}(), + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing) = + fit(GeneralizedLinearModel, formula, data, d, l; + offset, wts, dropcollinear, method, contrasts, + maxiter, atol, rtol, minstepfac, start) GLM.Link(r::GlmResp) = r.link GLM.Link(m::GeneralizedLinearModel) = Link(m.rr) diff --git a/src/lm.jl b/src/lm.jl index 79be796b..7c14044e 100644 --- a/src/lm.jl +++ b/src/lm.jl @@ -148,19 +148,13 @@ const FIT_LM_DOC = """ """ """ - fit(LinearModel, formula::FormulaTerm, data; - wts::AbstractWeights=uweights(0), - dropcollinear::Bool=true, method::Symbol=:qr, - contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}()) - fit(LinearModel, X::AbstractMatrix, y::AbstractVector; - wts::AbstractWeights=uweights(length(y)), - dropcollinear::Bool=true, method::Symbol=:qr) + fit(LinearModel, formula::FormulaTerm, data; ) + fit(LinearModel, X::AbstractMatrix, y::AbstractVector; ) Fit a linear model to data. $FIT_LM_DOC """ - function fit(::Type{LinearModel}, X::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}; wts::Union{AbstractWeights,AbstractVector{<:Real}}=uweights(length(y)), dropcollinear::Bool=true, method::Symbol=:qr) @@ -191,18 +185,26 @@ function fit(::Type{LinearModel}, f::FormulaTerm, data; end """ - lm(formula, data; - [wts::AbstractVector], dropcollinear::Bool=true, method::Symbol=:qr, - contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}()) - lm(X::AbstractMatrix, y::AbstractVector; - wts::AbstractVector=similar(y, 0), dropcollinear::Bool=true, method::Symbol=:cholesky) + lm(formula::FormulaTerm, data; ) + lm(X::AbstractMatrix, y::AbstractVector; ) Fit a linear model to data. -An alias for `fit(LinearModel, X, y; wts=wts, dropcollinear=dropcollinear, method=method)` +An alias for `fit(LinearModel, ...)`. $FIT_LM_DOC """ -lm(X, y; kwargs...) = fit(LinearModel, X, y; kwargs...) +lm(X::AbstractMatrix, y::AbstractVector; + wts::Union{AbstractWeights,AbstractVector{<:Real}}=uweights(length(y)), + dropcollinear::Bool=true, + method::Symbol=:qr) = + fit(LinearModel, X, y; wts, dropcollinear, method) + +lm(f::FormulaTerm, data; + wts::Union{AbstractWeights,AbstractVector{<:Real}}=uweights(0), + dropcollinear::Bool=true, + method::Symbol=:qr, + contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}()) = + fit(LinearModel, f, data; wts, dropcollinear, method, contrasts) dof(x::LinearModel) = linpred_rank(x.pp) + 1 diff --git a/src/negbinfit.jl b/src/negbinfit.jl index 9f78d388..2d72fbf0 100644 --- a/src/negbinfit.jl +++ b/src/negbinfit.jl @@ -44,14 +44,13 @@ function mle_for_θ(y::AbstractVector, μ::AbstractVector, wts::AbstractWeights; end """ - negbin(formula, data, [link::Link]; + negbin(formula::FormulaTerm, data, link::Union{Link,Nothing}=nothing; ) - negbin(X::AbstractMatrix, y::AbstractVector, [link::Link]; + negbin(X::AbstractMatrix, y::AbstractVector, link::Union{Link,Nothing}=nothing; ) Fit a negative binomial generalized linear model to data, while simultaneously -estimating the shape parameter θ. Extra arguments and keyword arguments will be -passed to [`glm`](@ref). +estimating the shape parameter θ. Arguments are the same as for [`glm`](@ref) except `initialθ`. In the first method, `formula` must be a [StatsModels.jl `Formula` object](https://juliastats.org/StatsModels.jl/stable/formula/) @@ -59,30 +58,72 @@ and `data` a table (in the [Tables.jl](https://tables.juliadata.org/stable/) def In the second method, `X` must be a matrix holding values of the independent variable(s) in columns (including if appropriate the intercept), and `y` must be a vector holding values of the dependent variable. -In both cases, `link` may specify the link function -(if omitted, it is taken to be `NegativeBinomial(θ)`). + +In both cases, `link` may specify the link function. If omitted, it is taken to be +`NegativeBinomialLink(θ)`. # Keyword Arguments - `initialθ::Real=Inf`: Starting value for shape parameter θ. If it is `Inf` then the initial value will be estimated by fitting a Poisson distribution. -- `dropcollinear::Bool=true`: See `dropcollinear` for [`glm`](@ref) -- `method::Symbol=:qr`: See `method` for [`glm`](@ref) -- `maxiter::Integer=30`: See `maxiter` for [`glm`](@ref) -- `atol::Real=1.0e-6`: See `atol` for [`glm`](@ref) -- `rtol::Real=1.0e-6`: See `rtol` for [`glm`](@ref) +$COMMON_FIT_KWARGS_DOCS +- `offset::Union{AbstractVector{<:Real},Nothing}=nothing,`: offset added to `Xβ` + to form `eta`. Can be of length 0. +- `maxiter::Integer=30`: Maximum number of iterations allowed to achieve convergence +- `atol::Real=1e-6`: Convergence is achieved when the relative change in + deviance is less than `max(rtol*dev, atol)`. +- `rtol::Real=1e-6`: Convergence is achieved when the relative change in + deviance is less than `max(rtol*dev, atol)`. +- `minstepfac::Real=0.001`: Minimum line step fraction. Must be between 0 and 1. +- `start::Union{AbstractVector,Nothing}=nothing`: Starting values for beta. Should have the + same length as the number of columns in the model matrix. """ -function negbin(F, - D, - args...; - wts::Union{AbstractWeights,AbstractVector{<:Real}}=uweights(0), - initialθ::Real=Inf, - dropcollinear::Bool=true, - method::Symbol=:qr, - maxiter::Integer=30, - minstepfac::Real=0.001, - atol::Real=1e-6, - rtol::Real=1.e-6, - kwargs...) +negbin(X::AbstractMatrix, y::AbstractVector, l::Union{Link,Nothing}=nothing; + initialθ::Real=Inf, + offset::Union{AbstractVector{<:Real},Nothing}=nothing, + wts::AbstractVector{<:Real}=uweights(length(y)), + dropcollinear::Bool=true, + method::Symbol=:qr, + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing, + kwargs...) = + _negbin(X, y, l; + initialθ, offset, wts, dropcollinear, method, contrasts=nothing, + maxiter, atol, rtol, minstepfac, start, kwargs...) +negbin(formula::FormulaTerm, data, l::Union{Link,Nothing}=nothing; + initialθ::Real=Inf, + offset::Union{AbstractVector{<:Real},Nothing}=nothing, + wts::AbstractVector{<:Real}=uweights(0), + dropcollinear::Bool=true, + method::Symbol=:qr, + contrasts::AbstractDict{Symbol}=Dict{Symbol,Any}(), + maxiter::Integer=30, + atol::Real=1e-6, + rtol::Real=1e-6, + minstepfac::Real=0.001, + start::Union{AbstractVector,Nothing}=nothing, + kwargs...) = + _negbin(formula, data, l; + initialθ, offset, wts, dropcollinear, method, contrasts, + maxiter, atol, rtol, minstepfac, start, kwargs...) + +function _negbin(F, + D, + l::Union{Link,Nothing}; + initialθ::Real, + offset::Union{AbstractVector{<:Real},Nothing}, + wts::AbstractVector, + dropcollinear::Bool, + method::Symbol, + contrasts::Union{AbstractDict{Symbol},Nothing}, + maxiter::Integer, + atol::Real, + rtol::Real, + minstepfac::Real, + start::Union{AbstractVector,Nothing}, + kwargs...) if haskey(kwargs, :verbose) Base.depwarn("""`verbose` argument is deprecated, use `ENV["JULIA_DEBUG"]=GLM` instead.""", :negbin) @@ -96,20 +137,18 @@ function negbin(F, rtol > 0 || throw(ArgumentError("rtol must be positive")) initialθ > 0 || throw(ArgumentError("initialθ must be positive")) + contrasts_kwarg = isnothing(contrasts) ? () : (contrasts=contrasts,) + # fit a Poisson regression model if the user does not specify an initial θ - if isinf(initialθ) - regmodel = glm(F, D, Poisson(), args...; - wts=wts, dropcollinear=dropcollinear, method=method, maxiter=maxiter, - atol=atol, rtol=rtol, kwargs...) - else - regmodel = glm(F, D, NegativeBinomial(initialθ), args...; - wts=wts, dropcollinear=dropcollinear, method=method, maxiter=maxiter, - atol=atol, rtol=rtol, kwargs...) - end + distr = isinf(initialθ) ? Poisson() : NegativeBinomial(initialθ) + regmodel = glm(F, D, distr, something(l, canonicallink(distr)); + offset, wts, dropcollinear, method, contrasts, + maxiter, atol, rtol, minstepfac, start, contrasts_kwarg...) μ = regmodel.rr.mu y = regmodel.rr.y wts = regmodel.rr.wts + lw, ly = length(wts), length(y) if lw != ly throw(ArgumentError("length of `wts` must be $ly but was $lw")) @@ -128,9 +167,9 @@ function negbin(F, break end @debug "NegativeBinomial dispersion optimization" iteration = i θ = θ - regmodel = glm(F, D, NegativeBinomial(θ), args...; - dropcollinear=dropcollinear, method=method, maxiter=maxiter, - atol=atol, rtol=rtol, kwargs...) + regmodel = glm(F, D, NegativeBinomial(θ), something(l, NegativeBinomialLink(θ)); + offset, wts, dropcollinear, method, contrasts, + maxiter, atol, rtol, minstepfac, start, contrasts_kwarg...) μ = regmodel.rr.mu prevθ = θ θ = mle_for_θ(y, μ, wts; maxiter=maxiter, tol=rtol) diff --git a/test/probability_weights.jl b/test/probability_weights.jl index dbcfd563..955b89a4 100644 --- a/test/probability_weights.jl +++ b/test/probability_weights.jl @@ -172,7 +172,7 @@ end @testset "GLM: NegativeBinomial(1) with LogLink link - ProbabilityWeights with $dmethod method with dropcollinear=$drop" for (dmethod, drop) in - itr + itr model = glm(@formula(Days ~ Eth + Sex + Age + Lrn), quine, @@ -221,7 +221,7 @@ end @test stderror(model) ≈ stderror(model_geom) rtol = 1e-06 end -@testset "GLM: NegaiveBinomial(2) with SqrtLink link - ProbabilityWeights with $dmethod method with dropcollinear=$drop" for (dmethod, +@testset "GLM: NegativeBinomial(2) with SqrtLink link - ProbabilityWeights with $dmethod method with dropcollinear=$drop" for (dmethod, drop) in itr diff --git a/test/runtests.jl b/test/runtests.jl index 86560d3c..f9a83bec 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -872,6 +872,23 @@ end end end end + @testset "negbin with contrasts" begin + gm20b = negbin(@formula(Days ~ Eth + Sex + Age + Lrn), quine, LogLink(), + contrasts=Dict(:Sex => EffectsCoding())) + test_show(gm20b) + @test dof(gm20b) == 8 + @test isapprox(deviance(gm20b), 167.9518430624193, rtol=1e-7) + @test isapprox(nulldeviance(gm20b), 195.28668602703388, rtol=1e-7) + @test isapprox(loglikelihood(gm20b), -546.57550938017, rtol=1e-7) + @test isapprox(nullloglikelihood(gm20b), -560.2429308624774, rtol=1e-7) + @test isapprox(aic(gm20b), 1109.15101876034) + @test isapprox(aicc(gm20b), 1110.202113650851) + @test isapprox(bic(gm20b), 1133.0198717340068) + @test isapprox(coef(gm20b)[1:7], + [2.9357217632468595, -0.5693411448715979, + 0.041194065435350515, -0.4484636623590206, + 0.08805060372902418, 0.3569553124412582, 0.2921383118842893]) + end end @testset "Weighted NegativeBinomial LogLink, θ to be estimated with Cholesky" begin @@ -881,17 +898,17 @@ end wts=fweights(wts)) test_show(gm20a) @test dof(gm20a) == 8 - @test isapprox(deviance(gm20a), 164.45910399188858, rtol=1e-7) - @test isapprox(nulldeviance(gm20a), 191.14269166948384, rtol=1e-7) - @test isapprox(loglikelihood(gm20a), -546.596822900127, rtol=1e-7) - @test isapprox(nullloglikelihood(gm20a), -559.9386167389254, rtol=1e-7) - @test isapprox(aic(gm20a), 1109.193645800254) - @test isapprox(aicc(gm20a), 1110.244740690765) - @test isapprox(bic(gm20a), 1133.0624987739207) + @test isapprox(deviance(gm20a), 168.40402933035944, rtol=1e-7) + @test isapprox(nulldeviance(gm20a), 196.50242701899307, rtol=1e-7) + @test isapprox(loglikelihood(gm20a), -537.7760823254398, rtol=1e-7) + @test isapprox(nullloglikelihood(gm20a), -551.8252811697569, rtol=1e-7) + @test isapprox(aic(gm20a), 1091.5521646508796) + @test isapprox(aicc(gm20a), 1092.6032595413906) + @test isapprox(bic(gm20a), 1115.4210176245463) @test isapprox(coef(gm20a)[1:7], - [2.894916710026395, -0.5694300339439156, - 0.08215779733345588, -0.44861865904551734, - 0.08783288494046998, 0.3568327292046044, 0.29190920267019166]) + [2.9531960459074083, -0.5991180517700114, + 0.09609289230162256, -0.48822247493829657, + 0.010246721844156351, 0.362004855234784, 0.24470519461989926]) end @testset "NegativeBinomial LogLink, θ to be estimated with QR" begin