1 Introducción

Este documento replica los ejemplos del artículo de Molina y Marhuenda (2015), que presenta el paquete sae de R para la estimación en áreas pequeñas (Small Area Estimation, SAE). Se cubren cinco enfoques:

Sección	Modelo	Nivel	Función principal	Datos
2	Fay-Herriot (FH)	Área	`mseFH()`	`milk`
3	FH espacial (SFH)	Área	`mseSFH()`	`grapes`, `grapesprox`
4	FH espacio-temporal (STFH)	Área	`pbmseSTFH()`	`spacetime`, `spacetimeprox`
5	Battese-Harter-Fuller (BHF)	Unidad	`pbmseBHF()`	`cornsoybean`, `cornsoybeanmeans`
6	EB para parámetros no lineales	Unidad	`ebBHF()`, `pbmseebBHF()`	`incomedata`, `Xoutsamp`

En todos los casos se comparan los estimadores directos con los estimadores basados en modelos (EBLUP o EB) mediante el coeficiente de variación (CV):

\[ \text{CV}(\hat{\delta}_d) = 100 \cdot \frac{\sqrt{\widehat{\text{MSE}}(\hat{\delta}_d)}}{\hat{\delta}_d}. \]

# Paquete con los modelos de estimación en áreas pequeñas y sus datos de ejemplo
library(sae)

# Función para graficar el estimador directo y el estimador del modelo
# sobre las mismas áreas (las áreas deben venir ya ordenadas).
#   directo         : valores del estimador directo (o su CV)
#   modelo          : valores del estimador del modelo (o su CV)
#   etiqueta_modelo : nombre del modelo que aparece en la leyenda
#   ylab, xlab      : títulos de los ejes
#   ylim            : rango del eje y (si no se da, se calcula solo)
#   tipo            : "b" = puntos unidos con líneas, "p" = solo puntos
#   etiquetas_x     : nombres para el eje x (opcional)
graficar_comparacion <- function(directo, modelo, etiqueta_modelo,
                                 ylab, xlab, ylim = NULL, tipo = "b",
                                 etiquetas_x = NULL) {
  # Si no se da el rango del eje y, se usa el mínimo y máximo de los datos
  if (is.null(ylim)) ylim <- range(c(directo, modelo), na.rm = TRUE)
  # Gráfico vacío (type = "n") para luego agregar las dos series
  plot(directo, type = "n", ylab = ylab, xlab = xlab, ylim = ylim,
       xaxt = if (is.null(etiquetas_x)) "s" else "n")
  # Si hay etiquetas, se dibuja el eje x con esos nombres en vertical (las = 2)
  if (!is.null(etiquetas_x)) axis(1, seq_along(etiquetas_x), etiquetas_x, las = 2, cex.axis = 0.8)
  # Estimador directo: negro, círculos, línea continua
  points(directo, type = tipo, col = 1, lwd = 2, pch = 1, lty = 1)
  # Estimador del modelo: azul, equis, línea punteada
  points(modelo,  type = tipo, col = 4, lwd = 2, pch = 4, lty = 2)
  # Leyenda en la esquina superior derecha, con fondo blanco
  legend("topright", legend = c("Directo", etiqueta_modelo), inset = 0.02,
         bg = "white", col = c(1, 4), lwd = 2, pch = c(1, 4),
         lty = if (tipo == "p") 0 else c(1, 2), bty = "n")
}

2 EBLUP basado en el modelo de Fay-Herriot

2.1 Ejemplo: gasto promedio en leche

El conjunto de datos milk contiene estimadores directos del gasto promedio en leche fresca en 1989 para 43 áreas pequeñas. Sus variables son:

SmallArea: áreas de interés inferencial.
ni: tamaño de muestra de cada área.
yi: estimador directo del gasto promedio en leche en 1989.
SD: desviación estándar estimada de los estimadores directos.
CV: coeficiente de variación de los estimadores directos.
MajorArea: áreas mayores, que se usan como variable auxiliar (efectos fijos).

# Cargar los datos de ejemplo del paquete
data("milk")
milk

2.2 Ajuste del modelo

Se ajusta el modelo FH con MajorArea como factor y las varianzas muestrales \(\psi_d = SD_d^2\) como varianzas conocidas de los errores de muestreo.

# La varianza muestral del estimador directo es la desviación estándar al cuadrado.
# Se guarda como columna porque, al usar `data`, la función la busca por nombre.
milk$var <- milk$SD^2

# mseFH() ajusta el modelo FH y calcula el EBLUP y su error cuadrático medio (MSE).
#   yi ~ as.factor(MajorArea) : fórmula; yi es el estimador directo y MajorArea
#                               la variable auxiliar (como categorías)
#   var                       : varianza muestral del estimador directo de cada área
#   data                      : base de datos donde están las variables
FH <- mseFH(yi ~ as.factor(MajorArea), var, data = milk)

# Coeficientes del modelo (beta), con su error estándar y valor p
FH$est$fit$estcoef

# Varianza del efecto aleatorio del área (sigma^2_u)
FH$est$fit$refvar

## [1] 0.01855022

# Medidas de bondad de ajuste (log-verosimilitud, AIC, BIC)
FH$est$fit$goodness

##    loglike        AIC        BIC        KIC 
##  12.677478 -15.354956  -6.548956 -10.354956

2.3 Resultados

# CV del EBLUP en porcentaje: raíz del MSE dividida entre la estimación
cv.FH <- 100 * sqrt(FH$mse) / FH$est$eblup

# Tabla con el estimador directo y el EBLUP, cada uno con su CV
results.FH <- data.frame(
  Area       = milk$SmallArea,  # área
  SampleSize = milk$ni,         # tamaño de muestra
  DIR        = milk$yi,         # estimador directo
  cv.DIR     = 100 * milk$CV,   # CV del directo (viene en proporción, se pasa a %)
  eblup.FH   = FH$est$eblup,    # EBLUP del modelo FH
  cv.FH      = cv.FH            # CV del EBLUP
)
results.FH

