Aula 08

Pós-estratificação

Leituras

Eckman e West (2016)

In the SRS case, we have a list of elements, which allows us to attribute to each element in the list a probability of being selected. In particular, “If the population is size \(N\) and the sample is size \(n\), then the probability of selection of each unit is \(\frac{n}{N}\).” (Eckman e West 2016, 478)

The mean is just:

\[ \bar{y} = \dfrac{1}{n} \sum^n_{i=1} y_i \]

The sample mean has itself a distribution: “We select just one sample, but there are many other samples that we could have selected, and each one would give an estimate of the mean.” The variance of the sample estimate of the mean is:

\[ \text{Var}(\bar{y}) = \dfrac{ \left( 1 - \frac{n}{N} \right) S^2 }{n} \]

where \(\frac{n}{N}\) is the sampling fraction and \(S^2\) is the variability of the variable \((Y)\) in the population. It gives the true variance of the estimate of the mean, but we usually cannot calculate it, because \(S^2\) is not known. We then estimate \(S^2\) by \(s^2\), the in-sample variability:

\[ s^2 = \dfrac{ \sum^n_1 (y_i - \bar{y})^2 }{n-1} \]

\(\left( 1 - \frac{n}{N} \right)\) is the finite population correction, which is often very small. We can estimate the variance as:

\[ \text{var}(\bar{y}) = \dfrac{s^2}{n} \]

Now, in the stratified context, we estimate the variance as:

\[ \text{var}(\bar{y}) = \sum_H \dfrac{W_h^2 s^2_h}{n_h} \]

The first step in deriving weights is to take the inverse of the probability of selection for each unit. In an SRS, the probability of selection of each unit is simply \(n/N\) and thus the weight for each case is \(N/n\). In stratified samples using SRS within strata, but disproportionate allocation to strata, the weight is different for each stratum, but constant within strata. In PPS samples, each unit may have a different probability of selection and thus a different weight. Weights that adjust for the probability of selection are called design weights. (Eckman e West 2016, 482)

Most surveys suffer from nonresponse, and weights are often adjusted to account for it. There are several different methods of adjusting weights for nonresponse using information that is available for respondents and nonrespondents. […]. If we have several variables available for each selected case, we can estimate a response and use the predicted propensities based on the model to adjust the weights. (Eckman e West 2016, 482–83)

Raking, post-stratification, and calibration are in this class: all involve adjusting the sample so that it matches the population based on known population information. For example, in many contries, women are more likely to respond to surveys than men and thus the set of respondents has more women than men. These adjustment techniques increase the weight of the responding men, and decrease the weight of women, so that the weighted sample reflects the population. (Eckman e West 2016, 483)

If these adjustment techniques increase the variance of our estimates, why do we use them? Because we believe that the point estimates (that is, the estimates of means, quantiles, correlation coefficients, regression coefficients, and so on) that we get with the final adjusted weights are more accurate (less biased) than those that we get wehn we use only the design weights. The cost of this reduction in bias is an increase in variance. (Eckman e West 2016, 483)

In complex surveys, weights, stratification, and cluster sampling interact to make the estimation of variances particularly difficult. When a sample caintains these elements, the variance estimation formulas given above do not work, and we need to turn to alternative methods. There are two broad classes of techniques to estimate variances and standard errors from complex samples: Taylor-series linearization and replication, which includes several similar techniques. (Eckman e West 2016, 484)

Replication techniques work on a different principle. Replication methods estimate variance by choosing subsets from the sample itself, and using the variability across all these subsets to estimate the variance of the sampling distribution. How exactly these subsamples of the sample are selected depends on the exact replication method in use. (Eckman e West 2016, 484)

Lumley (2010), cap. 7: Post-stratification, raking, and calibration

This chapter deals with techniques for using known population totals for a set of variables (auxiliary variables) to adjust the sampling weights and improve estimation for another set of variables. All of these techniques have the same idea: adjustments are made to the sampling weights so that estimated population totals for the auxiliary variables match the known population totals, making the sample more representative of the population. A second benefit is that the estimates are forced to be consistent with the population data, improving their credibility with people who may not understand the sampling process. (Lumley 2010, 135–36)

Post-stratification adjusts the sampling weights so that the estimated population group sizes are correct, as they would be in stratified sampling. The sampling weights \(\frac{1}{\pi_i}\) are replaced by weights \(\frac{g_i}{\pi_i}\), where \(g_i = \frac{N_i}{\hat{N}_i}\) for the group containing individual \(i\).

