Atualização (25/07/2022): eu e prof. Pedro Silva (ENCE/IBGE) fizemos uma apresentação no IPEA sobre Método dos Conglomerados Primários e estimação de variância do estimador de calibração usando a PNADC. Os slides estão disponíveis neste link.

Recentemente, a PNAD Contínua passou por duas mudanças:

  1. Mudança de pós-estratificação para raking;
  2. Mudança no método de estimação da variância dos estimadores, do Método do Conglomerado Primário para bootstrap.

Mas o que isso significa? E como considerar isso na estimação usando o R? Bom, esse é o tema desse post. Primeiro, vou explicar muito superficialmente o que é raking e porque a mudança no método de estimar a variância.

Calibração e raking

Calibração (Deville e Särndal, 1992; Deville, Särndal e Sautory, 1993) é o nome de uma classe de métodos que, por conseguinte, definem uma classe de estimadores. Estes estimadores incorporam informações adicionais sobre a população. Basicamente, esta classe de métodos procura estabelecer um sistema de ponderação tal que, para uma amostra \(S\) extraída de uma população finita \(U\) usando um plano amostral \(p(S)\), tenhamos

\[ \sum_{k \in S } w_k \mathbf{x}_k = \sum_{ k \in U } \mathbf{x}_k \]

Desta forma, os totais das variáveis \(\mathbf{x}\) são exatamente não-viesados e a variância é (praticamente) nula.

Existe uma infinidade de maneiras de estabelecer pesos \(w_k\) que satisfazem a equação acima. No entanto, é importante lembrar que a nossa amostra \(S\) já tem pesos básicos dados pelo inverso da probabilidade de inclusão – i.e., \(d_k = 1/ \pi_k\). Assim, para minimizar o risco de viés, devemos procurar pesos ajustados \(w_k\) que estejam o mais próximo possível dos pesos originais \(d_k\).

A proximidade entre os vetores \(w_k\) e \(d_k\) é dada por uma pseudo-distância1 \(G_k ( w_k , d_k )\). O problema em questão é um problema de otimização restrita: encontre os pesos \(w_k, (k \in S)\) que minimizem \(\sum_{k \in S} G_k( w_k , d_k ) / q_k\) sujeito à restrição \(\sum_{k \in S } w_k \mathbf{x}_k = \sum_{ k \in U } \mathbf{x}_k\).

Quando resolvemos este problema, temos uma amostra que: (1) possui informações coerentes com as informações externas nas variáveis calibradas, (2) reduz a variância do estimador de expansão, e (3) pode reduzir o erro de não-resposta e de cobertura.

Em geral, conhecemos poucos totais populacionais. Vamos pensar em uma tabela com o seguinte formato: nas linhas, temos faixas etárias, enquanto nas colunas temos o sexo. Quando conhecemos os valores dos totais na população para cada célula da tabela, podemos “corrigir” os pesos da amostra usando uma razão entre os totais na população e os totais estimados na amostra. Esta técnica é chamada de pós-estratificação: de certa forma, estamos criando (ou reajustando) estratos após a seleção da amostra.

Porém, em alguns casos, só conhecemos os totais nas margens desta tabela. O raking funciona neste caso: ajustamos iterativamente as marginais estimadas pela amostra às marginais da população. O problema é que este processo pode levar a pesos extremos, até negativos. É possível estabelecer restrições aos fatores de ajuste para evitar variações extremas. Além disso, em pesquisas domiciliares, costumamos trabalhar com pesos que não variam entre moradores do mesmo domicílio, restrição que também pode ser adicionada ao método.

Mesmo quando conhecemos muitas variáveis na população, devemos escolher variáveis associadas às variáveis de interesse. Quando as variáveis de calibração estão pouco relacionadas às variáveis de interesse, o método apenas aumenta a variabilidade dos pesos, aumentando a variância dos estimadores. As estimativas, no entanto, permanecem (aproximadamente) não-viesadas.

Até aqui, vimos que o processo de calibração gera novos pesos. Estes pesos reajustados são estimadores de calibração, cuja fórmula para estimação de variância é diferente daquela utilizada para o estimador de Horvitz-Thompson. Silva (2004) apresenta maiores detalhes sobre o procedimento em si e sobre a estimação de variância por linearização. Alternativamente, é possível utilizar estimadores de variância baseados em replicação, como bootstrap.