# Ordenar las áreas de mayor a menor tamaño de muestra
results.FH <- results.FH[order(results.FH$SampleSize, decreasing = TRUE), ]

# Dos gráficos lado a lado (1 fila, 2 columnas)
par(mfrow = c(1, 2))
# Izquierda: estimaciones
graficar_comparacion(results.FH$DIR, results.FH$eblup.FH, "EBLUP FH",
                     ylab = "Estimación", ylim = c(0.4, 1.6),
                     xlab = "Área (tamaño de muestra decreciente)")
# Derecha: coeficientes de variación
graficar_comparacion(results.FH$cv.DIR, results.FH$cv.FH, "EBLUP FH",
                     ylab = "CV (%)", ylim = c(5, 40),
                     xlab = "Área (tamaño de muestra decreciente)")

Figura 1. Estimaciones (izquierda) y CV (derecha) del estimador directo y del EBLUP FH, con las áreas ordenadas por tamaño de muestra decreciente.

# Volver a un gráfico por figura
par(mfrow = c(1, 1))

3 EBLUP basado en el modelo FH espacial

3.1 Ejemplo: superficie promedio dedicada al cultivo de uva

El conjunto de datos grapes contiene información de 274 municipios de la Toscana (Italia):

grapehect: estimador directo de la superficie promedio (en hectáreas) dedicada a la producción de uva en cada municipio.
area: superficie agraria utilizada para la producción.
workdays: número promedio de días de trabajo en el año de referencia.
var: varianza muestral de los estimadores directos de cada municipio.

El conjunto de datos grapesprox contiene la matriz de proximidad (estandarizada por filas) que representa la vecindad entre los municipios de la Toscana. Esta matriz permite modelar la correlación espacial de los efectos de área mediante un proceso SAR (Simultaneous Autoregressive).

# Datos de los municipios y matriz de vecindad entre ellos
data("grapes")
data("grapesprox")
grapes

# Primeras filas de la matriz de proximidad (cada fila suma 1)
head(grapesprox)

3.2 Ajuste del modelo

Se calculan los EBLUP espaciales de la superficie promedio dedicada a la producción de uva.

# mseSFH() ajusta el modelo FH espacial y calcula el EBLUP y su MSE.
#   grapehect ~ area + workdays - 1 : fórmula; grapehect es el estimador directo,
#                                     area y workdays son variables auxiliares
#                                     y "- 1" quita el intercepto
#   var        : varianza muestral del estimador directo de cada municipio
#   grapesprox : matriz de proximidad (quién es vecino de quién)
#   data       : base de datos donde están las variables
SFH <- mseSFH(grapehect ~ area + workdays - 1, var, grapesprox, data = grapes)

# Coeficientes del modelo (beta)
SFH$est$fit$estcoef

# Varianza del efecto aleatorio del área
SFH$est$fit$refvar

## [1] 69.74899

# Autocorrelación espacial (rho): qué tanto se parecen los municipios vecinos
SFH$est$fit$spatialcorr

## [1] 0.6142697

3.3 Resultados

# CV del EBLUP espacial en porcentaje
cv.SFH <- 100 * sqrt(SFH$mse) / SFH$est$eblup

# Tabla con el estimador directo y el EBLUP espacial, cada uno con su CV
results.SFH <- data.frame(
  DIR       = grapes$grapehect,                                   # estimador directo
  cv.DIR    = 100 * abs(sqrt(grapes$var) / grapes$grapehect),     # CV del directo
  eblup.SFH = SFH$est$eblup,                                      # EBLUP espacial
  cv.SFH    = cv.SFH                                              # CV del EBLUP
)
results.SFH

# Ordenar los municipios de menor a mayor CV del estimador directo
results.SFH <- results.SFH[order(results.SFH$cv.DIR), ]

par(mfrow = c(1, 2))
# Izquierda: estimaciones (solo puntos, porque son muchos municipios)
graficar_comparacion(results.SFH$DIR, results.SFH$eblup.SFH, "EBLUP SFH",
                     ylab = "Estimación", ylim = c(0, 400), tipo = "p",
                     xlab = "Área (CV creciente del estimador directo)")
# Derecha: coeficientes de variación
graficar_comparacion(results.SFH$cv.DIR, results.SFH$cv.SFH, "EBLUP SFH",
                     ylab = "CV (%)", ylim = c(0, 400), tipo = "p",
                     xlab = "Área (CV creciente del estimador directo)")

Figura 2. Estimaciones (izquierda) y CV (derecha) del estimador directo y del EBLUP SFH, con los municipios ordenados por CV creciente del estimador directo.

par(mfrow = c(1, 1))

4 EBLUP basado en el modelo FH espacio-temporal

4.1 Ejemplo: datos sintéticos

El conjunto de datos spacetime contiene datos sintéticos para \(D = 11\) áreas en \(T = 3\) periodos de tiempo. Sus variables son:

Area: código del área.
Time: periodo de tiempo.
X1 y X2: variables auxiliares para cada área y periodo.
Y: estimador directo de cada área y periodo.
Var: varianza muestral de los estimadores directos.

El conjunto spacetimeprox contiene la matriz de proximidad entre las 11 áreas.

# Datos por área y periodo, y matriz de vecindad entre las áreas
data("spacetime")
data("spacetimeprox")
spacetime

4.2 Ajuste del modelo

Se calculan los EBLUP de las medias de cada área y periodo. El MSE se estima mediante bootstrap paramétrico, por lo que se fija una semilla para garantizar la reproducibilidad.

D  <- nrow(spacetimeprox)            # Número de áreas
TT <- length(unique(spacetime$Time)) # Número de periodos (se evita `T`, que en R es TRUE)

