Stop considering that Negative Binomial has a dispersion parameter

docs/src/examples.md

@@ -268,34 +268,34 @@ GeneralizedLinearModel
 Days ~ 1 + Eth + Sex + Age + Lrn
 
 Coefficients:
-────────────────────────────────────────────────────────────────────────────
-                  Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
-────────────────────────────────────────────────────────────────────────────
-(Intercept)   2.88645      0.227144  12.71    <1e-36   2.44125     3.33164
-Eth: N       -0.567515     0.152449  -3.72    0.0002  -0.86631    -0.26872
-Sex: M        0.0870771    0.159025   0.55    0.5840  -0.224606    0.398761
-Age: F1      -0.445076     0.239087  -1.86    0.0627  -0.913678    0.0235251
-Age: F2       0.0927999    0.234502   0.40    0.6923  -0.366816    0.552416
-Age: F3       0.359485     0.246586   1.46    0.1449  -0.123814    0.842784
-Lrn: SL       0.296768     0.185934   1.60    0.1105  -0.0676559   0.661191
-────────────────────────────────────────────────────────────────────────────
+─────────────────────────────────────────────────────────────────────────────
+                  Coef.  Std. Error      z  Pr(>|z|)    Lower 95%   Upper 95%
+─────────────────────────────────────────────────────────────────────────────
+(Intercept)   2.88645      0.186483  15.48    <1e-53   2.52095      3.25195
+Eth: N       -0.567515     0.125159  -4.53    <1e-05  -0.812823    -0.322207
+Sex: M        0.0870771    0.130558   0.67    0.5048  -0.168812     0.342966
+Age: F1      -0.445076     0.196288  -2.27    0.0234  -0.829794    -0.0603587
+Age: F2       0.0927999    0.192524   0.48    0.6298  -0.284541     0.470141
+Age: F3       0.359485     0.202445   1.78    0.0758  -0.0372989    0.756269
+Lrn: SL       0.296768     0.15265    1.94    0.0519  -0.00242092   0.595956
+─────────────────────────────────────────────────────────────────────────────
 
 julia> nbrmodel = negbin(@formula(Days ~ Eth+Sex+Age+Lrn), quine, LogLink())
-GeneralizedLinearModel
+NegativeBinomialModel
 
 Days ~ 1 + Eth + Sex + Age + Lrn
 
 Coefficients:
 ────────────────────────────────────────────────────────────────────────────
                   Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
 ────────────────────────────────────────────────────────────────────────────
-(Intercept)   2.89453      0.227415  12.73    <1e-36   2.4488      3.34025
-Eth: N       -0.569341     0.152656  -3.73    0.0002  -0.868541   -0.270141
-Sex: M        0.0823881    0.159209   0.52    0.6048  -0.229655    0.394431
-Age: F1      -0.448464     0.238687  -1.88    0.0603  -0.916281    0.0193536
-Age: F2       0.0880506    0.235149   0.37    0.7081  -0.372834    0.548935
-Age: F3       0.356955     0.247228   1.44    0.1488  -0.127602    0.841513
-Lrn: SL       0.292138     0.18565    1.57    0.1156  -0.0717297   0.656006
+(Intercept)   2.89453      0.228424  12.67    <1e-36   2.44683     3.34223
+Eth: N       -0.569341     0.153333  -3.71    0.0002  -0.869869   -0.268813
+Sex: M        0.0823881    0.159915   0.52    0.6064  -0.23104     0.395816
+Age: F1      -0.448464     0.239746  -1.87    0.0614  -0.918357    0.0214296
+Age: F2       0.0880506    0.236193   0.37    0.7093  -0.374879    0.55098
+Age: F3       0.356955     0.248325   1.44    0.1506  -0.129752    0.843663
+Lrn: SL       0.292138     0.186474   1.57    0.1172  -0.0733444   0.657621
 ────────────────────────────────────────────────────────────────────────────
 
 julia> println("Estimated theta = ", round(nbrmodel.rr.d.r, digits=5))

src/GLM.jl

@@ -55,8 +55,9 @@ export
 # Model types
       GeneralizedLinearModel,
       LinearModel,
+      NegativeBinomialModel,
 
-# functions
+# Functions
       canonicallink,  # canonical link function for a distribution
       deviance,       # deviance of fitted and observed responses
       devresid,       # vector of squared deviance residuals