Adjusting the weights to incorporate the auxiliary information is straightforward, but estimating the resulting standard errors is more difficult. For replicate-weight designs it is sufficient to post-stratify each set of replicate weights (Valliant [178]). Since the estimated group size will then match the known population group size for every set of replicates, between-group variation will not produce differences between replicates and so the replicate-weight standard errors will correctly omit the between-group differences. (Lumley 2010, 137)

The function postStratify() creates a post-stratified survey design object. In addition to adjusting the sampling weights, it adds information to allow the standard errors to be adjusted. For a replicate-weight design this involves computing a post-stratification on each set of replicate weights; for linearization estimators it involves storing the group identifier and weight information so that between-group contributions to variance can be removed.

The first step is to set up the information about population group sizes. This information can be in a data frame (as here) or in a table object as is produced by the table() function. In the data frame format, one or more columns give the values of the grouping variables and the final column, which must be named Freq, gives the population counts. The call to postStratify() specifies the grouping variables as a model formula in the strata argument, and gives the table of population counts as the population argument. The function returns a new, post-stratified survey design object.

Raking allows multiple grouping variables to be used without constructing a complete cross-classification. The process involves post-stratifying on each set of variables in turn, and repeating this process until the weights stop changing. The name arises from the image of raking a garden bed alternately in each direction to smooth out the soil. Raking can also be applied when partial cross-classifications of the variables are available.

[…]

Standard error computations after raking are based on the iterative post-stratification algorithm. For standard errors based on replicate weights, raking each set of replicate weights ensures that the auxiliary information is correctly incorporated into between-replicate differences. […].

The rake() function performs the computations for raking by repeatedly calling postStratify(). Each of these calls to postStratify() accumulates the necessary information for standard error computations. (Lumley 2010, 139)

As in section 7.2, post-stratification can be seen as making the smallest possible changes to the weights that result in the estimated population totals matching the known totals. This view of post-stratification leads to calibration estimators. An alternative is to note that the estimated population total after post-stratification is a regression estimator, for a working regression model that uses indicatorvariables for each post-stratum as predictors. This view leads to generalized regression or GREG estimators. (Lumley 2010, 141)

“Regression thinking” is useful in understanding why large increases in precision are possible – it is surprising to many biostatisticians that small changes in sampling weights can increase precision dramatically, but it is easier to understand that a good regression model can reduce the unexplained variation. “Calibration thinking” often leads to simpler formulas and simpler software implementation, since an explicit model is not needed. (Lumley 2010, 141)

\[ \hat{T}_\text{reg} = \sum^N_{i = 1} x_i \hat{\beta}_i = \left( \sum^N_{i = 1} X_i \right) \hat{\beta} = T_X \hat{\beta} \]

Anotações de aula

Em amostragem, cada unidade amostral tem um peso. Lemos isso como “quantas unidades na população cada unidade amostrada representa”. Se o desenho amostral tiver equiprobabilidade, cada unidade amostrada representa a mesma quantidade de unidades da população e pesos não são necessários.

O estimador ponderado é:

\[ \hat{Y} = \sum_{i \in s} w_i y_i \]

onde \(w_i\) é o peso para a unidade \(i\) e \(y_i\) é a medida extraída do survey.

Considere o seguinte exemplo:

Área	Idade	Quantas pessoas ela representa
ZS	25	349
ZS	30	482

Nesse caso, indicamos quantas pessoas cada unidade amostral representa. Esse valor é calculado a partir da probabilidade de inclusão da unidade amostral. O peso é o inverso da probabilidade de inclusão. A ponderação pelo desenho ajusta amostras quando temos probabilidades desiguais de inclusão: \(w_i = \frac{1}{\pi_i}\), onde \(\pi_i\) é a probabilidade de inclusão da unidade amostral \(i\). Supondo que a probabilidade de inclusão é de 0.25, o peso da unidade amostral é 4. Isso significa que ela representa 4 unidades na população.

Esse estimador também é chamado de estimador de expansão, o estimador de Horvitz-Thompson. O peso absoluto me diz quantas pessoas de perfil similar uma única pessoa da amostra representa.

\[ \hat{T} = \sum_{i = 1}^n w_i y_i \]

Note que esse peso é útil para calcular totais populacionais, mas ele é diferente do peso usado para fazer ponderação. Por que ponderar?

• Ponderar pela probabilidade de inclusão quando não temos equiprobabilidade é importante para não incorporar viés na amostra (claro, isso depende de uma amostragem probabilística);
• Dar maior peso a subgrupos que responderam menos ao survey (não-resposta);
• Ajustar estimativas com totais populacionais conhecidos (aka calibração)

