Fix handling of zero weights by nalimilan · Pull Request #638 · JuliaStats/GLM.jl

nalimilan · 2026-01-11T22:47:22Z

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 LinPredModels (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.

Improves on #487. Cc: @gragusa

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.

nalimilan · 2026-01-11T22:51:40Z

+    @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]


@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)

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

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]

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.

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?

I'll try .... 🤞🏽

PR open: #648

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

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.

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.

codecov · 2026-01-11T22:55:42Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 97.21%. Comparing base (e1438c5) to head (19e10d4).

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #638      +/-   ##
==========================================
+ Coverage   97.20%   97.21%   +0.01%     
==========================================
  Files           8        8              
  Lines        1216     1222       +6     
==========================================
+ Hits         1182     1188       +6     
  Misses         34       34

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

andreasnoack · 2026-04-24T21:21:25Z

There are some conflicts here that need to be resolved

nalimilan · 2026-04-26T20:00:12Z

Done.

gragusa · 2026-04-27T09:53:43Z

I have a PR that fixes these cases against the nl/zerowts branch. What should I do? I do not think I can push it.

nalimilan commented Jan 11, 2026

View reviewed changes

nalimilan added the bug label Jan 11, 2026

nalimilan added this to the Release 2.0 milestone Jan 11, 2026

nalimilan requested review from andreasnoack and devmotion January 11, 2026 22:55

devmotion reviewed Jan 11, 2026

View reviewed changes

Comment thread src/glmfit.jl Outdated

Comment thread src/glmfit.jl Outdated

Comment thread src/linpred.jl Outdated

Comment thread src/lm.jl Outdated

Review fixes

add9660

Merge branch 'master' into nl/zerowts

19e10d4

gragusa mentioned this pull request Apr 27, 2026

Improve documentation #632

Open

Uh oh!

Conversation

nalimilan commented Jan 11, 2026

Uh oh!

nalimilan Jan 11, 2026

Choose a reason for hiding this comment

Uh oh!

nalimilan Apr 26, 2026

Choose a reason for hiding this comment

Uh oh!

gragusa Apr 27, 2026

Choose a reason for hiding this comment

Uh oh!

gragusa Apr 27, 2026

Choose a reason for hiding this comment

Uh oh!

nalimilan Apr 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gragusa Apr 27, 2026

Choose a reason for hiding this comment

Uh oh!

gragusa Apr 27, 2026

Choose a reason for hiding this comment

Uh oh!

nalimilan Apr 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

gragusa Apr 27, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nalimilan Apr 27, 2026

Choose a reason for hiding this comment

Uh oh!

codecov Bot commented Jan 11, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

andreasnoack commented Apr 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

nalimilan commented Apr 26, 2026

Uh oh!

gragusa commented Apr 27, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

nalimilan Apr 27, 2026 •

edited

Loading

nalimilan Apr 27, 2026 •

edited

Loading

gragusa Apr 27, 2026 •

edited

Loading

codecov Bot commented Jan 11, 2026 •

edited

Loading

andreasnoack commented Apr 24, 2026 •

edited

Loading