Skip to content

Fix handling of zero weights #638

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
nalimilan wants to merge 3 commits into
base: master
Choose a base branch
from
Open

Fix handling of zero weights #638

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
a0e1815
Fix handling of zero weights
nalimilan Jan 11, 2026
add9660
Review fixes
nalimilan Jan 12, 2026
19e10d4
Merge branch 'master' into nl/zerowts
nalimilan Apr 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 3 additions & 8 deletions src/glmfit.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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(μ)
Expand All @@ -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
Expand Down
17 changes: 15 additions & 2 deletions src/linpred.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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

Expand Down
9 changes: 2 additions & 7 deletions src/lm.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,19 +66,14 @@ 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
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}})
Expand All @@ -105,7 +100,7 @@ function residuals(r::LmResp; weighted::Bool=false)
end
end

link(rr::LmResp) = IdentityLink()
distr(rr::LmResp) = Normal()

"""
LinearModel
Expand Down
46 changes: 46 additions & 0 deletions test/runtests.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]

@nalimilan nalimilan Jan 11, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gragusa I get slightly different standard errors with svyglm. Is this expected? Generally in my tests I got exactly the same values up to a least four decimals.

> X = 1:10
> y = c(1, 4, 6, 2, 3, 5, 6, 7, 1, 6)
> wts = c(1, 0, 4, 0, 5, 2, 6, 4, 2, 6)
> svyd <- svydesign(~1, weights=wts, data=data.frame(X=X, y=y, wts=wts))

> summary(svyglm(y ~ X, svyd))

Call:
svyglm(formula = y ~ X, design = svyd)

Survey design:
svydesign(~1, weights = wts, data = data.frame(X = X, y = y, 
    wts = wts))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   3.6707     1.7371   2.113    0.079 .
X             0.2073     0.2313   0.896    0.405  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 3.500116)

@nalimilan nalimilan Apr 26, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gragusa I'd appreciate your input on this one too.

@gragusa gragusa Apr 27, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry for the late reply - I don't get the notifications from this repo (although I am subscribed).

I think they match to 4 digits (the one used by R to print the coefficient):

## GLM is master/main
using GLM, DataFrames
X = 1:10
y = [1, 4, 6, 2, 3, 5, 6, 7, 1, 6]
wts = [1, 0, 4, 0, 5, 2, 6, 4, 2, 6]
df = DataFrame(y=y, X=X, wts=wts)
julia> lm1 = lm(@formula(y~X), df, weights=GLM.pweights(df.wts))
LinearModel

y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
                Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept)  3.67073     1.73707   2.11    0.0675  -0.334964   7.67643
X            0.207317    0.231276  0.90    0.3962  -0.326007   0.740641
───────────────────────────────────────────────────────────────────────

The variance returned by R on this example is

R> ss$cov.scaled
            (Intercept)           X
(Intercept)   3.0174195 -0.37802382
X            -0.3780238  0.05348872

which matches Julia's GMM:

julia> @show vcov(ll);
vcov(lm1) = [3.0174195202295895 -0.378023819721369; -0.3780238197213689 0.05348872288847028]

@gragusa gragusa Apr 27, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, the problem is that the deviance reported by Julia and R differs.

- `GLM`: deviance = 30.626, dispersion = 3.828
- `svyglm`: deviance = 24.501, dispersion = 3.500

I'm unable to view the code for syvglm for the license. However, the ratio between the two is exactly n_all/n_nz = 10/8. The difference seems to be that syvglm normalizes the weights by the mean over non-zero observations, while GLM normalizes by the mean over all observations.

Let me try to fix this.

@nalimilan nalimilan Apr 27, 2026 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. The "problem" is that this PR changes the result. I get this:

LinearModel

y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
                Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept)  3.67073     1.76171   2.08    0.0823  -0.640024   7.98149
X            0.207317    0.234557  0.88    0.4108  -0.366623   0.781257
───────────────────────────────────────────────────────────────────────

julia> vcov(lm1)
2×2 Matrix{Float64}:
  3.10363   -0.388825
 -0.388825   0.055017

julia> deviance(lm1)
24.50081300813008

This PR seems more correct to me as vcov uses the correct number of observations (skipping zero weights). But R's survey doesn't seem to do the same?

@gragusa gragusa Apr 27, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll try .... 🤞🏽

@gragusa gragusa Apr 27, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PR open: #648

@nalimilan nalimilan Apr 27, 2026 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm so this works but are we sure it's right? :-) Do you think svyglm does not exclude observations with zero weights on purpose?

EDIT: and svyglm prints a warning when fitting a model with zero-weight observations, so it doesn't seem results are supposed to be trusted in that case.

1: In summary.glm(g) :
  observations with zero weight not used for calculating dispersion

@gragusa gragusa Apr 27, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, difficult to say what syvglm does internally. This is GPL code, so I'm not able to review it. Having weights equal to zero is a pathological case (observations with zero weights should not exist). The fix now handles this case by giving the same SE and deviance/dispersion as one would obtain by running the same model, while omitting observations with zero weights, which I think is the right thing to do. Throwing a warning is something we should do, since 0 weights may indicate an issue on the model's side.

@nalimilan nalimilan Apr 27, 2026

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The fix now handles this case by giving the same SE and deviance/dispersion as one would obtain by running the same model, while omitting observations with zero weights, which I think is the right thing to do.

This is was my goal with this PR, but I don't think that's the case with #648 at the moment, right?
Current state of this PR:

julia> lm1 = lm(@formula(y~X), df, weights=GLM.pweights(df.wts))
[ Warning: pweights is defined in StatsBase and is not public in GLM
LinearModel

y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
                Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept)  3.67073     1.76171   2.08    0.0823  -0.640024   7.98149
X            0.207317    0.234557  0.88    0.4108  -0.366623   0.781257
───────────────────────────────────────────────────────────────────────

julia> df2 = df[df.wts .> 0,:];

julia> lm1 = lm(@formula(y~X), df2, weights=GLM.pweights(df2.wts))
LinearModel

y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
                Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept)  3.67073     1.76171   2.08    0.0823  -0.640024   7.98149
X            0.207317    0.234557  0.88    0.4108  -0.366623   0.781257
───────────────────────────────────────────────────────────────────────

Current state of #648:

julia> lm1 = lm(@formula(y~X), df, weights=GLM.pweights(df.wts))
[ Warning: pweights is defined in StatsBase and is not public in GLM
LinearModel

y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
                Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept)  3.67073     1.73707   2.11    0.0790  -0.579731    7.92119
X            0.207317    0.231276  0.90    0.4046  -0.358596    0.77323
───────────────────────────────────────────────────────────────────────

julia> df2 = df[df.wts .> 0,:];

julia> lm1 = lm(@formula(y~X), df2, weights=GLM.pweights(df2.wts))
LinearModel

y ~ 1 + X

Coefficients:
───────────────────────────────────────────────────────────────────────
                Coef.  Std. Error     t  Pr(>|t|)  Lower 95%  Upper 95%
───────────────────────────────────────────────────────────────────────
(Intercept)  3.67073     1.76171   2.08    0.0823  -0.640024   7.98149
X            0.207317    0.234557  0.88    0.4108  -0.366623   0.781257
───────────────────────────────────────────────────────────────────────

If that's useful I can look at the svyglm code for you and/or write to Thomas Lumley (I think he's responsive).

Zero weights should be supported without a warning IMO as they can arise sometimes. For example I have a survey where we have two different definitions of who is in scope, and it's convenient to have two different weighting columns, with zero weights for those who are out of scope in a given definition. We simply need to skip observations with zero weights everywhere and everything should be fine.

end
Loading