Diseño complejo de encuestas

Question

Actualmente estoy analizando un conjunto de datos resultante de un complejo diseño de encuesta. Los individuos se han seleccionado a partir de un diseño de muestreo por conglomerados en tres etapas con dos estratos (que se han combinado para dar lugar a un estrato grande). Además, se ha estimado un factor de ponderación para garantizar la representatividad y se han producido correcciones de población finita en cada nivel (obtenidas dividiendo el número de conglomerados (o unidades) incluidos por el número potencial total de conglomerados (o unidades)).

Dado que se trata de datos de salud pública y que nuestro artículo se publicará en acceso abierto, me gustaría realizar los análisis en R para que cualquiera pueda reproducirlos. Sin embargo, aunque es bastante sencillo tener en cuenta este diseño en Stata, tuve problemas para hacerlo en R.

En stata, la sintaxis es

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

En R, he intentado pasar a través del paquete "encuesta":

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

En esta fase, todo funciona bien. Sin embargo, cuando se sigue con:

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

\=> Se calculan las estimaciones, pero no el SE ni los valores p.

He probado sin el fpc y todo funcionó bien. Alguien tiene una idea de lo que está mal en mi código (o tiene una idea de otro paquete que se podría utilizar (he comprobado lme4 pero no parece diseñado para ello))?

Muchas gracias por su ayuda.

Answer 1

No hay un problema general con los diseños de tres etapas en svydesign . He aquí un ejemplo inventado en el que todo el análisis funciona, con tres etapas, con correcciones de población finita y con tres enfoques de la ponderación. Primero, ponderaciones simplemente calculadas a partir de las fpcs. Segundo, ponderaciones suministradas, pero iguales a las que se obtendrían de las fpcs. Tercero, pesos suministrados y diferentes de los que se obtendrían de las 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

¿Puede dar más información sobre el resultado obtenido? ¿Hubo alguna advertencia? (Veo que utilizó una opción en Stata para manejar las PSU únicas; ¿hizo lo equivalente en R con options("survey.lonely.psu")?