Ponderação de ajuste

Problemas com quadro amostral (cobertura, e.g.) e não-resposta de unidade podem introduzir viés se respondentes diferem sistematicamente de não-respondentes:

• Entrevistamos muito menos pessoas mais ricas do que deveríamos;
• Usamos estratificação com alocação não-proporcional (i.e., regiões pequenas e grandes recebem o mesmo número de entrevistas);
• Bolsonaristas ou lulistas se recusam mais a responder.

A ponderação duplica, triplica, às vezes quadruplica a informação que temos. Na prática, isso lembra a amostragem por conglomerados, na perspectiva de que estamos aumentando a margem de erro ao não incorporar informações efetivamente novas.

Ideia básica: ajustar os pesos dos respondentes para compensar as unidades ausentes ou sub-representadas. Pesquisa domiciliar probabilística com \(n=1000\) entrevistou \(380\) homens; na população, homens são \(50\%\):

\[ \frac{0.5}{0.38} \ \text{ou} \ \frac{500}{380} \approx 1.315789474 \]

Cada homem na amostra representa 1.315789474 homens na população. O mesmo vale para mulheres, mas o peso é menor.

Calibração vs. ponderação

Quando já temos pesos (geralmente do desenho), podemos forçá-los a corresponder a totais populacionais conhecidos em variáveis auxiliares.

Ponderar é pegar a amostra que eu tenho e tentar aproximá-la de algum benchmark (por exemplo, a população). Calibrar, por outro lado, é pegar um peso que eu já tenho (no mais das vezes, a probabilidade de inclusão): supondo que a gente calcule o total populacional e, no entanto, a população brasileira deu \(205\) milhões, mas sabemos que a população brasileira é de \(203\) milhões. A calibragem ajusta o esquema de pesos para fazer com que o total bata com algum total conhecido da população.

Nesse caso, precisamos encontrar novos pesos \(w_i^{\text{Cal}}\) próximos aos pesos iniciais de forma que:

\[ \sum_{i \in s} w_i^{\text{Cal}} x_i = \sum_{i \in U} x_i = X, \]

onde \(X\) é o total populacional (\(U\) é a população) conhecido (e.g., total de eleitores) para a variável auxiliar \(x_i\). Esse ajuste pode ser feito de várias maneiras: algoritmos iterativos, regressão, etc.

Pesos e inflação da variância

Pesos de desenho muito grandes ou desiguais inflam estimativas.

Intuição: uma unidade com peso relativo de \(2\) equivale a duplicar a sua entrevista na amostra; para isso, no entanto, precisamos retirar uma unidade de outras unidades da amostra. A solução geralmente é usar tetos ou ajustes para limitar pesos extremos.

Outros tipos de pós-ajustes em R

• Rake e iterative proportional fitting, que ajustam pesos iterativamente para que totais populacionais conhecidos sejam respeitados (pesos relativos);
• Generalized regression estimation (GREG), que usa modelos de regressão para calibrar pesos;
• Multilevel regression and poststratification (MrP), que projeta resposta em subgrupos da população (e.g., regiões, idade, sexo) usando modelos de regressão

Usando o pacote survey no R

library(tidyverse)
library(survey)

set.seed(123)

macae <- readRDS(here::here("trabalhos/tarefa-1/dados/macae.rds"))

verdade <- mean(macae$alfabetizada)
print(verdade)

[1] 0.8720686

De maneira naive, podemos estimar a média apenas tirando a média simples da amostra. Devemos esperar algo razoavelmente destoante do parâmetro populacional.

amostra <- macae %>%
  group_by(subdistrito) %>%
  mutate(pi = 100 / n()) %>%
  slice_sample(n = 100)

mean(amostra$alfabetizada)

[1] 0.8485714

Usando o estimador de Horvitz-Thompson, podemos recuperar os pesos absolutos usando a probabilidade de inclusão (i.e., ponderação pela probabilidade de inclusão):

amostra <- amostra %>%
  mutate(peso = 1 / pi)

Vamos fazer a mesma coisa usando o pacote survey, que nos permite declarar o desenho amostral. Vamos usar as probabilidades de inclusão para “acertar” o cálculo, assim como no exemplo anterior:

# todas as observacoes sao intercambiaveis dentro de um estrato, entao ignoramos id (~1)
metodo1 <- svydesign(~1,
                     strata = ~subdistrito, # estrato
                     probs = ~pi, # probabilidades de inclusao para ponderar
                     data = amostra)

svytotal(~alfabetizada, metodo1)

              total     SE
alfabetizada 176920 3347.8