# Semilla: el MSE se calcula por simulación (bootstrap); así siempre da lo mismo
set.seed(123)

# pbmseSTFH() ajusta el modelo FH espacio-temporal y calcula el EBLUP
# y su MSE por bootstrap paramétrico.
#   Y ~ X1 + X2   : fórmula; Y es el estimador directo, X1 y X2 variables auxiliares
#   D             : número de áreas
#   TT            : número de periodos de tiempo
#   vardir = Var  : varianza muestral del estimador directo de cada área y periodo
#   spacetimeprox : matriz de proximidad entre las áreas
#   data          : base de datos donde están las variables
STFH <- pbmseSTFH(Y ~ X1 + X2, D, TT, vardir = Var, spacetimeprox, data = spacetime)

## 
## Bootstrap procedure with B = 100 iterations starts.
## b = 1 
## b = 2 
## b = 3 
## b = 4 
## b = 5 
## b = 6 
## b = 7 
## b = 8 
## b = 9 
## b = 10 
## b = 11 
## b = 12 
## b = 13 
## b = 14 
## b = 15 
## b = 16 
## b = 17 
## b = 18 
## b = 19 
## b = 20 
## b = 21 
## b = 22 
## b = 23 
## b = 24 
## b = 25 
## b = 26 
## b = 27 
## b = 28 
## b = 29 
## b = 30 
## b = 31 
## b = 32 
## b = 33 
## b = 34 
## b = 35 
## b = 36 
## b = 37 
## b = 38 
## b = 39 
## b = 40 
## b = 41 
## b = 42 
## b = 43 
## b = 44 
## b = 45 
## b = 46 
## b = 47 
## b = 48 
## b = 49 
## b = 50 
## b = 51 
## b = 52 
## b = 53 
## b = 54 
## b = 55 
## b = 56 
## b = 57 
## b = 58 
## b = 59 
## b = 60 
## b = 61 
## b = 62 
## b = 63 
## b = 64 
## b = 65 
## b = 66 
## b = 67 
## b = 68 
## b = 69 
## b = 70 
## b = 71 
## b = 72 
## b = 73 
## b = 74 
## b = 75 
## b = 76 
## b = 77 
## b = 78 
## b = 79 
## b = 80 
## b = 81 
## b = 82 
## b = 83 
## b = 84 
## b = 85 
## b = 86 
## b = 87 
## b = 88 
## b = 89 
## b = 90 
## b = 91 
## b = 92 
## b = 93 
## b = 94 
## b = 95 
## b = 96 
## b = 97 
## b = 98 
## b = 99 
## b = 100

# Coeficientes del modelo (beta)
STFH$est$fit$estcoef

# Componentes de varianza: varianza del efecto de área, autocorrelación espacial,
# varianza del efecto área-tiempo y autocorrelación en el tiempo
STFH$est$fit$estvarcomp

4.3 Resultados

# CV del EBLUP espacio-temporal y CV del estimador directo, en porcentaje
cv.STFH <- 100 * sqrt(STFH$mse) / STFH$est$eblup
cv.DIR  <- 100 * sqrt(spacetime$Var) / spacetime$Y

# Tabla con ambos estimadores y su CV para cada área y periodo
results.STFH <- data.frame(
  Area       = spacetime$Area,   # área
  Time       = spacetime$Time,   # periodo
  DIR        = spacetime$Y,      # estimador directo
  eblup.STFH = STFH$est$eblup,   # EBLUP espacio-temporal
  cv.DIR     = cv.DIR,           # CV del directo
  cv.STFH    = cv.STFH           # CV del EBLUP
)

# Quedarse solo con el último periodo (t = 3)
results.lasttime <- results.STFH[results.STFH$Time == TT, ]
results.lasttime

# Ordenar las áreas de menor a mayor CV del estimador directo
results.lasttime <- results.lasttime[order(results.lasttime$cv.DIR), ]

par(mfrow = c(1, 2))
# Izquierda: estimaciones, con el código de cada área en el eje x
graficar_comparacion(results.lasttime$DIR, results.lasttime$eblup.STFH, "EBLUP STFH",
                     ylab = "Estimación", ylim = c(0.05, 0.45),
                     xlab = "Área (t = 3)", etiquetas_x = results.lasttime$Area)
# Derecha: coeficientes de variación
graficar_comparacion(results.lasttime$cv.DIR, results.lasttime$cv.STFH, "EBLUP STFH",
                     ylab = "CV (%)",
                     xlab = "Área (t = 3)", etiquetas_x = results.lasttime$Area)

Figura 3. Estimaciones (izquierda) y CV (derecha) del estimador directo y del EBLUP STFH en el último periodo, con las áreas ordenadas por CV creciente del estimador directo.

par(mfrow = c(1, 1))

5 EBLUP basado en el modelo BHF

5.1 Ejemplo: hectáreas de maíz en condados de Iowa

Se usan los datos de producción de maíz y soya en 12 condados de Iowa, distribuidos en dos conjuntos de datos. Los condados son los dominios y los segmentos muestrales son las unidades.

El conjunto cornsoybean contiene, para cada segmento muestral:

County: código del condado.
CornHec: hectáreas de maíz reportadas en la encuesta.
SoyBeansHec: hectáreas de soya reportadas en la encuesta.
CornPix: número de píxeles clasificados como maíz según datos satelitales.
SoyBeansPix: número de píxeles clasificados como soya según datos satelitales.

Para obtener los EBLUP basados en el modelo BHF se necesitan, además, las medias poblacionales de las variables auxiliares y el tamaño poblacional \(N_d\) de cada condado. Estos valores están en cornsoybeanmeans:

SampSegments: número de segmentos muestrales del condado (tamaño de muestra).
PopnSegments: número de segmentos poblacionales del condado (tamaño de la población).
MeanCornPixPerSeg: número promedio de píxeles de maíz por segmento.
MeanSoyBeansPixPerSeg: número promedio de píxeles de soya por segmento.

