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!