Estimação em Pequenas Áreas

Nessa parte do trabalho, utilizaremos o pacote sae, o qual pode ter sua documentação consultada em https://cran.r-project.org/web/packages/sae/sae.pdf.

# install.packages("sae")
library("sae")

Incomedata

Dados sintéticos sobre renda e outras variáveis relacionadas para províncias espanholas.

Data frame com 17.199 observações nas seguintes 21 variáveis.

---

• provlab: nome da província.

• prov: código da província.

• ac: região da província.

• gen: gender: 1: masculino, 2: feminino.

• age: faixat-etária: \((0: \leq 13), (1: 14 - 15), (2: 16 - 24), (3: 25 - 49), (4: 50 - 64), (5: \geq 65).\)

• nat: nacionalidade: 1: espanhol, 2: outros.

• educ: nível de escolaridade: (0:idade \(\geq16\)), 1:ensino primário (ensino obrigatório), 2:ensino secundária, 3:ensino pós-secundário..

• labor: situação no mercado de trabalho: 0:idade \(\geq16\), 1:empregado, 2:desempregado, 3:inativo.

• age2: indicador de grupo de idade 16-24.

• age3: indicador de grupo de idade 25-49.

• age4: indicador de grupo de idade 50-64.

• age5: indicador de grupo de idade \(\geq 65\).

• educ1: indicador de nível de educação 1 (educação primária).

• educ2: indicador de nível de educação 2 (educação secundária).

• educ3: indicador de nível de educação 3 (educação pós-secundária).

• nat1: indicador de nacionalidade espanhola..

• labor1: indicador de estar empregado.

• labor2: indicador de estar desempregado.

• labor3: indicador de estar inativo.

• income: renda normalizada.

• weight: peso amostral.

data("incomedata")
attach(incomedata)
incomedata |> rmarkdown::paged_table()

sizeprov

Identificadores e tamanhos populacionais para os domínios no conjunto de dados incomedata.

Data frame com 52 observações nas seguintes 3 variáveis:

• provlab: nome da província.

• prov: código da província.

• Nd: contagem da população da província

data(sizeprov)
data("sizeprov")
attach(sizeprov)
sizeprov |> rmarkdown::paged_table()

sizeprovage

Nomes, códigos e tamanhos populacionais por faixa etária para domínios no conjunto de dados incomedata.

Data frame com 52 observações nas seguintes 7 variáveis:

provlab: nome da província.

prov: código da província.

age1: contagem populacional da província para a faixa etária \(<16\).

age2: contagem populacional da província para a faixa etária \(16-24\).

age3: contagem populacional da província para a faixa etária \(25-49\).

age4: contagem populacional da província para a faixa etária \(50-64\).

age5: contagem populacional da província para a faixa etária \(\geq 65\).

data(sizeprovage)
data("sizeprovage")
attach(sizeprovage)
sizeprovage |> rmarkdown::paged_table()

sizeprovedu

Identificadores e tamanhos populacionais por nível de educação para domínios no conjunto de dados incomedata.

Data frame com 52 observações nas seguintes 6 variáveis:

provlab: nome da província.

prov: código da província.

educ0: contagem populacional da província para o nível de educação 0 (idade <16).

educ1: contagem populacional da província para o nível de educação 1 (educação primária).

educ2: contagem populacional da província para o nível de educação 2 (educação secundária).

educ3: contagem populacional da província para o nível de educação 3 (educação pós-secundária).

data(sizeprovedu)
data("sizeprovedu")
attach(sizeprovedu)
sizeprovedu |> rmarkdown::paged_table()

sizeprovlab

Nomes, códigos e tamanhos populacionais por status da força de trabalho para domínios no conjunto de dados incomedata.

Data frame com 52 observações nas seguintes 6 variáveis:

provlab: nome da província.

prov: código da província.

labor0: contagem populacional da província para o status da força de trabalho 0 (idade <16).

