`dof()` incorrect for `NegativeBinomial()` and `Geometric()`

[gragusa]

The dof() method for GeneralizedLinearModel incorrectly counts the dispersion parameter θ for NegativeBinomial(θ) and Geometric() distributions, even though θ is fixed (not estimated). This causes dof() to be off by 1 (too high) and the AIC/BIC to be inflated.

For distributions where dispersion IS estimated (Normal, Gamma, InverseGaussian), Julia's dof() correctly matches R, and AIC agrees. The bug only affects NegativeBinomial and Geometric.

Here a comparison between GLM and R (using MASS package)

Distribution	dof (Jl)	AIC (Jl)	loglik_df (R)	AIC (R)
Geometric	8	1112.74	7	1110.74
NegBinomial(2)	8	1122.52	7	1120.52

The cause is in src/glmfit.jl

function dof(x::GeneralizedLinearModel)
    modelrank = linpred_rank(x.pp)
    dispersion_parameter(x.rr.d) ? modelrank + 1 : modelrank
end

The issue is that dispersion_parameter() returns true for NegativeBinomial and Geometric, but their θ parameter is fixed by the user, not estimated from the data. The +1 should only apply when dispersion is actually being estimated.

Since #485 is almost ready, I propose tackling this issue after it is merged. In #485, dispersion_parameter() is used in the calculation of the dispersion parameter for weighted models, and the best option is a special handling in dof() to not count fixed θ for these distributions.

[andreasnoack]

I believe this is mostly a duplicate of #564

[nalimilan]

See #624.

[andreasnoack]

As mentioned in #624, I think the two cases here are separate. Geometric is just wrong, but I don't think our negative binomial is. The dof value should just be the number of estimated parameters. That also seems to be the definition of MASS::glm.nb. I suppose R values in the table above are based on a fixed, and not an estimated, shape paramter?

> ft.fixed <- glm(Days ~ Eth + Sex + Age + Lrn, quine, family=negative.binomial(2, link="log"))
> ft.estimated <- glm.nb(Days ~ Eth + Sex + Age + Lrn, quine, link="log")
> logLik(ft.fixed)
'log Lik.' -553.2596 (df=7)
> logLik(ft.estimated)
'log Lik.' -546.5755 (df=8)
> ft.fixed$aic
[1] 1120.519
> ft.estimated$aic
[1] 1109.151

Update: Never mind. I should have realized that you refer to the negative binomial with a fixed shape in GLM.jl.