# Datos de la muestra (un registro por segmento) y datos poblacionales por condado
data("cornsoybean")
data("cornsoybeanmeans")
cornsoybeanmeans

5.2 Ajuste del modelo

Se calculan los EBLUP de las medias de hectáreas de maíz por condado, usando como variables auxiliares el número de píxeles de maíz y de soya.

# Medias poblacionales de las variables auxiliares en cada condado
# (la primera columna debe ser el código del condado)
Xmean <- data.frame(cornsoybeanmeans[, c("CountyIndex", "MeanCornPixPerSeg", "MeanSoyBeansPixPerSeg")])
# Tamaño de la población (número total de segmentos) de cada condado
Popn  <- data.frame(cornsoybeanmeans[, c("CountyIndex", "PopnSegments")])

Siguiendo a Battese, Harter y Fuller (1988), se elimina la observación 33, que corresponde a un segmento atípico.

# Quitar el segmento atípico (fila 33)
cornsoybean <- cornsoybean[-33, ]

# Semilla: el MSE se calcula por bootstrap
set.seed(123)

# pbmseBHF() ajusta el modelo BHF (a nivel de unidad) y calcula el EBLUP
# de la media de cada condado y su MSE por bootstrap paramétrico.
#   CornHec ~ CornPix + SoyBeansPix : fórmula; hectáreas de maíz explicadas por
#                                     los píxeles de maíz y de soya
#   dom      : variable que indica el dominio (condado) de cada segmento
#   meanxpop : medias poblacionales de las variables auxiliares por condado
#   popnsize : tamaño de la población de cada condado
#   B        : número de réplicas bootstrap para estimar el MSE
#   data     : base de datos de la muestra
BHF <- pbmseBHF(CornHec ~ CornPix + SoyBeansPix, dom = County,
                meanxpop = Xmean, popnsize = Popn, B = 200, data = cornsoybean)

## 
## Bootstrap procedure with B = 200 iterations starts.
## b = 1 
## b = 2 
## b = 3 
## b = 4 
## b = 5 
## b = 6 
## b = 7 
## b = 8 
## b = 9 
## b = 10 
## b = 11 
## b = 12

## b = 13 
## b = 14 
## b = 15 
## b = 16 
## b = 17 
## b = 18 
## b = 19 
## b = 20 
## b = 21 
## b = 22 
## b = 23 
## b = 24 
## b = 25 
## b = 26 
## b = 27 
## b = 28 
## b = 29 
## b = 30 
## b = 31 
## b = 32 
## b = 33 
## b = 34 
## b = 35 
## b = 36 
## b = 37 
## b = 38 
## b = 39 
## b = 40 
## b = 41 
## b = 42 
## b = 43 
## b = 44 
## b = 45 
## b = 46 
## b = 47 
## b = 48 
## b = 49 
## b = 50 
## b = 51 
## b = 52 
## b = 53 
## b = 54 
## b = 55 
## b = 56 
## b = 57 
## b = 58 
## b = 59 
## b = 60 
## b = 61 
## b = 62 
## b = 63 
## b = 64 
## b = 65 
## b = 66 
## b = 67 
## b = 68 
## b = 69 
## b = 70 
## b = 71 
## b = 72 
## b = 73 
## b = 74 
## b = 75 
## b = 76 
## b = 77 
## b = 78 
## b = 79 
## b = 80 
## b = 81 
## b = 82 
## b = 83 
## b = 84 
## b = 85 
## b = 86 
## b = 87 
## b = 88 
## b = 89 
## b = 90

## b = 91 
## b = 92 
## b = 93 
## b = 94 
## b = 95 
## b = 96 
## b = 97 
## b = 98 
## b = 99 
## b = 100 
## b = 101

## b = 102 
## b = 103 
## b = 104 
## b = 105 
## b = 106 
## b = 107 
## b = 108 
## b = 109 
## b = 110 
## b = 111 
## b = 112 
## b = 113 
## b = 114 
## b = 115 
## b = 116 
## b = 117 
## b = 118 
## b = 119 
## b = 120 
## b = 121 
## b = 122 
## b = 123 
## b = 124 
## b = 125

## b = 126 
## b = 127 
## b = 128 
## b = 129 
## b = 130 
## b = 131 
## b = 132 
## b = 133 
## b = 134 
## b = 135 
## b = 136 
## b = 137 
## b = 138 
## b = 139 
## b = 140 
## b = 141 
## b = 142 
## b = 143 
## b = 144 
## b = 145 
## b = 146 
## b = 147

## b = 148 
## b = 149 
## b = 150 
## b = 151 
## b = 152 
## b = 153 
## b = 154 
## b = 155 
## b = 156 
## b = 157 
## b = 158 
## b = 159 
## b = 160 
## b = 161 
## b = 162 
## b = 163 
## b = 164 
## b = 165 
## b = 166 
## b = 167 
## b = 168 
## b = 169 
## b = 170 
## b = 171 
## b = 172 
## b = 173 
## b = 174 
## b = 175 
## b = 176 
## b = 177 
## b = 178

## b = 179 
## b = 180 
## b = 181 
## b = 182 
## b = 183 
## b = 184 
## b = 185 
## b = 186 
## b = 187 
## b = 188 
## b = 189 
## b = 190 
## b = 191 
## b = 192 
## b = 193 
## b = 194 
## b = 195 
## b = 196 
## b = 197 
## b = 198 
## b = 199 
## b = 200

# Coeficientes del modelo (beta)
BHF$est$fit$fixed

## XsXs(Intercept)     XsXsCornPix XsXsSoyBeansPix 
##      51.0703981       0.3287217      -0.1345684

# Varianza del efecto aleatorio del condado
BHF$est$fit$refvar