labor1: contagem populacional da província para o status da força de trabalho 1 (empregado).

labor2: contagem populacional da província para o status da força de trabalho (desempregado).

labor3: contagem populacional da província para o status da força de trabalho 3 (inativo).

data(sizeprovlab)
data("sizeprovlab")
attach(sizeprovlab)
sizeprovlab |> rmarkdown::paged_table()

n = dim(incomedata)[1];n # tamanho total da amostra

[1] 17199

D = length(unique(prov));D # numero de provincias (areas ou dominios) 

[1] 52

incomedata$prov = as.factor(incomedata$prov)
nd = as.vector(table(incomedata$prov));nd # tamanho da amostra de cada provincia

 [1]   96  173  539  198   58  494  634 1420  168  282  398  118  250  224  495
[16]   92  142  208   89  285  122  115  232  218  130  510  173  944  379  885
[31]  564  129  803   72  472  448  164  381  434   58  482   20  134   72  275
[46]  714  299  524  104  564  235  180

Nd = sizeprov$Nd; Nd # tamanho da populacao de cada provincia

 [1]  296558  381413 1717556  618150  163142  660203  984224 5149877  353099
[10]  403966 1165911  542448  495305  771665 1106660  206479  657764  866032
[19]  205620  677697  479951  213217  647120  479017  398604  298760  346779
[28] 5921832 1443395 1334159  581532  327805 1051532  168113 1008195  922947
[37]  342223  942587  553771  153559 1784244   90064  707683  138975  594895
[46] 2376984  506269 1122978  193664  895909   70981   65336

z = 6557.143 # linha de pobreza
pobreza = numeric(n)
pobreza[income < z] = 1
sum(pobreza == 1)

[1] 3841

Estimador Direto Horvitz Thompson

direct(y, dom, sweight, domsize, data, replace = FALSE)

HT = direct(y = pobreza, dom = incomedata$provlab, 
            sweight = incomedata$weight, domsize = sizeprov[, c("provlab", "Nd")])

HT |> rmarkdown::paged_table()

cve = HT$CV
sum(cve > 20)

[1] 15

sum(cve > 15)

[1] 29

sum(cve > 10)

[1] 44

sum(cve > 5)

[1] 52

Estimador Sintético Pós Estratificado

A seguir, iremos calcular estimativas sintéticas pós-estratificadas com os níveis de educação como pós-estratos. Para a função pssynt(), construímos o data frame domsizebyps, contendo os códigos de domínio provlab na primeira coluna e, nas colunas restantes, os tamanhos das províncias por nível de educação.

pssynt(y, sweight, ps, domsizebyps, data)

Popn.educ = sizeprovedu[,-2]
colnames(Popn.educ) = c("provincia","0", "1", "2", "3")
PSYN.educ = pssynt(y = pobreza, sweight = incomedata$weight, 
                   ps = incomedata$educ,domsizebyps = Popn.educ)

PSYN.educ |> rmarkdown::paged_table()

Estimador Composto

Calculamos as estimativas SSD pela composição das estimativas diretas e pós-estratificadas anteriores, utilizando o valor padrão delta=1 na função ssd().

ssd(dom, sweight, domsize, direct, synthetic, delta = 1, data)

SSD = ssd(dom = provlab, sweight = weight, domsize = sizeprov[, c("provlab", "Nd")],
          direct = HT[, c("Domain", "Direct")], 
          synthetic = PSYN.educ, delta = 1, data = incomedata)

SSD |> rmarkdown::paged_table()

resultados = data.frame(Province = HT$Domain, 
                        SampleSize = HT$SampSize, 
                        HT = HT$Direct * 100, 
                        PSYN.educ = PSYN.educ$PsSynthetic * 100, 
                        SSD = SSD$ssd * 100)

resultados |> rmarkdown::paged_table()

resultados = resultados[order(resultados$SampleSize,decreasing = TRUE), ]