Bootstrap + Raking

A estimação de variância por bootstrap consiste em extrair várias réplicas da amostra obtida utilizando amostragem com reposição. No entanto, a “reamostragem” deve considerar o plano amostral, respeitando a estratificação, conglomeração e probabilidades de inclusão, por exemplo. Além disso, é importante que a fração amostral seja pequena, de modo que se possa aproximar um plano amostral sem reposição por um plano com reposição. Existem variações do bootstrap em amostras complexas, mas os artigo principal é Rao, Wu e Yue (1992)2, que explica o método utilizado pelo IBGE.

Suponha que tenhamos extraído \(R\) réplicas da amostra original. Repetimos o procedimento de estimação em cada uma das \(R\) réplicas. Por exemplo, no caso de uma média, estimamos a estimativa pontual com o peso final da amostra completa. Repetimos este procedimento \(R\) vezes, utilizando as \(R\) réplicas. Utilizando a distribuição das réplicas, estimamos a variância. Note que mesmo no caso de um estimador de médias por razão de totais, não há nenhum procedimento de linearização: o bootstrap contorna a necessidade de linearização, evitando o uso de fórmulas mais complicadas.

Para incorporar o efeito do raking, o procedimento de raking deve ser repetido em cada réplica. Novamente, o uso de fórmulas complicadas é contornado.

Como fazer no R

As divulgações mais recentes do IBGE trazem os pesos finais e os pesos de replicação, todos ajustados por raking. As informações do plano amostral também estão presentes na divulgação, permitindo o uso do Método do Conglomerado Primário. Vamos mostrar os dois casos.

Bootstrap

# carrega libraries
library( survey )

### usando os pesos de replicação

# declarando o plano amostral
pnadc.boot <- svrepdesign( data = pnadc.df , type = "bootstrap" , weights = ~v1028 , repweights = "^v1028[0-9]{3}" , mse = TRUE )

# estima proporções de ocupados e desocupados no Brasil
txunemp.boot.br <- svymean( ~factor( vd4002 ) , pnadc.boot , na.rm = TRUE )

# estima proporções de ocupados e desocupados no Amazonas
txunemp.boot.uf <- svymean( ~factor( vd4002 ) , subset( pnadc.boot , uf == 13 ) , na.rm = TRUE )

A opção mse = TRUE na função svrepdesign é interessante. Temos dois estimadores de variância por Bootstrap:

\[ \begin{align} \widehat{Var}_{Boot1} ( \widehat{\theta} ) &= \frac{1}{R-1} \sum_{ r = 1 }^R ( \widehat{\theta}_{(r)} - \widehat{\overline{\theta}} )^2 , \quad \widehat{\overline{\theta}} = \frac{1}{R} \sum_{ r=1 }^R \widehat{\theta}_{(r)} \\ \widehat{Var}_{Boot2} ( \widehat{\theta} ) &= \frac{1}{R-1} \sum_{ r = 1 }^R ( \widehat{\theta}_{(r)} - \widehat{\theta} )^2 . \end{align} \]

onde \(\widehat{\theta}_{(r)}\) é a estimativa usando o vetor de pesos replicados \(r\) e \(\widehat{\theta}\) é a estimativa usando os pesos finais da amostra completa. Para (uma classe de) estimadores lineares, \(\widehat{\theta} = \widehat{\overline{\theta}}\). Porém, no geral, \(\widehat{\theta} \neq \widehat{\overline{\theta}}\). Ao fazer isso, estamos aproximando a estimação da variância pela estimação do Erro Quadrático Médio (em inglês: Mean Squared Error, MSE). Como os estimadores que usamos são (aproximadamente) não-viesados com amostras grandes, isso significa que, na prática, estamos estimando variâncias. Ou seja: ao estimar variâncias com o estimador \(\widehat{Var}_{Boot2}\), estamos sendo (um pouco) conservadores, já que estamos incluindo um termo quadrático de viés na conta; esse viés, no entanto, desaparece rapidamente com o crescimento do tamanho da amostra. A opção mse = TRUE faz com que o survey use o estimador \(\widehat{Var}_{Boot2}\), mais recomendável3.

Método do Conglomerado Primário

# ajusta formatos
pnadc.df[ , c( "posest" , "posest_sxi" ) ] <- lapply( pnadc.df[ , c( "posest" , "posest_sxi" ) ] , factor )