## [1] 140.0239

# Varianza del error individual (de cada segmento)
BHF$est$fit$errorvar

## [1] 147.2686

5.3 Resultados

# CV del EBLUP en porcentaje
cv.BHF <- 100 * sqrt(BHF$mse$mse) / BHF$est$eblup$eblup

# Tabla con el EBLUP de cada condado y su CV
results.BHF <- data.frame(
  CountyIndex = BHF$est$eblup$domain,     # código del condado
  CountyName  = cornsoybeanmeans$CountyName, # nombre del condado
  SampleSize  = BHF$est$eblup$sampsize,   # tamaño de muestra
  eblup.BHF   = BHF$est$eblup$eblup,      # EBLUP del modelo BHF
  cv.BHF      = cv.BHF                    # CV del EBLUP
)
results.BHF

5.4 Comparación con el estimador directo

Como referencia se calcula el estimador directo de la media de hectáreas de maíz en cada condado con direct(). Al no especificar pesos, la función supone muestreo aleatorio simple; con replace = FALSE se aplica la corrección por población finita. En los condados con un solo segmento muestral no es posible estimar la varianza, por lo que su CV directo queda como NA.

# direct() calcula el estimador directo (media muestral) de cada dominio.
#   y        : variable de interés (hectáreas de maíz de cada segmento)
#   dom      : condado al que pertenece cada segmento
#   domsize  : tamaño de la población de cada condado
#   replace  : FALSE = muestreo sin reemplazo (corrige por población finita)
# Como no se dan pesos (sweight), se supone muestreo aleatorio simple.
DIR.corn <- direct(y = cornsoybean$CornHec, dom = cornsoybean$County,
                   domsize = Popn, replace = FALSE)

# Pegar el estimador directo y su CV a la tabla de resultados,
# emparejando por código de condado (direct() ya da el CV en porcentaje)
results.BHF$DIR    <- DIR.corn$Direct[match(results.BHF$CountyIndex, DIR.corn$Domain)]
results.BHF$cv.DIR <- DIR.corn$CV[match(results.BHF$CountyIndex, DIR.corn$Domain)]
results.BHF[, c("CountyName", "SampleSize", "DIR", "cv.DIR", "eblup.BHF", "cv.BHF")]

# Ordenar los condados de menor a mayor tamaño de muestra
results.BHF <- results.BHF[order(results.BHF$SampleSize), ]

# mar: márgenes del gráfico; se agranda el de abajo para que quepan los nombres
par(mfrow = c(1, 2), mar = c(7, 4, 2, 1))
# Izquierda: estimaciones, con el nombre de cada condado en el eje x
graficar_comparacion(results.BHF$DIR, results.BHF$eblup.BHF, "EBLUP BHF",
                     ylab = "Hectáreas de maíz", xlab = "",
                     ylim = c(60, 190), etiquetas_x = results.BHF$CountyName)
# Derecha: coeficientes de variación (los condados con n = 1 no tienen CV directo)
graficar_comparacion(results.BHF$cv.DIR, results.BHF$cv.BHF, "EBLUP BHF",
                     ylab = "CV (%)", xlab = "", ylim = c(0, 30),
                     etiquetas_x = results.BHF$CountyName)

Figura 4. Hectáreas de maíz estimadas (izquierda) y CV (derecha) del estimador directo y del EBLUP BHF, con los condados ordenados por tamaño de muestra creciente.

# Volver a un gráfico por figura y a los márgenes por defecto
par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

6 Estimadores EB de parámetros no lineales basados en el modelo BHF

6.1 Ejemplo: mapeo de pobreza

En este ejemplo se ilustra cómo estimar la incidencia de pobreza en las provincias (áreas) de España. Se usa el conjunto de datos incomedata, que contiene datos sintéticos de ingreso y otras variables sociológicas de las provincias españolas. Estos datos se generaron por simulación con el único propósito de ilustrar el uso del paquete sae, de modo que los resultados no son reales. Las variables que se usan son:

provlab: nombre de la provincia.
prov: código de la provincia.
income: ingreso.
weight: peso muestral.
educ: nivel educativo.
labor: situación laboral.
educ1, educ3, labor1, labor2: variables indicadoras de las categorías de nivel educativo y situación laboral.

# Datos de la muestra (un registro por persona) y datos de las personas
# fuera de la muestra en las provincias seleccionadas
data("incomedata")
data("Xoutsamp")
head(incomedata)

6.2 Transformación de la variable respuesta

El método EB supone que la variable respuesta del modelo BHF tiene una distribución aproximadamente normal. Sin embargo, el histograma del ingreso muestra una fuerte asimetría a la derecha, por lo que se necesita una transformación para lograr una normalidad aproximada.

# Histograma del ingreso para ver si es simétrico (prob = TRUE: escala de densidad)
hist(incomedata$income, breaks = 50, prob = TRUE, main = "",
     xlab = "Ingreso", ylab = "Densidad", col = "grey85", border = "white")

Distribución del ingreso (escala original).

Se usa la transformación logarítmica \(\log(\text{income} + c)\), que pertenece tanto a la familia Box-Cox como a la familia de potencias (con \(\lambda = 0\)). Para elegir la constante \(c\) se probó una grilla de valores en el rango del ingreso: para cada valor de \(c\) se ajustó un modelo BHF a \(\log(\text{income} + c)\) y se seleccionó el valor con el que la distribución de los residuos es aproximadamente simétrica. El valor resultante fue \(c = 3500\).

6.3 Definición del parámetro objetivo

Las funciones del paquete que implementan el método EB pueden estimar cualquier parámetro de dominio, siempre que sea una función \(\delta_d = h(\mathbf{O}_d)\) de los valores \(O_{dj}\) cuya transformación \(Y_{dj} = T(O_{dj})\) es la variable respuesta del modelo BHF. Basta con definir ese parámetro como una función de R.