src/glmfit.jl

@@ -299,7 +299,7 @@ end
 
 deviance(m::AbstractGLM) = deviance(m.rr)
 
-function nulldeviance(m::GeneralizedLinearModel)
+function nulldeviance(m::AbstractGLM)
     r = m.rr
     wts = weights(r)
     y = r.y
@@ -359,7 +359,7 @@ function loglikelihood(m::AbstractGLM)
     return ll
 end
 
-function nullloglikelihood(m::GeneralizedLinearModel)
+function nullloglikelihood(m::AbstractGLM)
     r = m.rr
     wts = weights(m)
     sumwt = sum(wts)
@@ -394,7 +394,15 @@ function nullloglikelihood(m::GeneralizedLinearModel)
     return ll
 end
 
-dof(obj::GeneralizedLinearModel) = linpred_rank(obj) + dispersion_parameter(obj.rr.d)
+"""
+    dof(model::AbstractGLM)
+
+For generalized linear models, consumed degrees of freedom correspond
+to the number of estimated coefficients, plus one for the estimated
+dispersion parameter if the distribution is `Gamma`, `InverseGaussian` or `Normal`
+(see [GLM.dispersion_parameter](@ref)).
+"""
+dof(obj::AbstractGLM) = linpred_rank(obj) + dispersion_parameter(obj.rr.d)
 
 function _fit!(m::AbstractGLM, maxiter::Integer, minstepfac::Real,
                atol::Real, rtol::Real, start::Union{AbstractVector,Nothing})
