Match svyglm for ProbabilityWeights deviance, nulldeviance, and SE

Author: gragusa

src/glmfit.jl

@@ -80,7 +80,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 * length(r.y) / sum(wts) : d
 end
 
 weights(r::GlmResp) = r.weights
@@ -323,7 +323,7 @@ function nulldeviance(m::GeneralizedLinearModel)
             end
         end
         if wts isa ProbabilityWeights
-            dev /= sum(wts) / nobs(m)
+            dev /= sum(wts) / length(y)
         end
     else
         X = fill(1.0, length(y), hasint ? 1 : 0)

src/linpred.jl

@@ -313,8 +313,10 @@ working_weights(x::LinPredModel) = working_weights(x.rr)
 function vcov(x::LinPredModel)
     if weights(x) isa ProbabilityWeights
         ## n-1 degrees of freedom - This is coherent with the `R` package `survey`,
-        ## `STATA` uses n-k
-        s = nobs(x) / (nobs(x) - 1)
+        ## `STATA` uses n-k. Use the full sample size (including zero-weighted obs)
+        ## so the correction matches `survey::svyglm`.
+        n = length(response(x))
+        s = n / (n - 1)
         mm = momentmatrix(x)
         A = invloglikhessian(x)
         if distr(x) isa Union{Gamma,InverseGaussian}

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) / nobs(r)) : v
+    return weights isa ProbabilityWeights ? v ./ (sum(weights) / length(r.y)) : v
 end
 
 weights(r::LmResp) = r.weights
@@ -235,7 +235,7 @@ function nulldeviance(obj::LinearModel)
             v += abs2(y[i] - m) * weights[i]
         end
     end
-    return v
+    return weights isa ProbabilityWeights ? v / (sum(weights) / length(y)) : v
 end
 
 function nullloglikelihood(m::LinearModel)

test/runtests.jl

@@ -2743,7 +2743,8 @@ end
     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 deviance(lmod) ≈ deviance(glmod) ≈ 30.626016260162601
+    @test nulldeviance(lmod) ≈ nulldeviance(glmod) ≈ 33.288888888888884
     @test coef(lmod) ≈ coef(glmod) ≈ [3.6707317073170724, 0.20731707317073195]
-    @test stderror(lmod) ≈ stderror(glmod) ≈ [1.7617126628715203, 0.23455696986842048]
+    @test stderror(lmod) ≈ stderror(glmod) ≈ [1.7370721114074674, 0.2312762912372783]
 end