# Ajuste os parâmetros de margem para aumentar o espaço ao redor do gráfico
par(mar = c(6, 5, 4, 2) + 0.1)

plot(resultados$HT, type = "n",
     xlab = "Província", 
     ylab = "Estimativa", cex.axis = 1.5, cex.lab = 1.5,
     xlim = c(0, 52),
     ylim = c(0, 40))

# Adicione pontos ao gráfico com cores e símbolos específicos
points(resultados$HT, type = "b", col = 'blue', lwd = 2, pch = 1)
points(resultados$PSYN.educ, type = "b", col = 'green', lwd = 2, pch = 2)
points(resultados$SSD, type = "b", col = 'red', lwd = 2, pch = 5)

# Ajuste a legenda para refletir as cores e símbolos corretos
legend("bottom", legend = c("HT", "Post-strat", "Composto"),
       col = c('blue', 'green', 'red'), lwd = rep(2, 3),
       pch = c(1, 2, 5), cex = 0.8)

---

Modelo Fay Herriot

# número de áreas
D = 50 
set.seed(1612)

x = sapply(1:2, function (i) {runif (n = D, min = 10, max = 20)});x

          [,1]     [,2]
 [1,] 18.85623 13.88190
 [2,] 18.10490 16.21155
 [3,] 11.63356 15.53969
 [4,] 13.72649 17.06431
 [5,] 17.53557 16.78931
 [6,] 10.40188 16.05491
 [7,] 10.53416 14.14101
 [8,] 19.89670 15.04002
 [9,] 16.06765 15.91225
[10,] 10.37912 11.22895
[11,] 12.16259 14.74180
[12,] 14.28087 12.42848
[13,] 18.26524 13.31538
[14,] 12.65394 18.99317
[15,] 10.70181 11.24130
[16,] 10.64762 17.03197
[17,] 11.01198 14.97359
[18,] 16.57665 19.48965
[19,] 14.26830 19.83772
[20,] 12.46373 16.25510
[21,] 12.08773 14.76526
[22,] 15.57703 12.19107
[23,] 16.48102 16.50580
[24,] 15.71326 18.46680
[25,] 16.07614 10.73179
[26,] 10.80556 15.40777
[27,] 17.09642 10.21305
[28,] 13.79047 10.71636
[29,] 12.23637 10.96988
[30,] 15.58827 16.31309
[31,] 16.46522 17.58170
[32,] 15.08033 12.67266
[33,] 16.74110 13.74074
[34,] 11.68258 11.64713
[35,] 10.47391 10.52326
[36,] 10.23308 15.86040
[37,] 10.61581 18.24448
[38,] 12.40970 14.28638
[39,] 12.42063 13.25670
[40,] 10.10880 19.78122
[41,] 15.60544 16.38225
[42,] 13.26826 17.95389
[43,] 17.71834 17.43805
[44,] 15.23757 13.82779
[45,] 14.45048 14.36708
[46,] 18.34556 19.87367
[47,] 11.77179 19.43011
[48,] 12.59975 11.70115
[49,] 19.73666 14.95643
[50,] 19.67324 10.59646

x.com.intercepto = cbind(1, x);x.com.intercepto

      [,1]     [,2]     [,3]
 [1,]    1 18.85623 13.88190
 [2,]    1 18.10490 16.21155
 [3,]    1 11.63356 15.53969
 [4,]    1 13.72649 17.06431
 [5,]    1 17.53557 16.78931
 [6,]    1 10.40188 16.05491
 [7,]    1 10.53416 14.14101
 [8,]    1 19.89670 15.04002
 [9,]    1 16.06765 15.91225
