diff --git a/src/glmfit.jl b/src/glmfit.jl index 21e2aa6f..9e4b9307 100644 --- a/src/glmfit.jl +++ b/src/glmfit.jl @@ -25,7 +25,7 @@ struct GlmResp{V<:FPVector,D<:UnivariateDistribution,L<:Link,W<:AbstractWeights} wrkresid::V end -link(rr::GlmResp) = rr.d +distr(rr::GlmResp) = rr.d function GlmResp(y::V, d::D, l::L, η::V, μ::V, off::V, weights::W) where {V<:FPVector,D,L,W} n = length(y) @@ -343,7 +343,7 @@ function loglikelihood(m::AbstractGLM) y = r.y mu = r.mu wts = weights(r) - d = link(r) + d = distr(r) ll = zero(eltype(mu)) N = length(y) δ = deviance(r) @@ -860,11 +860,6 @@ function checky(y, d::Binomial) return nothing end -function nobs(r::GlmResp{V,D,L,W}) where {V,D,L,W<:AbstractWeights} - return oftype(sum(one(eltype(weights(r)))), length(r.y)) -end -nobs(r::GlmResp{V,D,L,W}) where {V,D,L,W<:FrequencyWeights} = sum(r.weights) - function residuals(r::GlmResp; weighted::Bool=false) y, η, μ = r.y, r.eta, r.mu dres = similar(μ) @@ -885,7 +880,7 @@ end function momentmatrix(m::GeneralizedLinearModel) X = modelmatrix(m; weighted=false) r = varstruct(m) - if link(m) isa Union{Gamma,InverseGaussian} + if distr(m) isa Union{Gamma,InverseGaussian} r .*= sum(working_weights(m)) / sum(abs2, r) end return Diagonal(r) * X diff --git a/src/linpred.jl b/src/linpred.jl index 50a51375..433446ac 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -317,7 +317,7 @@ function vcov(x::LinPredModel) s = nobs(x) / (nobs(x) - 1) mm = momentmatrix(x) A = invloglikhessian(x) - if link(x) isa Union{Gamma,InverseGaussian} + if distr(x) isa Union{Gamma,InverseGaussian} r = varstruct(x) A ./= sum(working_weights(x)) / sum(abs2, r) end @@ -327,7 +327,7 @@ function vcov(x::LinPredModel) end end -link(x::LinPredModel) = link(x.rr) +distr(x::LinPredModel) = distr(x.rr) function _vcov(pp::LinPred, Z::AbstractMatrix, A::AbstractMatrix) if linpred_rank(pp) < size(Z, 2) @@ -450,9 +450,22 @@ end For linear and generalized linear models, return the number of rows when the model is unweighted or uses analytical or probability weights. If the model uses frequency weights, return the sum of weights. +Rows with zero weights are not counted. """ nobs(obj::LinPredModel) = nobs(obj.rr) +function nobs(rr::ModResp) + wts = weights(rr) + n = if wts isa UnitWeights + length(wts) + elseif wts isa FrequencyWeights + sum(wts) + elseif wts isa Union{ProbabilityWeights,AnalyticWeights} + count(!iszero, wts) + end + return oftype(sum(one(eltype(wts))), n) +end + weights(m::LinPredModel) = weights(m.rr) weights(pp::LinPred) = pp.weights diff --git a/src/lm.jl b/src/lm.jl index cb74e976..6000f78e 100644 --- a/src/lm.jl +++ b/src/lm.jl @@ -66,7 +66,7 @@ function deviance(r::LmResp) v += abs2(y[i] - mu[i]) * weights[i] end end - return weights isa ProbabilityWeights ? v ./ (sum(weights) / length(y)) : v + return weights isa ProbabilityWeights ? v ./ (sum(weights) / nobs(r)) : v end weights(r::LmResp) = r.weights @@ -74,11 +74,6 @@ function isweighted(r::LmResp) return weights(r) isa Union{AnalyticWeights,FrequencyWeights,ProbabilityWeights} end -nobs(r::LmResp{<:Any,<:FrequencyWeights}) = sum(r.weights) -function nobs(r::LmResp{<:Any,<:AbstractWeights}) - return oftype(sum(one(eltype(r.weights))), length(r.y)) -end - working_weights(r::LmResp) = r.weights function loglikelihood(r::LmResp{<:Any,<:Union{UnitWeights,FrequencyWeights}}) @@ -105,7 +100,7 @@ function residuals(r::LmResp; weighted::Bool=false) end end -link(rr::LmResp) = IdentityLink() +distr(rr::LmResp) = Normal() """ LinearModel diff --git a/test/runtests.jl b/test/runtests.jl index c0b8c62b..97277384 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2701,3 +2701,49 @@ end x -0.502451 0.675377 -0.74 0.4982 -2.3776 1.3727 ─────────────────────────────────────────────────────────────────────────""" end + +@testset "nobs, deviance & stderror with zero weights" begin + X = [ones(10) 1:10] + y = [1, 4, 6, 2, 3, 5, 6, 7, 1, 6] + wts = [1, 0, 4, 0, 5, 2, 6, 4, 2, 6] + + lmod = lm(X, y) + glmod = glm(X, y, Normal()) + @test nobs(lmod) == nobs(glmod) == 10 + @test dof_residual(lmod) == dof_residual(glmod) == 8 + @test deviance(lmod) ≈ deviance(glmod) ≈ 39.29696969696969 + @test coef(lmod) ≈ coef(glmod) ≈ [2.6666666666666665, 0.2606060606060606] + @test stderror(lmod) ≈ stderror(glmod) ≈ [1.51404201801774, 0.24400996532360475] + + lmod = lm(X, y, weights=uweights(10)) + glmod = glm(X, y, Normal(), weights=uweights(10)) + @test nobs(lmod) == nobs(glmod) == 10 + @test dof_residual(lmod) == dof_residual(glmod) == 8 + @test deviance(lmod) ≈ deviance(glmod) ≈ 39.29696969696969 + @test coef(lmod) ≈ coef(glmod) ≈ [2.6666666666666665, 0.2606060606060606] + @test stderror(lmod) ≈ stderror(glmod) ≈ [1.51404201801774, 0.24400996532360475] + + lmod = lm(X, y, weights=fweights(wts)) + glmod = glm(X, y, Normal(), weights=fweights(wts)) + @test nobs(lmod) == nobs(glmod) == sum(wts) + @test dof_residual(lmod) == dof_residual(glmod) == sum(wts) - 2 + @test deviance(lmod) ≈ deviance(glmod) ≈ 91.87804878048782 + @test coef(lmod) ≈ coef(glmod) ≈ [3.6707317073170724, 0.20731707317073195] + @test stderror(lmod) ≈ stderror(glmod) ≈ [0.9538283539174621, 0.13286974786560674] + + lmod = lm(X, y, weights=aweights(wts)) + glmod = glm(X, y, Normal(), weights=aweights(wts)) + @test nobs(lmod) == nobs(glmod) == count(!iszero, wts) + @test dof_residual(lmod) == dof_residual(glmod) == count(!iszero, wts) - 2 + @test deviance(lmod) ≈ deviance(glmod) ≈ 91.87804878048782 + @test coef(lmod) ≈ coef(glmod) ≈ [3.6707317073170724, 0.20731707317073195] + @test stderror(lmod) ≈ stderror(glmod) ≈ [2.060504744176091, 0.2870314608599428] + + lmod = lm(X, y, weights=pweights(wts)) + glmod = glm(X, y, Normal(), weights=pweights(wts)) + @test nobs(lmod) == nobs(glmod) == count(!iszero, wts) + @test dof_residual(lmod) == dof_residual(glmod) == count(!iszero, wts) - 2 + @test deviance(lmod) ≈ deviance(glmod) ≈ 24.500813008130084 + @test coef(lmod) ≈ coef(glmod) ≈ [3.6707317073170724, 0.20731707317073195] + @test stderror(lmod) ≈ stderror(glmod) ≈ [1.7617126628715203, 0.23455696986842048] +end