En este ejemplo, el parámetro objetivo es la incidencia de pobreza (proporción de personas con ingreso por debajo de la línea de pobreza), que es una función del ingreso. Se define la función povertyincidence() con la línea de pobreza \(z = 6557.143\):

\[ F_{0d} = \frac{1}{N_d} \sum_{j=1}^{N_d} I(E_{dj} < z). \]

# Incidencia de pobreza: proporción de personas con ingreso menor que la
# línea de pobreza (6557.143). `y < 6557.143` da TRUE/FALSE y mean() cuenta
# la proporción de TRUE.
povertyincidence <- function(y) {
  mean(y < 6557.143)
}

6.4 Datos fuera de muestra

Al estimar parámetros no lineales se necesitan los valores de las variables auxiliares para cada unidad fuera de la muestra. Aunque el modelo se ajusta con la muestra de todas las provincias, para ahorrar tiempo de cómputo aquí se calculan los estimadores EB y su MSE solo para las 5 provincias con los tamaños de muestra más pequeños. Para estas provincias, el conjunto Xoutsamp contiene los valores de las variables auxiliares (categorías de nivel educativo y situación laboral) de cada individuo fuera de la muestra.

# Códigos de las 5 provincias que se van a estimar
provincecodes   <- unique(Xoutsamp$domain)
# Nombres de esas provincias
provincelabels  <- unique(incomedata$provlab)[provincecodes]
# Variables auxiliares de las personas fuera de la muestra
# (la primera columna debe ser el código de la provincia)
Xoutsamp_AuxVar <- Xoutsamp[, c("domain", "educ1", "educ3", "labor1", "labor2")]

data.frame(Codigo = provincecodes, Provincia = provincelabels)

6.5 Ajuste del modelo y estimadores EB

# Semilla: los estimadores EB se calculan por simulación (Monte Carlo)
set.seed(123)

# ebBHF() ajusta el modelo BHF a log(ingreso + constante) y calcula el
# estimador EB del parámetro definido en `indicator`.
#   income ~ educ1 + educ3 + labor1 + labor2 : fórmula; ingreso explicado por
#                                              educación y situación laboral
#   dom        : variable que indica la provincia de cada persona
#   selectdom  : provincias para las que se quiere la estimación
#   Xnonsample : variables auxiliares de las personas fuera de la muestra
#   MC         : número de simulaciones Monte Carlo para el estimador EB
#   constant   : constante c de la transformación log(income + c)
#   indicator  : función con el parámetro a estimar (incidencia de pobreza)
#   data       : base de datos de la muestra
EB <- ebBHF(income ~ educ1 + educ3 + labor1 + labor2, dom = prov,
            selectdom = provincecodes, Xnonsample = Xoutsamp_AuxVar, MC = 50,
            constant = 3500, indicator = povertyincidence, data = incomedata)

# Resumen del modelo ajustado (coeficientes y varianzas)
EB$fit$summary

## Linear mixed model fit by REML ['lmerMod']
## Formula: ys ~ -1 + Xs + (1 | dom)
## 
## REML criterion at convergence: 18966.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.2201 -0.6617  0.0204  0.6882  3.5797 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  dom      (Intercept) 0.008904 0.09436 
##  Residual             0.174676 0.41794 
## Number of obs: 17199, groups:  dom, 52
## 
## Fixed effects:
##                Estimate Std. Error t value
## Xs(Intercept)  9.505176   0.014385 660.781
## Xseduc1       -0.124043   0.007281 -17.036
## Xseduc3        0.291927   0.010366  28.161
## Xslabor1       0.145985   0.006916  21.108
## Xslabor2      -0.081624   0.017083  -4.778
## 
## Correlation of Fixed Effects:
##          Xs(In) Xsedc1 Xsedc3 Xslbr1
## Xseduc1  -0.212                     
## Xseduc3  -0.070  0.206              
## Xslabor1 -0.199  0.128 -0.228       
## Xslabor2 -0.079  0.039 -0.039  0.168

# Estimadores EB de la incidencia de pobreza en cada provincia
EB$eb

6.6 Diagnóstico de residuos

par(mfrow = c(1, 2))
# Izquierda: residuos en orden; deben estar alrededor de 0 sin patrones
plot(EB$fit$residuals, xlab = "Índice", ylab = "Residuos",
     ylim = c(-2, 2), col = 4)
abline(h = 0)  # línea de referencia en 0
# Derecha: histograma de los residuos; debe verse simétrico
hist(EB$fit$residuals, prob = TRUE, xlab = "Residuos", ylab = "Densidad",
     main = "", xlim = c(-2, 2), ylim = c(0, 1), col = "grey85", border = "white")

Figura 5. Índice de los residuos (izquierda) e histograma de los residuos (derecha) del modelo BHF ajustado a log(income + 3500).

par(mfrow = c(1, 1))

6.7 Estimación del MSE por bootstrap paramétrico

# Semilla: el MSE se calcula por bootstrap
set.seed(123)

# pbmseebBHF() calcula los estimadores EB y además su MSE por bootstrap.
# Usa los mismos argumentos que ebBHF() y agrega:
#   B : número de réplicas bootstrap para estimar el MSE
pbmse.EB <- pbmseebBHF(income ~ educ1 + educ3 + labor1 + labor2, dom = prov,
                       selectdom = provincecodes, Xnonsample = Xoutsamp_AuxVar,
                       B = 200, MC = 50, constant = 3500,
                       indicator = povertyincidence, data = incomedata)