[10,]    1 10.37912 11.22895
[11,]    1 12.16259 14.74180
[12,]    1 14.28087 12.42848
[13,]    1 18.26524 13.31538
[14,]    1 12.65394 18.99317
[15,]    1 10.70181 11.24130
[16,]    1 10.64762 17.03197
[17,]    1 11.01198 14.97359
[18,]    1 16.57665 19.48965
[19,]    1 14.26830 19.83772
[20,]    1 12.46373 16.25510
[21,]    1 12.08773 14.76526
[22,]    1 15.57703 12.19107
[23,]    1 16.48102 16.50580
[24,]    1 15.71326 18.46680
[25,]    1 16.07614 10.73179
[26,]    1 10.80556 15.40777
[27,]    1 17.09642 10.21305
[28,]    1 13.79047 10.71636
[29,]    1 12.23637 10.96988
[30,]    1 15.58827 16.31309
[31,]    1 16.46522 17.58170
[32,]    1 15.08033 12.67266
[33,]    1 16.74110 13.74074
[34,]    1 11.68258 11.64713
[35,]    1 10.47391 10.52326
[36,]    1 10.23308 15.86040
[37,]    1 10.61581 18.24448
[38,]    1 12.40970 14.28638
[39,]    1 12.42063 13.25670
[40,]    1 10.10880 19.78122
[41,]    1 15.60544 16.38225
[42,]    1 13.26826 17.95389
[43,]    1 17.71834 17.43805
[44,]    1 15.23757 13.82779
[45,]    1 14.45048 14.36708
[46,]    1 18.34556 19.87367
[47,]    1 11.77179 19.43011
[48,]    1 12.59975 11.70115
[49,]    1 19.73666 14.95643
[50,]    1 19.67324 10.59646

beta = c(3, 2, 4)
sigma2.ed = rep(2,D)
sigma2.u  = 2

e.d = sapply(1:D, function (d) { 
  rnorm(1, mean = 0, sd = sqrt(sigma2.ed[d]))   # erros de amostragem
}); e.d

 [1] -0.855265059  1.347414906 -0.803040712 -0.314362798  1.014075226
 [6] -0.778488942  2.489122735 -0.822802867  0.903084810 -0.097091544
[11] -0.292957343 -0.263253761 -4.390509015 -0.550808746  0.557036548
[16] -1.507622510 -3.163748598  0.549885112  0.426986351 -3.087092531
[21]  1.266963100  0.607454972 -0.204064561  2.302759640  2.572582828
[26]  0.248219479 -1.105059217 -0.560733057 -0.134716131 -1.632545314
[31]  1.564220778  3.392639534  0.737392606  1.675110486 -1.931030234
[36]  1.831067412 -0.357295519  3.015400012 -0.167365175  0.485608392
[41]  1.079764018 -0.091780332 -0.451288102 -2.878136161 -0.007060956
[46]  0.599907666  1.732854958 -0.484226942 -0.118353345  0.026801836

u.d = rnorm(D, mean = 0, sd = sqrt(sigma2.u)); u.d  # efeitos aleatórios

 [1]  1.16047148  0.87480755  1.29945118  1.10521252  0.50646047  0.47972070
 [7]  0.23453647 -0.19730232 -0.81232036  2.32259411  0.78594184 -0.26826850
[13]  1.17762900  1.04841874 -2.16878988 -0.89278836  0.76609900  1.40694357
[19]  0.71141104 -0.81703855  0.26937099 -0.94023380  0.97905188  0.40547655
[25]  0.89767967 -1.34365594  0.56016096  3.22892088  1.33758118 -2.66401636
[31]  0.74309796 -0.09075346  0.38279443 -0.97480403  1.69748446  1.95036806
[37]  2.25973636 -1.56030113 -1.45387241 -0.04081634  0.81514663  0.76633411
[43] -1.10081415 -0.33229212  0.70201217 -1.68816325  0.85934256 -0.19598090
[49] -0.89576418  1.26289568

mu = x.com.intercepto %*% beta + u.d # parametros de interesse verdadeiros
mu.dir = mu + e.d # estimativas diretas

mu

           [,1]
 [1,]  97.40054
 [2,] 104.93081
 [3,]  89.72531
 [4,]  99.81544
 [5,] 105.73486
 [6,]  88.50313
 [7,]  80.86690
 [8,] 102.75617
 [9,]  97.97200