svymean(~alfabetizada, metodo1)

                mean     SE
alfabetizada 0.85581 0.0162

Agora, vejamos o método 2: não temos a probabilidade de inclusão, mas tenho os totais populacionais extraídos a partir de alguma outra fonte, como, por exemplo, o censo demográfico. Eu sei, por exemplo, quantas pessoas moram em cada distrito.

# MÉTODO 2: pós-estratificação

# vamos ver o caso naive
desenho <- svydesign(~id_pessoa,
                     data = amostra)

# estimamos a média de maneira naive sem ponderar
svymean(~alfabetizada, desenho)

                mean     SE
alfabetizada 0.84857 0.0136

# agora vamos estratificar
estratos <- macae %>%
  group_by(subdistrito) %>%
  summarise(Freq = n())

# criamos um objeto do tipo postStratify
metodo2 <- postStratify(desenho, ~subdistrito, estratos)

# calculamos a media usando esse meotod, que estratificou os resultados
svymean(~alfabetizada, metodo2)

                mean     SE
alfabetizada 0.85581 0.0161

Agora passemos ao método 3, que usa Rake:

desenho <- svydesign(~id_pessoa,
                     data = amostra)

# se fossem varias variaveis desbalanceadas, criaria um banco pra cada e 
# passaria cada um deles na lista
metodo3 <- rake(desenho, 
                list(~subdistrito),
                list(estratos))

svymean(~alfabetizada, metodo3)

                mean     SE
alfabetizada 0.85581 0.0161

Podemos extrair os pesos calculados por cada desenho da seguinte maneira:

amostra <- amostra %>%
  ungroup() %>%
  mutate(
    peso_strat = weights(metodo1),
    peso_pos_strat = weights(metodo2),
    peso_rake = weights(metodo3)
  )

head(amostra)

# A tibble: 6 × 11
  codigo_setor  id_pessoa distrito subdistrito bairro alfabetizada      pi  peso
  <chr>         <glue>    <chr>    <chr>       <chr>         <dbl>   <dbl> <dbl>
1 330240305020… 33024030… Macaé    Aeroporto   Parqu…            1 0.00265  378.
2 330240305020… 33024030… Macaé    Aeroporto   Ajuda             1 0.00265  378.
3 330240305020… 33024030… Macaé    Aeroporto   Ajuda             0 0.00265  378.
4 330240305020… 33024030… Macaé    Aeroporto   Ajuda             1 0.00265  378.
5 330240305020… 33024030… Macaé    Aeroporto   Parqu…            1 0.00265  378.
6 330240305020… 33024030… Macaé    Aeroporto   Parqu…            1 0.00265  378.
# ℹ 3 more variables: peso_strat <dbl>, peso_pos_strat <dbl>, peso_rake <dbl>

Observe que chegamos \(\pm\) aos mesmos pesos com cada approach. De fato, são diferentes maneiras de resolver o mesmo problema.

Simulação

NotaSimulação e texto preparados com auxílio do Claude

Uma amostra de 1.000 pessoas sub-representa homens (38%, como no exemplo da aula) e jovens. Conhecemos as margens da população (sexo e faixa etária) e, para a pós-estratificação, também a tabela cruzada. O raking (iterative proportional fitting) ajusta os pesos alternando entre as margens: primeiro acerta o sexo, depois a idade, depois o sexo de novo, até convergir. Mova o número de iterações para ver as margens ponderadas chegando às da população. O resultado do raking acerta as margens, mas não necessariamente a tabela cruzada; por isso a estimativa da aprovação pode diferir da pós-estratificação quando o desfecho depende da combinação sexo × idade.

Interativo

import {scheme, palette, plotStyle, fmt} from "/_extensions/felipelmc/course-notes/ojs/notes.js"
mode = scheme()
pal = palette(mode)

viewof iter = Inputs.range([0, 8], {value: 0, step: 1, label: "Iterações do raking"})
viewof interacao = Inputs.range([-0.3, 0.3], {value: -0.2, step: 0.01, label: "Efeito extra entre homens jovens"})

celulas = [
  {sexo: "Homens", idade: "18–34", n: 70,  pop: 0.16},
  {sexo: "Homens", idade: "35–59", n: 170, pop: 0.21},
  {sexo: "Homens", idade: "60+",   n: 140, pop: 0.13},
  {sexo: "Mulheres", idade: "18–34", n: 150, pop: 0.14},
  {sexo: "Mulheres", idade: "35–59", n: 280, pop: 0.19},
  {sexo: "Mulheres", idade: "60+",   n: 190, pop: 0.17}
]
aprov = (c) => 0.4 + (c.sexo === "Mulheres" ? 0.15 : 0) + (c.idade === "18–34" ? -0.12 : c.idade === "60+" ? 0.08 : 0)
             + (c.sexo === "Homens" && c.idade === "18–34" ? interacao : 0)
