O Thomas Lumley já escreveu sobre isso; link aqui!

Uma das coisas mais interessantes que pesquisas domiciliares permitem fazer é filtrar domínios de análise. Por exemplo, além da renda média da população, é possível estimar as rendas médias das áreas urbanas e áreas rurais. Para avaliar a precisão dessas estimativas, estimamos a variância do estimador e o intervalo de confiança.

No entanto, Tem uma prática problemática em alguns exemplos deste tipo de análise com o pacote survey: deletar as linhas antes de descrever o plano amostral. Sempre que eu comento filtrar observações desta maneira é arriscado, ouço o algo nas linhas de “minha análise é sobre populações indígenas, deletar informações de pessoas de outra cor/raça não influencia na análise!”.

O problema é que elas influenciam.

O caso da POF

Vamos considerar o seguinte problema: tomando a população de referência da POF, qual é o tamanho da população de indígenas em São Paulo?

Deletando linhas na crianção do objeto de plano amostral:

# cria variável unitária
pof.df$one <- 1

# cria objeto de desenho amostral só com as linhas de pessoas indígenas em São Paulo
pof.des.rows <- svydesign( ids = ~cod_upa , strata = ~estrato_pof , data = pof.df[ uf == 35 & v0405 == 5 ] , weights = ~peso_final )

# # opção "conservadora" para tratar estratos unitários
options( survey.lonely.psu="adjust" )

# estima totais e variâncias
svytotal( ~one , pof.des.rows )
##      total    SE
## one 148046 41047

Usando a amostra completa e a função subset:

# cria objeto de desenho amostral com amostra completa
pof.des.survobj <- svydesign( ids = ~cod_upa , strata = ~estrato_pof , data = pof.df , weights = ~peso_final )

# filtra domínio de interesse: pessoas indígenas em São Paulo
pof.des.survobj <- subset( pof.des.survobj , uf == 35 & v0405 == 5 )

# estima totais e variâncias
svytotal( ~one , pof.des.survobj )
##      total    SE
## one 148046 44622

Como podemos notar, as estimativas dos totais são idênticas, mas as variâncias são diferentes. Mas temos o mesmo estimador, (aparentemente) a mesma amostra… Os valores deveriam ser idênticos. Então, fazemos a pergunta:

Qual é a forma correta?

Resposta: a segunda, que considera o plano amostral com a amostra completa.

Por quê?

Compare as contagem de PSUs1 nas saídas abaixo.

summary( pof.des.rows )
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (28) clusters.
## svydesign(ids = ~cod_upa, strata = ~estrato_pof, data = pof.df[uf == 
##     35 & v0405 == 5], weights = ~peso_final)
## Probabilities:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0001323 0.0002264 0.0002532 0.0006395 0.0003734 0.0080181 
## Stratum Sizes: 
##            3501 3502 3505 3508 3510 3517 3522 3523 3526 3530 3531 3535 3538
## obs           2    1   12    2    6    1    1    1    2    1    3    2    1
## design.PSU    2    1    7    2    1    1    1    1    2    1    1    2    1
## actual.PSU    2    1    7    2    1    1    1    1    2    1    1    2    1
##            3540 3541 3542 3547 3551
## obs           1    1    2    2    1
## design.PSU    1    1    1    1    1
## actual.PSU    1    1    1    1    1
## Data variables:
##  [1] "cod_upa"     "uf"          "estrato_pof" "num_dom"     "num_uc"     
##  [6] "v0306"       "peso"        "peso_final"  "v0405"       "pos_estrato"
## [11] "npes"        "one"
summary( pof.des.survobj )
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (28) clusters.
## subset(pof.des.survobj, uf == 35 & v0405 == 5)
## Probabilities:
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0001323 0.0002264 0.0002532 0.0006395 0.0003734 0.0080181 
## Stratum Sizes: 
##            3501 3502 3505 3508 3510 3517 3522 3523 3526 3530 3531 3535 3538
## obs           2    1   12    2    6    1    1    1    2    1    3    2    1
## design.PSU   23   13   19   12   25    7   20   10   10    3   14   13    3
## actual.PSU    2    1    7    2    1    1    1    1    2    1    1    2    1
##            3540 3541 3542 3547 3551
## obs           1    1    2    2    1
## design.PSU    3    3    3    3    5
## actual.PSU    1    1    1    1    1
## Data variables:
##  [1] "cod_upa"     "uf"          "estrato_pof" "num_dom"     "num_uc"     
##  [6] "v0306"       "peso"        "peso_final"  "v0405"       "pos_estrato"
## [11] "npes"        "one"

Na primeira, a contagem de UPAs do plano amostral coincide com a contagem de UPAs no domínio. Já na segunda, estas contagens divergem. O que isso significa? Quando vamos estimar a variância em cada estrato, a primeira considera apenas as UPAs onde há pelo menos uma pessoa indígenas na amostra. Já na segunda, todas as UPAs são consideradas. Ou seja: não é porque uma pessoa indígena não foi selecionada naquela UPA que não existem pessoas indígenas nela. Se fizermos isso, estamos subestimando a variabilidade dessa característica entre as UPAs. De outra maneira, ignorar as outras UPAs seria similar a dizer que \(\sum ( x - \bar{x} )^2 = \sum_{x>0} ( x - \bar{x} )^2\), o que não faz o menor sentido.

Assim, como regra2, o melhor é deixar que a função subset faça o trabalho.

Como isso afeta a convey

No caso da convey, o modo correto de declarar o plano amostral é rodando todo e qualquer subset() depois da convey_prep(). O motivo é o mesmo: não fazer isso declara incorretamente o plano amostral. Assim, o certo é fazer:

# cria objeto de desenho amostral com amostra completa
pof.des.survobj <- svydesign( ids = ~cod_upa , strata = ~estrato_pof , data = pof.df , weights = ~peso_final )

# cria objeto convey
pof.des.survobj <- convey_prep( pof.des.survobj )

# filtra domínio de interesse: pessoas indígenas em São Paulo
pof.des.survobj <- subset( pof.des.survobj , uf == 35 & v0405 == 5 )

  1. Unidades Primárias de Amostragem (UPA) em português.↩︎

  2. Existem casos onde isso não é um problema.↩︎