[10,]  70.99665
[11,]  87.07835
[12,]  81.00741
[13,]  93.96965
[14,] 105.32898
[15,]  67.20002
[16,]  91.53032
[17,]  85.68441
[18,] 115.51885
[19,] 111.59888
[20,]  92.13083
[21,]  86.50587
[22,]  81.97809
[23,] 102.96428
[24,] 108.69920
[25,]  78.97712
[26,]  84.89856
[27,]  78.60520
[28,]  76.67529
[29,]  72.68986
[30,]  96.76488
[31,] 107.00032
[32,]  83.76054
[33,]  91.82795
[34,]  71.97887
[35,]  67.73834
[36,]  88.85814
[37,]  99.46928
[38,]  83.40462
[39,]  79.41419
[40,] 102.30169
[41,] 100.55503
[42,] 102.11843
[43,] 107.08805
[44,]  88.45402
[45,]  90.07130
[46,] 117.49763
[47,] 105.12336
[48,]  74.80813
[49,] 101.40328
[50,]  85.99520

mu.dir

           [,1]
 [1,]  96.54528
 [2,] 106.27823
 [3,]  88.92227
 [4,]  99.50108
 [5,] 106.74894
 [6,]  87.72464
 [7,]  83.35602
 [8,] 101.93336
 [9,]  98.87508
[10,]  70.89956
[11,]  86.78539
[12,]  80.74416
[13,]  89.57914
[14,] 104.77817
[15,]  67.75706
[16,]  90.02270
[17,]  82.52067
[18,] 116.06873
[19,] 112.02587
[20,]  89.04374
[21,]  87.77283
[22,]  82.58555
[23,] 102.76021
[24,] 111.00196
[25,]  81.54970
[26,]  85.14678
[27,]  77.50014
[28,]  76.11456
[29,]  72.55514
[30,]  95.13233
[31,] 108.56454
[32,]  87.15318
[33,]  92.56534
[34,]  73.65398
[35,]  65.80731
[36,]  90.68920
[37,]  99.11199
[38,]  86.42002
[39,]  79.24683
[40,] 102.78729
[41,] 101.63480
[42,] 102.02665
[43,] 106.63676
[44,]  85.57588
[45,]  90.06424
[46,] 118.09754
[47,] 106.85622
[48,]  74.32390
[49,] 101.28493
[50,]  86.02200

library(sae)

FH = mseFH(mu.dir ~ x, vardir = sigma2.ed, method = "REML")
FH

$est
$est$eblup
        [,1]
1   96.46004
2  105.23301
3   88.91198
4   99.30715
5  106.08826
6   88.16149
7   82.19305
8  102.52894
9   98.98352
10  69.96994
11  86.75712
12  81.18045
13  91.33684
14 104.78055
15  68.82295
16  91.56097
17  84.02834
18 115.22913
19 111.66732
20  91.31906
21  87.20198
22  82.88739
23 102.50869
24 109.77855
25  79.84105
26  85.98389
27  77.85915
28  74.87419
29  72.11597
30  97.54484
31 107.52482
32  85.57412
33  92.09953
34  73.49612
35  66.15158
36  88.99853
37  98.41805
38  85.87960
39  80.30091
40 102.88090
41 100.81670
42 101.91258
43 107.58220
44  87.40367
45  89.87678
46 118.81232
47 105.78931
48  74.86699
49 101.88002
50  85.37135

$est$fit
$est$fit$method
[1] "REML"

$est$fit$convergence
[1] TRUE

$est$fit$iterations
[1] 2

$est$fit$estcoef
                 beta  std.error    tvalue       pvalue
X(Intercept) 3.656570 2.00555690  1.823219 6.827021e-02
Xx1          1.958174 0.09440393 20.742502 1.432766e-95
Xx2          4.019737 0.09905509 40.580820 0.000000e+00