## 
## Bootstrap procedure with B = 200 iterations starts.
## b = 1 
## b = 2 
## b = 3 
## b = 4 
## b = 5 
## b = 6 
## b = 7 
## b = 8 
## b = 9 
## b = 10 
## b = 11 
## b = 12 
## b = 13 
## b = 14 
## b = 15 
## b = 16 
## b = 17 
## b = 18 
## b = 19 
## b = 20 
## b = 21 
## b = 22 
## b = 23 
## b = 24 
## b = 25 
## b = 26 
## b = 27 
## b = 28 
## b = 29 
## b = 30 
## b = 31 
## b = 32 
## b = 33 
## b = 34 
## b = 35 
## b = 36 
## b = 37 
## b = 38 
## b = 39 
## b = 40 
## b = 41 
## b = 42 
## b = 43 
## b = 44 
## b = 45 
## b = 46 
## b = 47 
## b = 48 
## b = 49 
## b = 50 
## b = 51 
## b = 52 
## b = 53 
## b = 54 
## b = 55 
## b = 56 
## b = 57 
## b = 58 
## b = 59 
## b = 60 
## b = 61 
## b = 62 
## b = 63 
## b = 64 
## b = 65 
## b = 66 
## b = 67 
## b = 68 
## b = 69 
## b = 70 
## b = 71 
## b = 72 
## b = 73 
## b = 74 
## b = 75 
## b = 76 
## b = 77 
## b = 78 
## b = 79 
## b = 80 
## b = 81 
## b = 82 
## b = 83 
## b = 84 
## b = 85 
## b = 86 
## b = 87 
## b = 88 
## b = 89 
## b = 90 
## b = 91 
## b = 92 
## b = 93 
## b = 94 
## b = 95 
## b = 96 
## b = 97 
## b = 98 
## b = 99 
## b = 100 
## b = 101 
## b = 102 
## b = 103 
## b = 104 
## b = 105 
## b = 106 
## b = 107 
## b = 108 
## b = 109 
## b = 110 
## b = 111 
## b = 112 
## b = 113 
## b = 114 
## b = 115 
## b = 116 
## b = 117 
## b = 118 
## b = 119 
## b = 120 
## b = 121 
## b = 122 
## b = 123 
## b = 124 
## b = 125 
## b = 126 
## b = 127 
## b = 128 
## b = 129 
## b = 130 
## b = 131 
## b = 132 
## b = 133 
## b = 134 
## b = 135 
## b = 136 
## b = 137 
## b = 138 
## b = 139 
## b = 140 
## b = 141 
## b = 142 
## b = 143 
## b = 144 
## b = 145 
## b = 146 
## b = 147 
## b = 148 
## b = 149 
## b = 150 
## b = 151 
## b = 152 
## b = 153 
## b = 154 
## b = 155 
## b = 156 
## b = 157 
## b = 158 
## b = 159 
## b = 160 
## b = 161 
## b = 162 
## b = 163 
## b = 164 
## b = 165 
## b = 166 
## b = 167 
## b = 168 
## b = 169 
## b = 170 
## b = 171 
## b = 172 
## b = 173 
## b = 174 
## b = 175 
## b = 176 
## b = 177 
## b = 178 
## b = 179 
## b = 180 
## b = 181 
## b = 182 
## b = 183 
## b = 184 
## b = 185 
## b = 186 
## b = 187 
## b = 188 
## b = 189 
## b = 190 
## b = 191 
## b = 192 
## b = 193 
## b = 194 
## b = 195 
## b = 196 
## b = 197 
## b = 198 
## b = 199 
## b = 200

6.8 Resultados

# CV del estimador EB en porcentaje
pbcv.EB <- 100 * sqrt(pbmse.EB$mse$mse) / pbmse.EB$est$eb$eb

# Tabla con el estimador EB de cada provincia y su CV
results.EB <- data.frame(
  ProvinceIndex = pbmse.EB$est$eb$domain,    # código de la provincia
  ProvinceName  = provincelabels,            # nombre de la provincia
  SampleSize    = pbmse.EB$est$eb$sampsize,  # tamaño de muestra
  EB            = pbmse.EB$est$eb$eb,        # estimador EB de la incidencia de pobreza
  cv.EB         = pbcv.EB                    # CV del estimador EB
)
results.EB

6.9 Comparación con el estimador directo

Para medir la ganancia del método EB se calcula el estimador directo de la incidencia de pobreza con la función direct(). Este estimador es la media ponderada (con los pesos muestrales) de la variable indicadora \(I(E_{dj} < z)\). Los tamaños poblacionales \(N_d\) de las provincias están en el conjunto de datos sizeprov.

# Tamaño de la población de cada provincia
data("sizeprov")

# Variable indicadora de pobreza: 1 si el ingreso está bajo la línea, 0 si no
pobre <- as.integer(incomedata$income < 6557.143)

# direct() calcula el estimador directo de cada provincia.
#   y       : variable de interés (indicadora de pobreza); su media es la incidencia
#   dom     : provincia de cada persona
#   sweight : peso muestral de cada persona
#   domsize : tamaño de la población de cada provincia (código y N_d)
DIR.pov <- direct(y = pobre, dom = incomedata$prov, sweight = incomedata$weight,
                  domsize = sizeprov[, c("prov", "Nd")])

# Quedarse solo con las 5 provincias seleccionadas, en el mismo orden de results.EB
DIR.pov <- DIR.pov[match(results.EB$ProvinceIndex, DIR.pov$Domain), ]

# Tabla con ambos estimadores y su CV
comparacion.EB <- data.frame(
  ProvinceName = results.EB$ProvinceName,  # nombre de la provincia
  SampleSize   = results.EB$SampleSize,    # tamaño de muestra
  DIR          = DIR.pov$Direct,           # estimador directo
  cv.DIR       = DIR.pov$CV,               # CV del directo (direct() ya lo da en %)
  EB           = results.EB$EB,            # estimador EB
  cv.EB        = results.EB$cv.EB          # CV del EB
)
comparacion.EB

# Ordenar las provincias de menor a mayor tamaño de muestra
comparacion.EB <- comparacion.EB[order(comparacion.EB$SampleSize), ]