# declarando o plano amostral
pnadc.ucm <- svydesign( ids = ~upa+v1008 , strata = ~estrato , data = pnadc.df , weights = ~v1028 )

# coleta tabela com os totais das marginais na população
pop.posest <- pnadc.df[ !duplicated( pnadc.df$posest ) , c( "posest" , "v1029" ) ]
pop.posest_sxi <- pnadc.df[ !duplicated( pnadc.df$posest_sxi ) , c( "posest_sxi" , "v1033" ) ]
pop.posest <- pop.posest[ order( pop.posest$posest ) , ]
pop.posest_sxi <- pop.posest_sxi[ order( pop.posest_sxi$posest_sxi ) , ]
colnames( pop.posest )[2] <-"Freq"
colnames( pop.posest_sxi )[2] <-"Freq"

# aplica calibração via raking com limites e pesos fixos no domicílio
pnadc.ucm <- calibrate( pnadc.ucm ,
                        formula  = list( ~posest , ~posest_sxi ), population = list( pop.posest, pop.posest_sxi ) ,
                        bounds = c(.2,5) , bounds.const = FALSE ,
                        calfun = "raking" , aggregate.stage = 2 ,
                        sparse = TRUE )

# estima proporções de ocupados e desocupados no Brasil
txunemp.ucm.br <- svymean( ~factor( vd4002 ) , pnadc.ucm , na.rm = TRUE )

# estima proporções de ocupados e desocupados no Amazonas
txunemp.ucm.uf <- svymean( ~factor( vd4002 ) , subset( pnadc.ucm , uf == 13 ) , na.rm = TRUE )

Comparando resultados

É importante lembrar o que muda em relação aos resultados: o métodos de estimação das variâncias. Ou seja: as estimativas das proporções devem ser idênticas, mas as estimativas de variâncias podem diferir4.

Os resultados para o Brasil são praticamente idênticos:

txunemp.boot.br
##                    mean     SE
## factor(vd4002)1 0.87359 0.0014
## factor(vd4002)2 0.12641 0.0014
txunemp.ucm.br
##                    mean     SE
## factor(vd4002)1 0.87359 0.0014
## factor(vd4002)2 0.12641 0.0014

Já os resultados para o Amazonas diferem um pouco nas estimativas dos erros-padrões:

txunemp.boot.uf
##                    mean     SE
## factor(vd4002)1 0.86553 0.0068
## factor(vd4002)2 0.13447 0.0068
txunemp.ucm.uf
##                    mean     SE
## factor(vd4002)1 0.86553 0.0072
## factor(vd4002)2 0.13447 0.0072

Referências

DEVILLE, J.-C.; SÄRNDAL, C.-E. Calibration Estimators in Survey Sampling. Journal of the American Statistical Association, v. 87, n. 418, p. 376–382, 1992.
DEVILLE, J.-C.; SÄRNDAL, C.-E.; SAUTORY, O. Generalized Raking Procedures in Survey Sampling. Journal of the American Statistical Association, v. 88, n. 423, p. 1013–1020, 1993.
RAO, J. N. K.; WU, C. F. J. Resampling Inference With Complex Survey Data. Journal of the American Statistical Association, v. 83, n. 401, p. 231–241, 1988.
RAO, J. N. K.; WU, C. F. J.; YUE, K. Some Recent work on Resampling Methods for Complex Surveys. Survey Methodology, v. 18, n. 2, p. 209–217, dez. 1992.
SILVA, P. L. N. Calibration Estimation: When and Why, How Much and How: Textos para discussão. Rio de Janeiro: IBGE, Diretoria de Pesquisas, 2004. Disponível em: <https://biblioteca.ibge.gov.br/visualizacao/livros/liv66414.pdf>.
WOLTER, K. M. Introduction to Variance Estimation. Nova York: Springer, 2007.

  1. Simetria e desigualdade triangular — propriedades de distâncias — podem ser desnecessárias neste caso.↩︎

  2. O método anterior de Rao e Wu (1988) tinha um problema: ele não poderia ser aplicado para estimar a variância de estimadores não-regulares, como quantis. O método de Rao, Wu e Yue (1992) funciona para estes casos e é mais fácil de implementar.↩︎

  3. Vide Wolter (2007, p. 214–217).↩︎

  4. Não muito, afinal os métodos tentam estimar a mesma quantidade↩︎