$est$fit$refvar
[1] 1.832603

$est$fit$goodness
  loglike       AIC       BIC       KIC 
-103.0355  214.0711  221.7192  218.0711 



$mse
 [1] 1.117027 1.101065 1.077597 1.072587 1.095571 1.099029 1.095161 1.139295
 [9] 1.071140 1.133267 1.070966 1.078547 1.108943 1.109394 1.127381 1.102900
[17] 1.085431 1.127057 1.122911 1.072359 1.071691 1.086828 1.078950 1.097686
[25] 1.118777 1.089269 1.143629 1.110255 1.113094 1.069713 1.090264 1.077332
[33] 1.080875 1.105686 1.146939 1.101182 1.120831 1.069813 1.076392 1.163795
[41] 1.070313 1.086179 1.105682 1.067121 1.061959 1.163820 1.128105 1.095793
[49] 1.134980 1.187304

cve.dir = sqrt(sigma2.ed)/mu.dir*100
cve.fh  = sqrt(FH$mse)/FH$est$eblup*100

cve.dir

          [,1]
 [1,] 1.464819
 [2,] 1.330671
 [3,] 1.590393
 [4,] 1.421305
 [5,] 1.324803
 [6,] 1.612105
 [7,] 1.696594
 [8,] 1.387390
 [9,] 1.430303
[10,] 1.994672
[11,] 1.629553
[12,] 1.751475
[13,] 1.578731
[14,] 1.349722
[15,] 2.087183
[16,] 1.570952
[17,] 1.713769
[18,] 1.218428
[19,] 1.262399
[20,] 1.588224
[21,] 1.611220
[22,] 1.712423
[23,] 1.376227
[24,] 1.274044
[25,] 1.734174
[26,] 1.660913
[27,] 1.824788
[28,] 1.858007
[29,] 1.949157
[30,] 1.486575
[31,] 1.302648
[32,] 1.622676
[33,] 1.527800
[34,] 1.920078
[35,] 2.149022
[36,] 1.559407
[37,] 1.426884
[38,] 1.636442
[39,] 1.784568
[40,] 1.375864
[41,] 1.391466
[42,] 1.386122
[43,] 1.326197
[44,] 1.652584
[45,] 1.570228
[46,] 1.197496
[47,] 1.323473
[48,] 1.902771
[49,] 1.396272
[50,] 1.644014

cve.fh

        [,1]
1  1.0956818
2  0.9971362
3  1.1675294
4  1.0428831
5  0.9866269
6  1.1891197
7  1.2732212
8  1.0410500
9  1.0455872
10 1.5214395
11 1.1928411
12 1.2792875
13 1.1529454
14 1.0052227
15 1.5427732
16 1.1469847
17 1.2398675
18 0.9213202
19 0.9489573
20 1.1339887
21 1.1871581
22 1.2577430
23 1.0133043
24 0.9543806
25 1.3247854
26 1.2138093
27 1.3735138
28 1.4072758
29 1.4629671
30 1.0603014
31 0.9710847
32 1.2129205
33 1.1288346
34 1.4307097
35 1.6189367
36 1.1790893
37 1.0757104
38 1.2043812
39 1.2920066
40 1.0485844
41 1.0261786
42 1.0226401
43 0.9774054
44 1.1818901
45 1.1465853
46 0.9079909
47 1.0039982
48 1.3982146
49 1.0456951
50 1.2763471

library(ggplot2) 

df = data.frame("CVE" = c(cve.dir, cve.fh), "Estimador" = c(rep("Direto", D),
                rep("FH", D)), "Área" = c(1:D, 1:D))

# gráfico
ggplot(df, aes(x = Área, y = CVE, lty = Estimador, color = Estimador)) +
  geom_point(pch = 19) + 
  geom_line() +
  theme(legend.position = "bottom") +
  scale_y_continuous("Coeficiente de Variação Estimada (CVE)")