# Margen inferior más grande para que quepan los nombres de las provincias
par(mfrow = c(1, 2), mar = c(7, 4, 2, 1))
# Izquierda: incidencia de pobreza estimada
graficar_comparacion(comparacion.EB$DIR, comparacion.EB$EB, "EB",
                     ylab = "Incidencia de pobreza", xlab = "",
                     ylim = c(0, 0.45), etiquetas_x = comparacion.EB$ProvinceName)
# Derecha: coeficientes de variación
graficar_comparacion(comparacion.EB$cv.DIR, comparacion.EB$cv.EB, "EB",
                     ylab = "CV (%)", xlab = "", ylim = c(0, 110),
                     etiquetas_x = comparacion.EB$ProvinceName)

Figura 6. Incidencia de pobreza estimada (izquierda) y CV (derecha) del estimador directo y del estimador EB, con las provincias ordenadas por tamaño de muestra creciente.

par(mfrow = c(1, 1), mar = c(5, 4, 4, 2) + 0.1)

7 Referencias

Battese, G. E., Harter, R. M. y Fuller, W. A. (1988). An error-components model for prediction of county crop areas using survey and satellite data. Journal of the American Statistical Association, 83(401), 28–36.
Fay, R. E. y Herriot, R. A. (1979). Estimates of income for small places: an application of James-Stein procedures to census data. Journal of the American Statistical Association, 74(366), 269–277.
Molina, I. y Rao, J. N. K. (2010). Small area estimation of poverty indicators. Canadian Journal of Statistics, 38(3), 369–385.
Molina, I. y Marhuenda, Y. (2015). sae: An R package for small area estimation. The R Journal, 7(1), 81–98.
Rao, J. N. K. y Molina, I. (2015). Small Area Estimation (2.ª ed.). Wiley.

8 Información de la sesión

# Versión de R y de los paquetes usados (para reproducibilidad)
sessionInfo()

## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.6.2
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] sae_1.3      lme4_2.1-0   Matrix_1.7-5 MASS_7.3-65 
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.7.3      nlme_3.1-169     cli_3.6.6        knitr_1.52      
##  [5] rlang_1.3.0      xfun_0.61        reformulas_0.4.4 jsonlite_2.0.0  
##  [9] minqa_1.2.8      htmltools_0.5.9  sass_0.4.10      rmarkdown_2.32  
## [13] grid_4.6.1       evaluate_1.0.5   jquerylib_0.1.4  fastmap_1.2.0   
## [17] yaml_2.3.12      lifecycle_1.0.5  compiler_4.6.1   Rcpp_1.1.2      
## [21] lattice_0.22-9   digest_0.6.39    nloptr_2.2.1     R6_2.6.1        
## [25] Rdpack_2.6.6     splines_4.6.1    rbibutils_2.4.1  bslib_0.12.0    
## [29] tools_4.6.1      boot_1.3-32      cachem_1.1.0

Estimación en áreas pequeñas con el paquete `sae`

Réplica de los ejemplos de Molina y Marhuenda (2015)

Yeison Buitrago

2026-10-27

1 Introducción

2 EBLUP basado en el modelo de Fay-Herriot

2.1 Ejemplo: gasto promedio en leche

2.2 Ajuste del modelo

2.3 Resultados

3 EBLUP basado en el modelo FH espacial

3.1 Ejemplo: superficie promedio dedicada al cultivo de uva

3.2 Ajuste del modelo

3.3 Resultados

4 EBLUP basado en el modelo FH espacio-temporal

4.1 Ejemplo: datos sintéticos

4.2 Ajuste del modelo

4.3 Resultados

5 EBLUP basado en el modelo BHF

5.1 Ejemplo: hectáreas de maíz en condados de Iowa

5.2 Ajuste del modelo

5.3 Resultados

5.4 Comparación con el estimador directo

6 Estimadores EB de parámetros no lineales basados en el modelo BHF

6.1 Ejemplo: mapeo de pobreza

6.2 Transformación de la variable respuesta

6.3 Definición del parámetro objetivo

6.4 Datos fuera de muestra

6.5 Ajuste del modelo y estimadores EB

6.6 Diagnóstico de residuos

6.7 Estimación del MSE por bootstrap paramétrico

6.8 Resultados

6.9 Comparación con el estimador directo

7 Referencias

8 Información de la sesión

Estimación en áreas pequeñas con el paquete sae

Réplica de los ejemplos de Molina y Marhuenda (2015)

Yeison Buitrago

2026-10-27

1 Introducción

2 EBLUP basado en el modelo de Fay-Herriot

2.1 Ejemplo: gasto promedio en leche

2.2 Ajuste del modelo

2.3 Resultados

3 EBLUP basado en el modelo FH espacial

3.1 Ejemplo: superficie promedio dedicada al cultivo de uva

3.2 Ajuste del modelo

3.3 Resultados

4 EBLUP basado en el modelo FH espacio-temporal

4.1 Ejemplo: datos sintéticos

4.2 Ajuste del modelo

4.3 Resultados

5 EBLUP basado en el modelo BHF

5.1 Ejemplo: hectáreas de maíz en condados de Iowa

5.2 Ajuste del modelo

5.3 Resultados

5.4 Comparación con el estimador directo

6 Estimadores EB de parámetros no lineales basados en el modelo BHF

6.1 Ejemplo: mapeo de pobreza

6.2 Transformación de la variable respuesta

6.3 Definición del parámetro objetivo

6.4 Datos fuera de muestra

6.5 Ajuste del modelo y estimadores EB

6.6 Diagnóstico de residuos

6.7 Estimación del MSE por bootstrap paramétrico

6.8 Resultados

6.9 Comparación con el estimador directo

7 Referencias

8 Información de la sesión

Estimación en áreas pequeñas con el paquete `sae`