@@ -668,14 +676,20 @@ function glm(formula::FormulaTerm, data, d::UnivariateDistribution,
 end
 
 GLM.Link(r::GlmResp) = r.link
-GLM.Link(m::GeneralizedLinearModel) = Link(m.rr)
+GLM.Link(m::AbstractGLM) = Link(m.rr)
 
 """
+    dispersion(m::LinearModel, sqr::Bool=false)
     dispersion(m::AbstractGLM, sqr::Bool=false)
 
 Return the estimated dispersion (or scale) parameter for a model's distribution,
 generally written σ for linear models and ϕ for generalized linear models.
-It is, by definition, equal to 1 for the Bernoulli, Binomial, and Poisson families.
+
+It is, by definition, equal to 1 for the Bernoulli, Binomial, Geometric, Negative
+Binomial and Poisson families (see [`GLM.dispersion_parameter`](@ref)).
+For other distributions, it is estimated as the square root of the sum of Pearson
+residuals (i.e. for a linear model, the residual residual sum of squares)
+divided by the residual degrees of freedom.
 
 If `sqr` is `true`, the squared dispersion parameter is returned.
 """
@@ -881,7 +895,7 @@ function residuals(r::GlmResp; weighted::Bool=false)
     return dres
 end
 
-function momentmatrix(m::GeneralizedLinearModel)
+function momentmatrix(m::AbstractGLM)
     X = modelmatrix(m; weighted=false)
     r = varstruct(m)
     if link(m) isa Union{Gamma,InverseGaussian}
@@ -890,7 +904,7 @@ function momentmatrix(m::GeneralizedLinearModel)
     return Diagonal(r) * X
 end
 
-function varstruct(x::GeneralizedLinearModel)
+function varstruct(x::AbstractGLM)
     wrkwts = working_weights(x)
     wts = weights(x)
     wrkres = working_residuals(x)
@@ -901,20 +915,20 @@ function varstruct(x::GeneralizedLinearModel)
     end
 end
 
-function invloglikhessian(m::GeneralizedLinearModel)
+function invloglikhessian(m::AbstractGLM)
     r = varstruct(m)
     wts = weights(m)
     return inverse(m.pp) * sum(wts) / nobs(m)
 end
 
-function StatsBase.cooksdistance(m::GeneralizedLinearModel)
+function StatsBase.cooksdistance(m::AbstractGLM)
     h = leverage(m)
     hh = h ./ (1 .- h) .^ 2
     Rp = pearson_residuals(m)
     return (Rp .^ 2) .* hh ./ (dispersion(m)^2 * dof(m))
 end
 
-function pearson_residuals(m::GeneralizedLinearModel)
+function pearson_residuals(m::AbstractGLM)
     y = m.rr.y
     μ = predict(m)
     v = glmvar.(m.rr.d, μ)

src/glmtools.jl

@@ -452,8 +452,8 @@ devresid(::Poisson, y, μ::Real) = 2 * (xlogy(y, y / μ) - (y - μ))
 
 Does distribution `D` have a separate dispersion parameter, ϕ?
 
-Returns `true` for `Gamma`, `InverseGaussian`, `NegativeBinomial`
-and `Normal` distributions, and false for other known distributions.
+Returns `true` for `Gamma`, `InverseGaussian` and `Normal` distributions,
+and false for other known distributions.
 
 # Examples
 ```jldoctest; setup = :(using GLM)
@@ -464,8 +464,8 @@ julia> GLM.dispersion_parameter(Bernoulli())
 false
 ```
 """
-dispersion_parameter(::Union{Gamma,InverseGaussian,NegativeBinomial,Normal}) = true
-dispersion_parameter(::Union{Bernoulli,Binomial,Geometric,Poisson}) = false
+dispersion_parameter(::Union{Gamma,InverseGaussian,Normal}) = true
+dispersion_parameter(::Union{Bernoulli,Binomial,Geometric,NegativeBinomial,Poisson}) = false
 
 """
     _safe_int(x::T)

src/lm.jl

@@ -208,6 +208,13 @@ function lm(f::FormulaTerm, data;
     return fit(LinearModel, f, data; wts, dropcollinear, method, contrasts)
 end
 
+"""
+    dof(model::GeneralizedLinearModel)
+
+For linear models, consumed degrees of freedom correspond
+to the number of estimated coefficients, plus one for the estimated
+dispersion parameter σ.
+"""
 dof(x::LinearModel) = linpred_rank(x.pp) + 1
 
 """

src/negbinfit.jl

@@ -1,3 +1,25 @@
+mutable struct NegativeBinomialModel{G<:GlmResp,L<:LinPred} <: AbstractGLM
+    rr::G
+    pp::L
+    formula::Union{FormulaTerm,Nothing}
+    fit::Bool
+    maxiter::Int
+    minstepfac::Float64
+    atol::Float64
+    rtol::Float64
+end
+
+"""
+    dof(model::NegativeBinomialModel)
+
+For negative binomial models (fitted with [`negbin`](@ref)), consumed degrees
+of freedom correspond to the number of estimated coefficients plus one for the
+estimated shape parameter θ.
+"""
+function dof(obj::NegativeBinomialModel)
+    return linpred_rank(obj) + 1
+end
+
 function mle_for_θ(y::AbstractVector, μ::AbstractVector, wts::AbstractWeights;
                    maxiter=30, tol=1.e-6)
     function first_derivative(θ::Real)
@@ -180,5 +202,8 @@ function _negbin(F,
         ll = loglikelihood(regmodel)
     end
     converged || throw(ConvergenceException(maxiter))
-    return regmodel
+
+    return NegativeBinomialModel(regmodel.rr, regmodel.pp, regmodel.formula, regmodel.fit,
+                                 regmodel.maxiter, regmodel.minstepfac, regmodel.atol,
+                                 regmodel.rtol)
 end

test/analytic_weights.jl

@@ -334,16 +334,16 @@ end
     @test coef(model) ≈ [3.02411915515531, -0.4641576651688563, 0.0718560942992554,
                          -0.47848540911607984, 0.09677889908013552, 0.3562972562034356,
                          0.3480161821981514] rtol = 1e-07
-    @test stderror(model) ≈ [0.1950707397084349, 0.13200639191036218, 0.1373161597645507,
-                             0.2088476016141468, 0.20252412726336674, 0.21060778935484836,
-                             0.16126722793064027] rtol = 1e-07
-    # Tests below are broken because dof(model)==8 instead of 7.
-    # GLM.jl computes dof = linpred_rank + dispersion_parameter, where dispersion_parameter
-    # is 1 for NegativeBinomial (counting the shape parameter θ). This differs from R's
-    # calculation for weighted models. Needs investigation to determine if the discrepancy
-    # is in GLM.jl's dof calculation or in the expected R values.
-    @test_broken aic(model) ≈ 4023.1878928645556 rtol = 1e-07
-    @test_broken bic(model) ≈ 4044.073139216514 rtol = 1e-07
+    # Values match summary(model, dispersion=1) in R with
+    # model = glm(..., family=negative.binomial(1))
+    # but by default summary.glm estimates the dispersion instead of fixing
+    # it to 1 as it should
+    @test stderror(model) ≈ [0.09794735860165615, 0.06628199301176854, 0.06894809115074668,
+                             0.10486488624060987, 0.10168979390861807, 0.10574869762156038,
+                             0.08097420980935782] rtol = 1e-07
+    @test dof(model) == 7
+    @test aic(model) ≈ 4023.1878928645556 rtol = 1e-07
+    @test bic(model) ≈ 4044.073139216514 rtol = 1e-07
     @test GLM.momentmatrix(model) ≈
           [-3.866780529709063 -0.0 -3.866780529709063 -0.0 -0.0 -0.0 -3.866780529709063
            -4.370085797122667 -0.0 -4.370085797122667 -0.0 -0.0 -0.0 -4.370085797122667
@@ -505,11 +505,17 @@ end
     @test coef(model) ≈ [3.02411915515531, -0.4641576651688563, 0.0718560942992554,
                          -0.47848540911607984, 0.09677889908013552, 0.3562972562034356,
                          0.3480161821981514] rtol = 1e-07
-    @test stderror(model) ≈ [0.1950707397084349, 0.13200639191036218, 0.1373161597645507,
-                             0.2088476016141468, 0.20252412726336674, 0.21060778935484836,
-                             0.16126722793064027] rtol = 1e-07
-    @test_broken aic(model) ≈ 4023.1878928645556 rtol = 1e-07
-    @test_broken bic(model) ≈ 4044.073139216514 rtol = 1e-07
+    @test dof(model) == 7
+    # Values match summary(model, dispersion=1) in R with
+    # model = glm(..., family=negative.binomial(2))
+    # but by default summary.glm estimates the dispersion instead of fixing
+    # it to 1 as it should
+    @test stderror(model) ≈ [0.09794735860165615, 0.06628199301176854,
+                             0.06894809115074668, 0.10486488624060987,
+                             0.10168979390861807, 0.10574869762156038,
+                             0.08097420980935782] rtol = 1e-07
+    @test aic(model) ≈ 4023.1878928645556 rtol = 1e-07
+    @test bic(model) ≈ 4044.073139216514 rtol = 1e-07
     @test GLM.momentmatrix(model) ≈
           [-3.866780529709063 -0.0 -3.866780529709063 -0.0 -0.0 -0.0 -3.866780529709063
            -4.370085797122667 -0.0 -4.370085797122667 -0.0 -0.0 -0.0 -4.370085797122667
@@ -671,11 +677,17 @@ end
     @test coef(model) ≈ [4.733877229152363, -1.007977895471349, 0.02522392818548873,
                          -0.9859743168046422, 0.2132095063819721, 0.7456070470961186,
                          0.5840284357554036] rtol = 1e-07
-    @test stderror(model) ≈ [0.42307979153860564, 0.286636744566765, 0.29612422536777805,
-                             0.42042723748229144, 0.45565954626859695, 0.4766324296069839,
-                             0.3235019638755972] rtol = 1e-06
-    @test_broken aic(model) ≈ 4025.0711662068925 rtol = 1e-07
-    @test_broken bic(model) ≈ 4045.956412558851 rtol = 1e-07
+    @test dof(model) == 7
+    # Values match summary(model, dispersion=1) in R with
+    # model = glm(..., family=negative.binomial(2))
+    # but by default summary.glm estimates the dispersion instead of fixing
+    # it to 1 as it should
+    @test stderror(model) ≈ [0.211629897728351, 0.14337934865016283,
+                             0.14812510732683284, 0.21030305642114705,
+                             0.2279267057995862, 0.23841760903991838,
+                             0.16181980065967336] rtol = 1e-06
+    @test aic(model) ≈ 4025.0711662068925 rtol = 1e-07
+    @test bic(model) ≈ 4045.956412558851 rtol = 1e-07
     @test GLM.momentmatrix(model) ≈
           [-1.4294351675636041 -0.0 -1.4294351675636041 -0.0 -0.0 -0.0 -1.4294351675636041
            -1.5410055711037194 -0.0 -1.5410055711037194 -0.0 -0.0 -0.0 -1.5410055711037194