From a0e181559baab2b6f5950b5041846bfa3d147c12 Mon Sep 17 00:00:00 2001 From: Milan Bouchet-Valat Date: Sun, 11 Jan 2026 23:31:41 +0100 Subject: [PATCH 1/2] Fix handling of zero weights Observations with zero weights have to be skipped for probability and analytic weights. Values have been checked against R. Avoid duplication by using a single method for all `LinPredModel`s (though this means we cannot call `nobs` from functions which don't get passed the full model). Also fix name of internal function `link`, which actually returns the distribution. --- src/glmfit.jl | 11 ++++++----- src/linpred.jl | 17 ++++++++++++++--- src/lm.jl | 11 +++-------- test/runtests.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 69 insertions(+), 16 deletions(-) diff --git a/src/glmfit.jl b/src/glmfit.jl index 0296af2e..ab88324a 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, wts::W) where {V<:FPVector,D,L,W} n = length(y) @@ -79,7 +79,7 @@ end function deviance(r::GlmResp) wts = weights(r) d = sum(r.devresid) - return wts isa ProbabilityWeights ? d * nobs(r) / sum(wts) : d + return wts isa ProbabilityWeights ? d * sum(!iszero, wts) / sum(wts) : d end weights(r::GlmResp) = r.wts @@ -342,7 +342,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,7 +860,8 @@ function checky(y, d::Binomial) 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)) + n = W<:ProbabilityWeights ? sum(!iszero, weights(r)) : length(r.y) + return oftype(sum(one(eltype(weights(r)))), n) end nobs(r::GlmResp{V,D,L,W}) where {V,D,L,W<:FrequencyWeights} = sum(r.wts) @@ -884,7 +885,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 b229e355..44a4ece4 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -316,7 +316,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 @@ -326,7 +326,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) @@ -434,8 +434,19 @@ 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(m::LinPredModel) + wts = weights(m) + n = if wts isa UnitWeights + length(wts) + elseif wts isa FrequencyWeights + sum(wts) + elseif wts isa Union{ProbabilityWeights,AnalyticWeights} + sum(!iszero, wts) + end + return oftype(sum(one(eltype(wts))), n) +end weights(m::LinPredModel) = weights(m.rr) weights(pp::LinPred) = pp.wts diff --git a/src/lm.jl b/src/lm.jl index 48d4034d..b94e32b0 100644 --- a/src/lm.jl +++ b/src/lm.jl @@ -67,7 +67,7 @@ function deviance(r::LmResp) v += abs2(y[i] - mu[i]) * wts[i] end end - return wts isa ProbabilityWeights ? v ./ (sum(wts) / length(y)) : v + return wts isa ProbabilityWeights ? v ./ (sum(wts) / sum(!iszero, wts)) : v end weights(r::LmResp) = r.wts @@ -75,15 +75,10 @@ function isweighted(r::LmResp) return weights(r) isa Union{AnalyticWeights,FrequencyWeights,ProbabilityWeights} end -nobs(r::LmResp{<:Any,<:FrequencyWeights}) = sum(r.wts) -function nobs(r::LmResp{<:Any,<:AbstractWeights}) - return oftype(sum(one(eltype(r.wts))), length(r.y)) -end - working_weights(r::LmResp) = r.wts function loglikelihood(r::LmResp{<:Any,<:Union{UnitWeights,FrequencyWeights}}) - n = nobs(r) + n = sum(weights(r)) return -n / 2 * (log(2π * deviance(r) / n) + 1) end @@ -106,7 +101,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 e699039d..ad9b15fa 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2659,3 +2659,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, wts=uweights(10)) + glmod = glm(X, y, Normal(), wts=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, wts=fweights(wts)) + glmod = glm(X, y, Normal(), wts=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, wts=aweights(wts)) + glmod = glm(X, y, Normal(), wts=aweights(wts)) + @test nobs(lmod) == nobs(glmod) == sum(!iszero, wts) + @test dof_residual(lmod) == dof_residual(glmod) == sum(!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, wts=pweights(wts)) + glmod = glm(X, y, Normal(), wts=pweights(wts)) + @test nobs(lmod) == nobs(glmod) == sum(!iszero, wts) + @test dof_residual(lmod) == dof_residual(glmod) == sum(!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 From add96604815f0822d622fcc736111008d5bd69b4 Mon Sep 17 00:00:00 2001 From: Milan Bouchet-Valat Date: Mon, 12 Jan 2026 17:25:11 +0100 Subject: [PATCH 2/2] Review fixes --- src/glmfit.jl | 8 +------- src/linpred.jl | 8 +++++--- src/lm.jl | 4 ++-- test/runtests.jl | 8 ++++---- 4 files changed, 12 insertions(+), 16 deletions(-) diff --git a/src/glmfit.jl b/src/glmfit.jl index ab88324a..8c8fd143 100644 --- a/src/glmfit.jl +++ b/src/glmfit.jl @@ -79,7 +79,7 @@ end function deviance(r::GlmResp) wts = weights(r) d = sum(r.devresid) - return wts isa ProbabilityWeights ? d * sum(!iszero, wts) / sum(wts) : d + return wts isa ProbabilityWeights ? d * nobs(r) / sum(wts) : d end weights(r::GlmResp) = r.wts @@ -859,12 +859,6 @@ function checky(y, d::Binomial) return nothing end -function nobs(r::GlmResp{V,D,L,W}) where {V,D,L,W<:AbstractWeights} - n = W<:ProbabilityWeights ? sum(!iszero, weights(r)) : length(r.y) - return oftype(sum(one(eltype(weights(r)))), n) -end -nobs(r::GlmResp{V,D,L,W}) where {V,D,L,W<:FrequencyWeights} = sum(r.wts) - function residuals(r::GlmResp; weighted::Bool=false) y, η, μ = r.y, r.eta, r.mu dres = similar(μ) diff --git a/src/linpred.jl b/src/linpred.jl index 44a4ece4..0e21e314 100644 --- a/src/linpred.jl +++ b/src/linpred.jl @@ -436,14 +436,16 @@ 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. """ -function nobs(m::LinPredModel) - wts = weights(m) +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} - sum(!iszero, wts) + count(!iszero, wts) end return oftype(sum(one(eltype(wts))), n) end diff --git a/src/lm.jl b/src/lm.jl index b94e32b0..5a3a794a 100644 --- a/src/lm.jl +++ b/src/lm.jl @@ -67,7 +67,7 @@ function deviance(r::LmResp) v += abs2(y[i] - mu[i]) * wts[i] end end - return wts isa ProbabilityWeights ? v ./ (sum(wts) / sum(!iszero, wts)) : v + return wts isa ProbabilityWeights ? v ./ (sum(wts) / nobs(r)) : v end weights(r::LmResp) = r.wts @@ -78,7 +78,7 @@ end working_weights(r::LmResp) = r.wts function loglikelihood(r::LmResp{<:Any,<:Union{UnitWeights,FrequencyWeights}}) - n = sum(weights(r)) + n = nobs(r) return -n / 2 * (log(2π * deviance(r) / n) + 1) end diff --git a/test/runtests.jl b/test/runtests.jl index ad9b15fa..0ab8820c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2691,16 +2691,16 @@ end lmod = lm(X, y, wts=aweights(wts)) glmod = glm(X, y, Normal(), wts=aweights(wts)) - @test nobs(lmod) == nobs(glmod) == sum(!iszero, wts) - @test dof_residual(lmod) == dof_residual(glmod) == sum(!iszero, wts) - 2 + @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, wts=pweights(wts)) glmod = glm(X, y, Normal(), wts=pweights(wts)) - @test nobs(lmod) == nobs(glmod) == sum(!iszero, wts) - @test dof_residual(lmod) == dof_residual(glmod) == sum(!iszero, wts) - 2 + @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]