De acuerdo con Skinner, Holt, y Smith (1989) y Kish (1987), una cuestión fundamental en el análisis de encuestas es identificar cuál es la medida de probabilidad que rige la inferencia. En la inferencia estadística clásica se suele asumir que las observaciones disponibles son variables aleatorias independientes e idénticamente distribuidas (IID). Bajo ese enfoque, la aleatoriedad proviene del modelo que genera los datos. Sin embargo, este supuesto no describe adecuadamente la información proveniente de encuestas con diseños muestrales complejos, donde las observaciones son seleccionadas mediante estratificación, conglomeración y probabilidades de inclusión desiguales.
Según Särndal, Swensson, y Wretman (2003), en encuestas de hogares conviene distinguir tres objetos. Primero, la población finita, denotada como \(U = \{1,2,\ldots,N\}\), la cual contiene todas las unidades de interés. Segundo, la variable de estudio \(y_k\), observada o definida para cada unidad \(k \in U\). Tercero, el diseño muestral \(p(s)\), que asigna probabilidades a las posibles muestras \(s \subset U\). En la notación utilizada a lo largo de este documento, una unidad puede identificarse como \(k\) dentro del conglomerado o UPM \(i\), perteneciente al estrato \(h\), de modo que los estimadores ponderados suelen escribirse como sumas de la forma
\[ \hat{t}_y = \sum_h \sum_i \sum_k w_{hik} \ y_{hik} \]
donde \(w_{hik}\) representa el factor de expansión o peso muestral asociado con la unidad observada. En la teoría de muestreo, los valores \(y_k\) se consideran fijos una vez definida la población. La aleatoriedad no está en los valores de la variable, sino en los indicadores de inclusión
\[ I_k = \begin{cases} 1, & \text{si } k \in s,\\ 0, & \text{si } k \notin s, \end{cases} \]
que inducen la probabilidad de inclusión \(\pi_k = \Pr_p(I_k=1)\) y, a su vez, el correspondiente peso de muestreo, definido como \(w_k = \frac{1}{\pi_k}\).
Un estimador que ignora los pesos de muestreo puede no representar adecuadamente a la población objetivo, especialmente cuando la muestra proviene de un diseño complejo. En la inferencia basada en el diseño, la población se considera fija y la fuente de aleatoriedad proviene del mecanismo de selección de la muestra; por ello, sus elementos centrales son las probabilidades de inclusión \(\pi_k\) y los pesos de muestreo \(w_k\).
Como demuestran Horvitz y Thompson (1952), por el contrario, la inferencia basada en modelos considera que los valores de la variable de interés en toda la población se conciben como realizaciones de un proceso estocástico, usualmente denotado por \(\xi\). Así, mientras la inferencia basada en los modelos estudia propiedades como \(E_{\xi}(Y_k)\) y \(Var_{\xi}(Y_k)\), la inferencia basada en el diseño se orienta a incorporar las características del plan muestral en los estimadores. La inferencia basada en el diseño no garantiza que todo estimador sea exactamente insesgado; esa propiedad debe verificarse para cada estimador y diseño. En particular, el total de Horvitz–Thompson satisface \(E_{p}(\hat{t}_y) = t_y\) bajo probabilidades de inclusión positivas.
De acuerdo con David A. Binder (2011), esta distinción permite entender por qué, aun cuando exista un modelo razonable para la variable de interés, el diseño muestral debe incorporarse explícitamente si se desea obtener inferencias válidas para la población objetivo. Para comprender la distinción entre inferencia basada en el modelo e inferencia basada en el diseño, partimos de un ejemplo sencillo. Supóngase que se generan \(N=100\) realizaciones independientes de una variable Bernoulli con parámetro \(\theta=0.3\), donde \(\theta\) representa la proporción esperada de personas desempleadas en una población generada por un modelo. Si \(Y_k \sim Bernoulli(\theta)\), entonces la esperanza bajo el modelo es \(E_{\xi}(Y_k) = \theta\), mientras que la varianza bajo el modelo es \(Var_{\xi}(Y_k) = \theta(1-\theta)\).
En la literatura especializada, el modelo \(\xi\) recibe el nombre de modelo de superpoblación, pues genera poblaciones finitas. De esta forma, para una población finita generada por ese modelo, la media poblacional es \(\bar{Y}_U = \frac{1}{N}\sum_{k \in U}Y_k\). Nótese que, bajo el modelo \(\xi\), esta media es insesgada para \(\theta\), puesto que
\[ E_{\xi}(\bar{Y}_U) = E_{\xi}\left(\frac{1}{N}\sum_{k \in U}Y_k\right) = \frac{1}{N}\sum_{k \in U}E_{\xi}(Y_k) = \theta. \]
La siguiente simulación de Monte Carlo reproduce este proceso. En cada repetición se genera una población completa de tamaño \(N=100\) y se calcula su media. El promedio de esas medias poblacionales debe aproximarse a \(\theta\).
set.seed(2026)
population_size <- 100
theta <- 0.3
n_sim_model <- 1000
model_estimates <- rep(NA, n_sim_model)
for (sim in seq_len(n_sim_model)) {
outcome <- rbinom(population_size, 1, theta)
model_estimates[sim] <- mean(outcome)
}
cbind(
theta,
expected_model_mean = mean(model_estimates),
bias = mean(model_estimates) - theta
) theta expected_model_mean bias
[1,] 0.3 0.302 0.00225El código fija una semilla para que el ejercicio sea reproducible. Luego genera n_sim_model poblaciones independientes mediante rbinom(). La salida reporta el parámetro verdadero theta, el promedio de las medias simuladas expected_model_mean y el sesgo Monte Carlo. Como se espera bajo el modelo Bernoulli IID, expected_model_mean queda muy cerca de 0.3 y el sesgo es cercano a cero. Las pequeñas diferencias se deben a que se usan 1000 simulaciones, no un número infinito de repeticiones.
En la teoría de muestreo, en cambio, las características de interés son parámetros fijos de la población. Si una persona está desempleada, su estado se considera un valor fijo de la población finita. Lo aleatorio es el mecanismo mediante el cual se selecciona la muestra. El total poblacional se define como \(t_y = \sum_{k \in U} y_k\), mientras que la media poblacional es \(\bar{y}_U = \frac{t_y}{N}\). Si se observa una muestra \(s\), se propone estimar el total como
\[ \hat{t}_{y} = \sum_{k \in s}\frac{y_k}{\pi_k} = \sum_{k \in s}w_k y_k. \]
Como se mencionó en los capítulos anteriores, cuando el tamaño poblacional \(N\) es desconocido, la media poblacional puede estimarse mediante el estimador de razón ponderado, definido como
\[ \hat{\bar{y}} = \frac{\sum_{k \in s}w_k y_k} {\sum_{k \in s}w_k} \]
Según Cochran (1977), estos estimadores incorporan las probabilidades de inclusión en su forma funcional y esa es la diferencia central frente al promedio simple no ponderado de una muestra. Para ilustrar el impacto del diseño, supóngase que la población anterior se divide en \(N_I\) UPMs. La UPM \(i\) tiene tamaño \(N_i\), total \(t_{y_i}=\sum_{k \in U_i}y_{ik}\) y media \(\bar{y}_i=t_{y_i}/N_i\). Bajo un diseño de tamaño fijo que asigna probabilidades de inclusión de primer orden proporcionales a \(N_i\), siempre que todas sean menores o iguales que uno, estas probabilidades toman la siguiente forma:
\[ \pi_{Ii} = \frac{n_I N_i}{N} \]
Cuando todos los miembros de las UPM seleccionadas son observados, el estimador de Horvitz–Thompson de la media poblacional es exactamente insesgado bajo las probabilidades de inclusión anteriores y puede escribirse como
\[ \hat{\bar{y}} = \frac{1}{N}\sum_{i \in S_I}\frac{t_{y_i}}{\pi_{Ii}} = \frac{1}{N}\sum_{i \in S_I}\frac{N_i\bar{y}_i}{n_I N_i/N} = \frac{1}{n_I}\sum_{i \in S_I}\bar{y}_i. \]
Esta expresión ejemplifica por qué, bajo este diseño específico de selección proporcional al tamaño, el promedio de las medias de los conglomerados brinda insesgamiento. En contraste, el promedio simple de las personas observadas en las UPMs seleccionadas, sería:
\[ \bar{y}_s = \frac{\sum_{i \in S_I}t_{y_i}} {\sum_{i \in S_I}N_i}, \]
el cual trata la muestra como si fuera autoponderada. Este estimador puede ser sesgado cuando las probabilidades de selección son desiguales, especialmente si el tamaño del conglomerado se relaciona con la variable de interés. La siguiente simulación fija una población finita y repite muchas veces la selección de hogares con probabilidades proporcionales al tamaño mediante S.piPS() del paquete TeachingSampling. En cada muestra se calculan dos estimadores: design_estimates, que corresponde a \(\hat{\bar{y}}\), y simple_estimates, que corresponde al promedio simple de las personas observadas en las UPMs seleccionadas, \(\bar{y}_s\). Además, estos estimadores se comparan con el parámetro del modelo theta_population, correspondiente a \(\theta = 0.3\).
library(TeachingSampling)
set.seed(2026)
population_size <- 100
theta <- 0.3
n_sim_sampling <- 1000
design_estimates <- simple_estimates <- rep(NA, n_sim_sampling)
cluster_size <- rep(2:6, each = 5)
n_clusters <- length(cluster_size)
cluster_id <- rep(seq_len(n_clusters), cluster_size)
size_center <- weighted.mean(cluster_size, cluster_size)
prob_person <- theta + rep((cluster_size - size_center) / 12, cluster_size)
outcome <- rbinom(population_size, 1, prob_person)
theta_population <- mean(outcome)
cluster_totals <- tapply(outcome, cluster_id, sum)
cluster_means <- tapply(outcome, cluster_id, mean)
sampled_cluster_count <- floor(n_clusters * 0.3)
for (sim in seq_len(n_sim_sampling)) {
pps_sample <- S.piPS(sampled_cluster_count, cluster_size)
sampled_clusters <- pps_sample[, 1]
sampled_cluster_size <- cluster_size[sampled_clusters]
sampled_cluster_total <- cluster_totals[sampled_clusters]
sampled_cluster_mean <- cluster_means[sampled_clusters]
design_estimates[sim] <- mean(sampled_cluster_mean)
simple_estimates[sim] <-
sum(sampled_cluster_total) / sum(sampled_cluster_size)
}
cbind(
theta_population,
expected_pps = mean(design_estimates),
bias_pps = mean(design_estimates) - theta_population,
expected_simple = mean(simple_estimates),
bias_simple = mean(simple_estimates) - theta_population
) theta_population expected_pps bias_pps expected_simple bias_simple
[1,] 0.3 0.304 0.00368 0.358 0.0585El primer bloque del código genera hogares de tamaños desiguales, entre 2 y 6 personas. Para hacer visible el efecto del diseño, la probabilidad individual de desempleo varía suavemente con el tamaño del hogar, pero se centra de manera que el promedio esperado de la población siga siendo cercano a theta. Luego se calcula theta_population, que es la proporción real de desempleo en esa población específica. En cada repetición, S.piPS() selecciona UPMs con probabilidades proporcionales al tamaño; después se calculan la media de las medias de las UPMs seleccionadas y el promedio simple de las personas observadas. La salida compara ambos estimadores con theta_population. El sesgo de Monte Carlo de expected_pps es cercano a cero, mientras que el de expected_simple no lo es, mostrando resultado de ignorar el diseño en este ejercicio.
No obstante, en muchas aplicaciones ambos marcos inferenciales no se presentan como enfoques excluyentes, sino como componentes complementarios de un mismo problema. En particular, puede considerarse que la muestra ha sido obtenida mediante un diseño probabilístico complejo y que, al mismo tiempo, los valores de la variable de interés responden a un modelo de superpoblación. David A. Binder (2011) plantea esta formulación como un esquema de inferencia doble, en el cual la incertidumbre no proviene de una única fuente, sino de la combinación entre el mecanismo que genera los valores poblacionales (modelo de superpoblación) y el mecanismo que determina qué unidades son observadas (diseño de muestreo).
Formalmente, intervienen dos medidas de probabilidad: \(\xi\), asociada al proceso generador de los valores \(Y_k\), y \(p(\cdot)\), asociada al diseño muestral que selecciona la muestra \(s\). Por tanto, las propiedades de los estimadores deben evaluarse considerando simultáneamente ambas fuentes de variación. Bajo este marco, la esperanza combinada puede expresarse como
\[ E_{\xi p}(\hat{\theta}) = E_{\xi}\{E_p(\hat{\theta}\mid U)\}, \]
Donde primero se evalúa el comportamiento del estimador bajo el diseño, condicionado a la población finita generada, y luego se promedia sobre el modelo de superpoblación. En este contexto, la verosimilitud conjunta completa puede resultar intratable, especialmente cuando las probabilidades de inclusión dependen de variables relacionadas con la respuesta. Las ecuaciones de estimación ponderadas de David A. Binder (1983) fundamentan la máxima pseudo-verosimilitud (MPV), que pondera las contribuciones de cada unidad por el inverso de su probabilidad de inclusión; Pfeffermann (1993) examina las condiciones y los objetivos bajo los cuales resulta pertinente incorporar esos pesos en modelos para datos de encuestas.
El método de máxima verosimilitud, presentado formalmente por Casella y Berger (2002), parte de una muestra aleatoria \(y_1,\ldots,y_n\) generada por una distribución con densidad o función de probabilidad \(f(y;\theta)\). Si las observaciones son independientes e idénticamente distribuidas (IID), la función de verosimilitud se define como
\[ L(\theta) = \prod_{k=1}^{n} f(y_k;\theta). \]
El límite superior correcto es el tamaño muestral \(n\), no el tamaño poblacional \(N\). Como los productos pueden ser difíciles de manipular, se utiliza la log-verosimilitud:
\[ \ell(\theta) = \sum_{k=1}^{n}\log f(y_k;\theta) \]
El estimador de máxima verosimilitud \(\hat{\theta}\) es el valor de \(\theta\) que maximiza \(\ell(\theta)\). Si la función es diferenciable, este valor satisface las ecuaciones de score:
\[ U(\theta) = \frac{\partial \ \ell(\theta)}{\partial \ \theta} = \sum_{k=1}^{N} u_k(\theta) = 0 \]
donde \(u_k(\theta) = \frac{\partial}{\partial\theta}\log f(y_k;\theta)\) es la contribución de la unidad \(k\) al score total. Para una distribución Bernoulli con parámetro \(\theta\), la función de probabilidad es \(f(y_k;\theta) = \theta^{y_k}(1-\theta)^{1-y_k}\), para \(y_k \in \{0,1\}\). La log-verosimilitud es
\[ \ell(\theta) = \sum_{k=1}^{N} \left[ y_k\log(\theta) +(1-y_k)\log(1-\theta) \right] \]
Al derivar e igualar a cero y despejar, se obtiene el estimador MV, definido por:
\[ \hat{\theta}_{MV} = \frac{1}{n}\sum_{k=1}^{n}y_k \]
Así, la proporción muestral es el estimador natural del parámetro Bernoulli. El punto clave es que este resultado depende del supuesto de que las observaciones tienen igual distribución y no provienen de un diseño con probabilidades desiguales.
Por otro lado, en un modelo de regresión lineal múltiple con errores normales, se tiene que \(\mathbf{y} \sim N(\mathbf{X}\boldsymbol{\beta},\sigma^2 I)\). En este caso, la log-verosimilitud, omitiendo constantes, puede escribirse como
\[ \ell(\boldsymbol{\beta},\sigma^2) = -\frac{n}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})'(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \]
Maximizar esta expresión respecto de \(\beta\) equivale a minimizar la suma de cuadrados de los residuos. Por ello, el estimador de máxima verosimilitud de \(\beta\) coincide con el estimador de mínimos cuadrados ordinarios \(\hat{\boldsymbol{\beta}}_{MV} = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{y}\). Este resultado es importante porque muestra que muchos procedimientos estadísticos usuales pueden interpretarse como máxima verosimilitud bajo supuestos de modelo específicos. Sin embargo, si los datos provienen de una encuesta compleja, la matriz \(\mathbf{X}'\mathbf{X}\) no refleja por sí sola el mecanismo de selección de la muestra.
Cuando el diseño muestral es complejo, los datos observados no satisfacen necesariamente los supuestos de independencia e idéntica distribución. Algunas unidades pueden representar a muchas personas de la población y otras a pocas. Además, dos unidades seleccionadas dentro del mismo conglomerado pueden estar correlacionadas. En este contexto, usar directamente la log-verosimilitud ordinaria puede producir estimaciones que describen la muestra, pero no la población objetivo.
La máxima pseudo-verosimilitud, desarrollada a partir de los resultados de David A. Binder (1983), incorpora los pesos de muestreo en la log-verosimilitud. Para unidades observadas \(k \in s\), se define
\[ \ell_p(\theta) = \sum_{k \in s} w_k \log f(y_k;\theta) \]
En la notación de esta sección, el peso \(w_k=1/\pi_k\) indica cuántas unidades de la población representa la observación \(y_k\). Por eso, la contribución de una unidad con baja probabilidad de inclusión se amplifica en la pseudo-verosimilitud. Las ecuaciones de estimación resultantes son
\[ U_p(\theta) = \sum_{k \in s} w_k \ u_k(\theta) = 0 \]
La estimación MPV conserva la lógica de máxima verosimilitud, pero reemplaza el score ordinario por un score ponderado. La solución \(\hat{\theta}_{MPV}\) es el valor que hace que la suma ponderada de contribuciones individuales sea igual a cero. Para una distribución Bernoulli, la pseudo-log-verosimilitud es
\[ \ell_p(\theta) = \sum_{k \in s} w_k \left[ y_k\log(\theta) +(1-y_k)\log(1-\theta) \right]. \]
Al derivar e igualar a cero, se obtiene que:
\[ \frac{\partial \ell_p(\theta)}{\partial \theta} = \sum_{k \in s} w_k \left[ \frac{y_k}{\theta} -\frac{1-y_k}{1-\theta} \right] =0 \]
Cuya solución es:
\[ \hat{\theta}_{MPV} = \frac{\sum_{k \in s}w_k y_k} {\sum_{k \in s}w_k} = \hat{p}_d \]
Por tanto, para una variable binaria, la MPV conduce al estimador ponderado de la proporción, equivalente al estimador de Hájek. Este resultado conecta directamente la teoría de pseudo-verosimilitud con los estimadores de proporciones presentados en los capítulos de variables categóricas.
Asimismo, para un modelo de regresión lineal múltiple, la pseudo-log-verosimilitud ponderada implica minimizar una suma de cuadrados ponderada:
\[ (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})'\mathbf{W}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}), \]
donde
\[ \mathbf{W}= \begin{pmatrix} w_1 & 0 & \cdots & 0\\ 0 & w_2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & w_n \end{pmatrix}. \]
La solución es el estimador de mínimos cuadrados ponderados:
\[ \hat{\boldsymbol{\beta}}_{MPV} = (\mathbf{X}'\mathbf{W}\mathbf{X})^{-1}\mathbf{X}'\mathbf{W}\mathbf{y} \]
De acuerdo con Heeringa, West, y Berglund (2017), este resultado extiende de manera natural el estimador ordinario al contexto de encuestas complejas: cada observación contribuye al ajuste del modelo en proporción a su peso muestral. No obstante, para la inferencia no basta con ponderar la estimación puntual; también se requiere estimar correctamente la varianza considerando estratos, conglomerados y pesos.
La siguiente simulación de Monte Carlo anidada ilustra la inferencia doble. En el ciclo externo se genera una nueva población bajo el modelo Bernoulli. En el ciclo interno se seleccionan muchas muestras por diseño y se comparan el estimador que incorpora el diseño y el promedio simple de la muestra.
set.seed(2026)
population_size <- 100
theta <- 0.3
n_sim_population <- 100
expected_design_estimates <- expected_simple_estimates <- population_means <-
rep(NA, n_sim_population)
cluster_size <- rep(2:6, each = 5)
n_clusters <- length(cluster_size)
cluster_id <- rep(seq_len(n_clusters), cluster_size)
size_center <- weighted.mean(cluster_size, cluster_size)
prob_person <- theta + rep((cluster_size - size_center) / 12, cluster_size)
sampled_cluster_count <- floor(n_clusters * 0.3)
for (population_sim in seq_len(n_sim_population)) {
outcome <- rbinom(population_size, 1, prob_person)
population_means[population_sim] <- mean(outcome)
cluster_totals <- tapply(outcome, cluster_id, sum)
cluster_means <- tapply(outcome, cluster_id, mean)
design_estimates <- simple_estimates <- rep(NA, 100)
for (sample_sim in seq_len(100)) {
pps_sample <- S.piPS(sampled_cluster_count, cluster_size)
sampled_clusters <- pps_sample[, 1]
sampled_cluster_size <- cluster_size[sampled_clusters]
sampled_cluster_total <- cluster_totals[sampled_clusters]
sampled_cluster_mean <- cluster_means[sampled_clusters]
design_estimates[sample_sim] <- mean(sampled_cluster_mean)
simple_estimates[sample_sim] <-
sum(sampled_cluster_total) / sum(sampled_cluster_size)
}
expected_design_estimates[population_sim] <- mean(design_estimates)
expected_simple_estimates[population_sim] <- mean(simple_estimates)
}
cbind(
theta,
mean_population = mean(population_means),
expected_pps = mean(expected_design_estimates),
bias_pps = mean(expected_design_estimates) - theta,
expected_simple = mean(expected_simple_estimates),
bias_simple = mean(expected_simple_estimates) - theta
) theta mean_population expected_pps bias_pps expected_simple bias_simple
[1,] 0.3 0.304 0.302 0.00248 0.332 0.0317El código reproduce las dos fuentes de aleatoriedad. El ciclo externo representa el modelo \(\xi\), porque genera nuevas poblaciones finitas con rbinom(). El ciclo interno representa el diseño \(p\), porque para cada población se seleccionan múltiples muestras con S.piPS(). La salida incluye mean_population, que es el promedio de las medias poblacionales generadas; expected_pps, que resume el comportamiento del estimador que incorpora el diseño; y expected_simple, que resume el promedio simple de la muestra seleccionada. Cuando el diseño se incorpora correctamente, el sesgo Monte Carlo respecto de theta debe ser pequeño. La comparación con expected_simple muestra por qué la inferencia basada en encuestas requiere pesos y varianzas de diseño, incluso cuando se parte de un modelo estadístico sencillo.
Como explica Pfeffermann (1993), en síntesis, la máxima verosimilitud ordinaria es adecuada cuando las observaciones pueden tratarse como IID bajo un modelo probabilístico. En encuestas complejas, la máxima pseudo-verosimilitud ofrece una extensión natural: conserva la estructura de las ecuaciones de estimación, pero pondera cada contribución individual por su peso muestral. La ponderación protege frente a ciertos diseños informativos o a determinadas especificaciones incorrectas, pero su necesidad y sus propiedades dependen del diseño y del modelo. Así, los modelos estadísticos se ajustan de manera más coherente con la población que la encuesta busca representar.