Complex survey design

Question

I am currently analyzing a dataset resulting from a complex survey design. Individuals have been selected from a three-stage cluster sampling design with two strates (that have been combined to result in one large strate). Moreover, a weighting factor was estimated to ensure representativeness and finite population corrections were produced at each level (obtained by dividing the number of clusters (or units) included by the total potential number of clusters (or units)).

Since this is public health data and that our paper will be published in open access, I would like to perform the analyses in R so that anyone can reproduce them. However, while it is quite straightforward to account for this design in Stata, I had trouble to do so in R.

In stata, the syntax is:

svyset level3cluster [pweight=WEIGHTS], strata(stratum) fpc(fpc3) 
vce(linearized) singleunit(scaled) || level2cluster, fpc(fpc2) || level1unit, fpc(fpc1) 

In R, I tried to pass throught the "survey" package:

Data <- svydesign(ids=~level3cluster+level2cluster+level1unit,
                  strata = ~stratum, weights= ~WEIGHTS,
                  fpc = ~fpc3+fpc2+fpc1, data=data_tot)

At this stage, everything works well. However, when following with:

summary(svyglm(VD~VI, design=Data, family=gaussian))

=> Estimates are calculated but SE and p-values are not.

I have tried without the fpc and everything worked well. Someone has an idea of what is wrong in my code (or has an idea of another package that could be used (I have checked lme4 but it does not seem designed for this))?

Thank you very much for your help!

Answer 1

There isn't a general problem with three-stage designs in svydesign. Here's a made-up example where the analysis all works, with three stages, with finite population corrections, and with three approaches to weighting. First, weights just calculated from the fpcs. Second, weights supplied, but equal to those you'd get from the fpcs. Third, weights supplied and different from those you'd get from the fpcs.

> library(survey)
> 
> set.seed(2020-6-1)
> one<-rep(1:10,each=100)
> two<-rep(1:50,each=20)
> three<-1:1000
> y<-rnorm(1000)
> x<-rnorm(1000)
> 
> d<-data.frame(one=one,two=two,three=three,fpc1=10/420,fpc2=5/27,fpc3=20/123,y=y,x=x)
> 
> des<-svydesign(id=~one+two+three, fpc=~fpc1+fpc2+fpc3,data=d)
> 
> summary(svyglm(y~x,design=des))

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

Survey design:
svydesign(id = ~one + two + three, fpc = ~fpc1 + fpc2 + fpc3, 
    data = d)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.02510    0.04030  -0.623  0.55075   
x            0.10811    0.03117   3.468  0.00847 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.054681)

Number of Fisher Scoring iterations: 2

> 
> d$weights<-with(d,(1/fpc1)*(1/fpc2)*(1/fpc3))
> 
> des2<-svydesign(id=~one+two+three, fpc=~fpc1+fpc2+fpc3,data=d,weights=~weights)
> 
> summary(svyglm(y~x,design=des2))

Call:
svyglm(formula = y ~ x, design = des2)

Survey design:
svydesign(id = ~one + two + three, fpc = ~fpc1 + fpc2 + fpc3, 
    data = d, weights = ~weights)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.02510    0.04030  -0.623  0.55075   
x            0.10811    0.03117   3.468  0.00847 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.054681)

Number of Fisher Scoring iterations: 2

> 
> d$weights3<-with(d, (1/fpc1)*(1/fpc2)*(1/fpc3)*(5+rexp(1000)))
> 
> des3<-svydesign(id=~one+two+three, fpc=~fpc1+fpc2+fpc3,data=d,weights=~weights3)
> 
> summary(svyglm(y~x,design=des3))

Call:
svyglm(formula = y ~ x, design = des3)

Survey design:
svydesign(id = ~one + two + three, fpc = ~fpc1 + fpc2 + fpc3, 
    data = d, weights = ~weights3)

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.01823    0.03893  -0.468  0.65207   
x            0.11031    0.03110   3.547  0.00754 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 1.059041)

Number of Fisher Scoring iterations: 2

Can you give more information on the output you got? Were there any warnings? (I see you used an option in Stata to handle singleton PSUs; did you do the equivalent in R with options("survey.lonely.psu")?