# Equivalence of svyglm and glm for simple random surveys

I have been exploring the use of the svyglm function in R's survey package to analyse surveys with both equal and unequal sampling probabilites.

For an unequal sampling probability design (where the sampling units have different probabilities of being included in the sample) I understand that the use of survey weights in svyglm is preferable to a weighted least squares (WLS) regression with weights = 1/sampling probability.

But when sampling probabilities are equal (as with a simple random sampling design), I expected the glm and svyglm models to produce identical results, but this appears not to be the case...

data <- data.frame(y=rnorm(20,10,2),x=rep(c("A","B"),each=10),z=rnorm(20,10,2), wts=rep(1,20))
summary(glm(y~x+z, data=data), weights=wts)

Call:
glm(formula = y ~ x + z, data = data)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.3396  -1.0414   0.0092   0.9199   3.2734

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  9.38740    2.66991   3.516  0.00265 **
xB           0.26102    0.73897   0.353  0.72827
z            0.08557    0.26961   0.317  0.75481
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.497642)

Null deviance: 43.254  on 19  degrees of freedom
Residual deviance: 42.460  on 17  degrees of freedom
AIC: 79.814

Number of Fisher Scoring iterations: 2


whereas svyglm gives the same coefficients but with much smaller standard errors:

des <- svydesign(ids=~1, weights=~wts, data=data)
des

Independent Sampling design (with replacement)
svydesign(ids = ~1, weights = ~wts, data = data)

summary(svyglm(y~x+z, design=des))

Call:
svyglm(formula = y ~ x + z, design = des)

Survey design:
svydesign(ids = ~1, weights = ~wts, data = data)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  9.38740    1.78979   5.245 6.58e-05 ***
xB           0.26102    0.71788   0.364    0.721
z            0.08557    0.18729   0.457    0.654
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 2.234732)

Number of Fisher Scoring iterations: 2


Why are the standard errors so different when the weights are even? Why doesn't the svyglm approach give equivalent results to glm in this situation?

Grateful for any help that can be provided.

You'd get a closer match using model-agnostic 'sandwich' variance estimators. With equal weights, glm assumes that the error variances for all observations are equal, but svyglm makes no such assumption. If the assumption is true, the glm standard errors will be more accurate (the svyglm standard error estimates will be noisier); if the assumption is false, the svyglm standard error estimates will be valid and the glm ones will be biased.