rake = {
  const N = d3.sum(celulas, (c) => c.n);
  const alvoSexo = d3.rollup(celulas, (v) => d3.sum(v, (c) => c.pop), (c) => c.sexo);
  const alvoIdade = d3.rollup(celulas, (v) => d3.sum(v, (c) => c.pop), (c) => c.idade);
  let w = celulas.map(() => 1);
  const passos = [w.slice()];
  for (let k = 0; k < iter; k++) {
    const dim = k % 2 === 0 ? "sexo" : "idade", alvo = dim === "sexo" ? alvoSexo : alvoIdade;
    const atual = d3.rollup(celulas.map((c, i) => ({c, wn: w[i] * c.n})), (v) => d3.sum(v, (d) => d.wn) / N, (d) => d.c[dim]);
    w = w.map((wi, i) => wi * alvo.get(celulas[i][dim]) / atual.get(celulas[i][dim]));
  }
  const margens = [];
  for (const dim of ["sexo", "idade"]) {
    const alvo = dim === "sexo" ? alvoSexo : alvoIdade;
    const pond = d3.rollup(celulas.map((c, i) => ({c, wn: w[i] * c.n})), (v) => d3.sum(v, (d) => d.wn) / N, (d) => d.c[dim]);
    for (const [cat, p] of alvo) margens.push({categoria: cat, fonte: "População", p}, {categoria: cat, fonte: "Amostra ponderada", p: pond.get(cat)});
  }
  const est = (pesos) => d3.sum(celulas, (c, i) => pesos[i] * c.n * aprov(c)) / d3.sum(celulas, (c, i) => pesos[i] * c.n);
  const ps = celulas.map((c) => c.pop * N / c.n);
  return {margens, bruta: est(celulas.map(() => 1)), raking: est(w), pos: est(ps), verdade: d3.sum(celulas, (c) => c.pop * aprov(c))};
}

Plot.plot({
  width, height: 220, style: plotStyle, marginLeft: 70,
  x: {label: "Proporção →", domain: [0, 0.65], tickFormat: ".0%"}, y: {label: null, domain: ["Homens", "Mulheres", "18–34", "35–59", "60+"]},
  color: {domain: ["População", "Amostra ponderada"], range: [pal.award, pal.accent], legend: true},
  marks: [
    Plot.link(d3.groups(rake.margens, (d) => d.categoria).map(([c, v]) => ({c, x1: v[0].p, x2: v[1].p})), {y1: "c", y2: "c", x1: "x1", x2: "x2", stroke: pal.lineStrong}),
    Plot.dot(rake.margens, {y: "categoria", x: "p", fill: "fonte", r: 5})
  ]
})

html`<p class="cn-readout">Aprovação estimada: sem pesos <b>${fmt(100 * rake.bruta, 1)}%</b> · raking (${iter} it.) <b>${fmt(100 * rake.raking, 1)}%</b> · pós-estratificação <b>${fmt(100 * rake.pos, 1)}%</b> · população <b>${fmt(100 * rake.verdade, 1)}%</b></p>`

Em R

library(survey)
amostra <- data.frame(
  sexo  = rep(c("H", "H", "H", "M", "M", "M"), c(70, 170, 140, 150, 280, 190)),
  idade = rep(c("18-34", "35-59", "60+"), 2)[rep(1:6, c(70, 170, 140, 150, 280, 190))]
)
d <- svydesign(ids = ~1, data = amostra)
pop_sexo  <- data.frame(sexo = c("H", "M"), Freq = c(500, 500))
pop_idade <- data.frame(idade = c("18-34", "35-59", "60+"), Freq = c(300, 400, 300))
d_rake <- rake(d, list(~sexo, ~idade), list(pop_sexo, pop_idade))
svymean(~sexo + idade, d_rake)

Eckman, Stephanie, e Brady T. West. 2016. “Analysis of Data from Stratified and Clustered Surveys”. Cap. 31 em The SAGE Handbook of Survey Methodology, organizado por Christof Wolf, Dominique Joye, Tom W. Smith, e Yang-chih Fu. SAGE. https://doi.org/10.4135/9781473957893.n31.

Lumley, Thomas. 2010. Complex Surveys: A Guide to Analysis Using R. John Wiley & Sons. https://doi.org/10.1002/9780470580066.