5.6 Estimación de los parámetros en un modelo de regresión con muestras complejas
Una vez se establecen los supuestos del modelo y las características distribucionales de los errores, el paso siguiente es el proceso de estimación de los parámetros. A modo ilustrativo, si en lugar de observar una muestra de tamaño \(n\) de los \(N\) elementos de la población se hubiera realizado un censo completo, el parámetro de regresión de población finita \(\beta_{1}\) podría calcularse como sigue (Téllez et al., 2016):
\[ \beta_{1} = \frac{\sum_{i=1}^{N}\left(X_{i}-\bar{X}\right)\left(Y_{i}-\bar{Y}\right)}{\sum_{i=1}^{N}\left(X_{i}-\bar{X}\right)^{2}} \]
Sin embargo, cuando se desea estimar los parámetros de un modelo de regresión lineal utilizando información proveniente de encuestas con diseños complejos, el enfoque estándar cambia. La principal razón es que los datos no cumplen con los supuestos clásicos de independencia e idéntica distribución, ya que el diseño muestral introduce elementos como estratificación, conglomerados o probabilidades de selección desiguales. En consecuencia, los estimadores tradicionales (como los obtenidos por máxima verosimilitud o mínimos cuadrados) pueden resultar sesgados y sus errores estándar poco confiables.
En este contexto, Wolter (2007) propone el uso de métodos no paramétricos robustos, como la linealización de Taylor o los procedimientos de replicación (Jackknife, Bootstrap, Balanced Repeated Replication, entre otros), para estimar varianzas y obtener inferencias válidas.
5.6.1 Estimación de parámetros
En el caso de un modelo de regresión lineal simple, la estimación de la pendiente \((\beta_1)\) bajo un esquema de encuesta compleja se realiza mediante un estimador ponderado, que puede expresarse de la siguiente forma:
\[ \hat{\beta_{1}} = \frac{\sum_{h=1}^H \sum_{\alpha=1}^{a_h} \sum_{i=1}^{n_{h\alpha}} \omega_{h\alpha i}\,(y_{h\alpha i}-\bar{y}_{\omega})(x_{h\alpha i}-\bar{x}_{\omega})} {\sum_{h=1}^H \sum_{\alpha=1}^{a_h} \sum_{i=1}^{n_{h\alpha}} \omega_{h\alpha i}\,(x_{h\alpha i}-\bar{x}_{\omega})^{2}} = \frac{t_{xy}}{t_{x^{2}}} \]
donde \(\omega_{h\alpha i}\) representa los pesos de muestreo. La principal diferencia respecto al estimador clásico es la incorporación explícita de dichos pesos, que corrigen las probabilidades desiguales de selección y aseguran la representatividad de las estimaciones.
La varianza del estimador \(\hat{\beta_1}\) puede expresarse como:
\[ var\left(\hat{\beta_{1}}\right) = \frac{var(t_{xy})+\hat{\beta}_{1}^{2}var(t_{x^{2}})-2\hat{\beta}_{1}cov(t_{xy},t_{x^{2}})}{(t_{x^{2}})^{2}} \]
Este enfoque permite cuantificar la precisión del coeficiente considerando tanto los pesos como la estructura del diseño muestral.
En el caso de la regresión múltiple, el cálculo de la varianza de cada coeficiente se realiza considerando su interdependencia con los demás. Esto se refleja en la construcción de una matriz de varianza-covarianza, que recoge tanto la variabilidad individual como las covarianzas entre todos los parámetros. De acuerdo con Kish y Frankel (1974), esta estimación requiere utilizar totales ponderados de cuadrados y productos cruzados de todas las combinaciones entre la variable dependiente \(y\) y el conjunto de predictores \(x = {1, x_1, …, x_p}\). En términos generales:
\[ var(\hat{\beta}) = \hat{\Sigma}(\hat{\beta}) = \begin{bmatrix} var(\hat{\beta}_{0}) & cov(\hat{\beta}_{0},\hat{\beta}_{1}) & \cdots & cov(\hat{\beta}_{0},\hat{\beta}_{p}) \\ cov(\hat{\beta}_{0},\hat{\beta}_{1}) & var(\hat{\beta}_{1}) & \cdots & cov(\hat{\beta}_{1},\hat{\beta}_{p}) \\ \vdots & \vdots & \ddots & \vdots \\ cov(\hat{\beta}_{0},\hat{\beta}_{p}) & cov(\hat{\beta}_{1},\hat{\beta}_{p}) & \cdots & var(\hat{\beta}_{p}) \end{bmatrix} \]
Este marco permite evaluar la estabilidad y precisión de los parámetros en modelos más complejos, garantizando que las inferencias reflejen adecuadamente tanto el diseño de la encuesta como la relación entre las variables.
Para ejemplificar los conceptos trabajados hasta este momento, se tomará la misma base que se ha venido trabajando durante todo el desarrollo de este libro. Se inicia con el cargue de las librerías, la base de datos y la definición del diseño de muestreo:
knitr::opts_chunk$set(warning = FALSE, message = FALSE, error = FALSE)
options(digits = 4)
options(tinytex.verbose = TRUE)
library (survey)
library(srvyr)
library(convey)
library(TeachingSampling)
library(printr)
library(stargazer)
library(jtools)
library(broom)Cargue de la base y definición del diseño muestral:
data(BigCity, package = "TeachingSampling")
library(tidyverse)
encuesta <- readRDS("Data/encuesta.rds")
head(encuesta)| HHID | Stratum | NIh | nIh | dI | PersonID | PSU | Zone | Sex | Age | MaritalST | Income | Expenditure | Employment | Poverty | dki | dk | wk | Region | CatAge |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| idHH00031 | idStrt001 | 9 | 2 | 4.5 | idPer01 | PSU0003 | Rural | Male | 68 | Married | 409.9 | 346.3 | Employed | NotPoor | 8 | 36 | 34.50 | Norte | Más de 60 |
| idHH00031 | idStrt001 | 9 | 2 | 4.5 | idPer02 | PSU0003 | Rural | Female | 56 | Married | 409.9 | 346.3 | Employed | NotPoor | 8 | 36 | 33.64 | Norte | 46-60 |
| idHH00031 | idStrt001 | 9 | 2 | 4.5 | idPer03 | PSU0003 | Rural | Female | 24 | Married | 409.9 | 346.3 | Employed | NotPoor | 8 | 36 | 33.64 | Norte | 16-30 |
| idHH00031 | idStrt001 | 9 | 2 | 4.5 | idPer04 | PSU0003 | Rural | Male | 26 | Married | 409.9 | 346.3 | Employed | NotPoor | 8 | 36 | 34.50 | Norte | 16-30 |
| idHH00031 | idStrt001 | 9 | 2 | 4.5 | idPer05 | PSU0003 | Rural | Female | 3 | NA | 409.9 | 346.3 | NA | NotPoor | 8 | 36 | 33.64 | Norte | 0-5 |
| idHH00041 | idStrt001 | 9 | 2 | 4.5 | idPer01 | PSU0003 | Rural | Female | 61 | Widowed | 823.8 | 392.2 | Employed | NotPoor | 8 | 36 | 33.64 | Norte | Más de 60 |
library(srvyr)
diseno <- encuesta %>%
as_survey_design(
strata = Stratum,
ids = PSU,
weights = wk,
nest = T
)Para efectos de los ejemplos y como se ha hecho en anteriores ocasiones, se divide la muestra en sub-grupos de la encuesta como sigue:
sub_Urbano <- diseno %>% filter(Zone == "Urban")
sub_Rural <- diseno %>% filter(Zone == "Rural")
sub_Mujer <- diseno %>% filter(Sex == "Female")
sub_Hombre <- diseno %>% filter(Sex == "Male")En este capítulo se ajustarán los modelos de regresión usando la base de datos de ejemplo que se ha venido trabajando en capítulos anteriores. Puesto que, en modelos de regresión, se utiliza muy frecuente el recurso gráfico. A continuación, se define un tema estándar que la CEPAL tiene para generar sus gráficos el cual se utilizará en este capítulo.
Para observar que existe una correlación entre el ingreso y el gasto, las cuales son las variables que se utilizarán para el ajuste de los modelos, se construye un scatterplot usando la librería ggplot. Cabe resaltar que, como la base de datos encuesta, la cual se usa para ejemplificar es una muestra de la base BigCity, analizaremos de manera gráfica si poblacionalmente las dos variables mencionadas anteriormente tienen correlación como se muestra a continuación:
library(ggplot2); library(ggpmisc)
plot_BigCity <- ggplot(data = BigCity,
aes(x = Expenditure, y = Income)) +
geom_point() + geom_smooth(method = "lm",
se = FALSE,
formula = y ~ x) + theme_cepal()
plot_BigCity + stat_poly_eq(formula = y~x, aes(label = paste(..eq.label..,
..rr.label.., sep = "~~~"),size = 3), parse = TRUE)
Si bien, existen unas observaciones por fuera de la nube de punto, el comportamiento general de la relación ingresos vs gastos mantiene una tendencia lineal.
Una vez hecho el análisis gráfico de las variables a utilizar en los modelos a trabajar, se realizará primero un ajuste del modelo con los datos poblacionales y con esto poder analizar qué tan bueno serán los ajustes que se realizarán posteriormente. A continuación, se muestra el ajuste del modelo con los datos poblacionales:
Ahora bien, para observar los parámetros poblacionales del modelo se utilizará la función modelsummary de la librería modelsummary de la siguiente manera:
library(knitr)
tabla <- data.frame(
Term = c("(Intercept)", "Expenditure", "Num.Obs.", "R2", "R2 Adj.", "RMSE"),
Pob = c(123.337, 1.229, 150266, 0.359, 0.359, 461.74)
)
kable(tabla,
format = "latex",
booktabs = TRUE,
digits = 3,
col.names = c("", "Pob"))De la anterior salida se puede observar que, el intercepto es igual a 123.337 y el parámetro \(\beta_{1}\) asociado al gasto es 1.229. La demás información relacionada a esta salida se analizará más adelante.
Una vez revisada la información poblacional, se utilizará la información obtenida de la muestra para estimar los parámetros y con ello analizar qué tan buenas son las estimaciones. A continuación, se presenta una sintaxis similar a la anterior que permite construir el scatterplot pero para los datos de la muestra.
plot_sin <- ggplot(data = encuesta,
aes(x = Expenditure, y = Income)) +
geom_point() +
geom_smooth(method = "lm",
se = FALSE, formula = y ~ x) + theme_cepal()
plot_sin + stat_poly_eq(formula = y~x, aes(label = paste(..eq.label..,
..rr.label.., sep = "~~~"), size = 5), parse = TRUE)
Como se puede observar, los datos de la muestra tienen una tendencia lineal aunque un poco dispersa a medida que crecen los gastos en las familias.
Una vez hecho los análisis gráficos se procede a ajustar los modelos de regresión lineal. A modo de comparar el efecto que tiene hacer un correcto uso de los factores de expansión del diseño, primero, se ajustará un modelo sin tener encuesta dichos factores como se muestra a continuación:
fit_sinP <- lm(Income ~ Expenditure, data = encuesta)
stargazer(fit_sinP, type = "latex",
title = "Modelo sin factores de expansion",
label = "tab:sin_factores",
header = FALSE)Para el modelo ajustado sin factores de expansión, el \(\hat{\beta}_{0}\) es 121.52 y el \(\hat{\beta}_{1}\) asociado a la variable gastos es 1.22.
Ahora, haciendo un Scatterplot con los datos encuesta pero utilizando los factores de expansión del diseño se debe agregar mapping = aes(weight = wk) en la función geom_smoothcomo sigue:
plot_Ponde <- ggplot(data = encuesta,
aes(x = Expenditure, y = Income)) +
geom_point(aes(size = wk)) +
geom_smooth(method = "lm", se = FALSE, formula = y ~ x, mapping = aes(weight = wk)) + theme_cepal()
plot_Ponde + stat_poly_eq(formula = y~x, aes(weight = wk,
label = paste(..eq.label..,..rr.label.., sep = "~~~")), parse = TRUE,size = 5)
En este sentido, para ajustar modelos teniendo en cuenta los factores de expansión existen 2 formas, la primera es usando la función lm y la segunda es usando la función svyglm de la librería survey. A continuación. se ajusta el modelo usando la función lm:
fit_Ponde <- lm(Income ~ Expenditure, data = encuesta, weights = wk)
stargazer(fit_Ponde, header = FALSE,
title = "Modelo encuesta ponderada",
type = "latex")Para el modelo ajustado con factores de expansión usando la función lm, el \(\hat{\beta}_{0}\) es 103.14 y el \(\hat{\beta}_{1}\) asociado a la variable gastos es 1.26. Ahora, haciendo el mismo ajuste pero usando la función svyglm:
## Stratified 1 - level Cluster Sampling design (with replacement)
## With (238) clusters.
## Called via srvyr
## Sampling variables:
## - ids: PSU
## - strata: Stratum
## - weights: wk
##
## Call: svyglm(formula = Income ~ Expenditure, design = diseno, family = stats::gaussian())
##
## Coefficients:
## (Intercept) Expenditure
## 103.14 1.26
##
## Degrees of Freedom: 2604 Total (i.e. Null); 118 Residual
## Null Deviance: 6.35e+08
## Residual Deviance: 3.11e+08 AIC: 38300Obteniendo estimaciones para el \(\hat{\beta}_{0}\) es 103.14 y el \(\hat{\beta}_{1}\) asociado a la variable gastos es 1.26. Siendo exactamente las mismas que con la función lm ya que, como se definió en los argumentos de la función, la función de enlace es Gausiana.
Por último y a modo de resumen se muestra un gráfico donde se encuentran depositados todos los modelos estimados anteriormente y así poder comparar de manera gráfica su ajuste:
df_model <- data.frame(
intercept = c(coefficients(fit)[1],
coefficients(fit_sinP)[1],
coefficients(fit_Ponde)[1],
coefficients(fit_svy)[1]),
slope = c(coefficients(fit)[2],
coefficients(fit_sinP)[2],
coefficients(fit_Ponde)[2],
coefficients(fit_svy)[2]),
Modelo = c("Población", "Sin ponderar",
"Ponderado(lm)", "Ponderado(svyglm)"))
plot_BigCity + geom_abline( data = df_model,
mapping = aes( slope = slope,
intercept = intercept, linetype = Modelo,
color = Modelo ), size = 2
)
5.6.2 Uso de ponderaciones
En el análisis de datos provenientes de encuestas con diseños complejos surge un interrogante clave: ¿cómo deben incorporarse los pesos muestrales en los modelos de regresión? La respuesta no es trivial, ya que los pesos reflejan tanto la probabilidad de selección como los ajustes por no respuesta y calibración, pero su uso directo en el modelado puede generar problemas de eficiencia y aumentar la varianza de los estimadores.
Para enfrentar estas limitaciones, se han propuesto diferentes estrategias de ajuste que buscan lograr un balance entre precisión y eficiencia en la estimación. Entre los procedimientos más utilizados destacan los siguientes:
Pesos tipo Senado Este procedimiento ajusta los pesos de manera que su suma coincida con el tamaño de la muestra, en lugar del tamaño de la población. El objetivo es mantener la representatividad relativa de las unidades, pero reduciendo la dispersión de los pesos originales, lo que resulta particularmente ventajoso en encuestas donde existe alta variabilidad entre los factores de expansión:
\[ w_k^{Senate} = w_k \times \frac{n}{\sum w_k} \]
Pesos normalizados En este enfoque los pesos originales se reescalan para que su suma sea igual a uno, lo que evita un incremento innecesario de la varianza en los modelos. Esta técnica es especialmente útil cuando se trabaja con diferentes subconjuntos de datos (por ejemplo, modelos estimados en subpoblaciones) o cuando se busca minimizar la inflación de la varianza:
\[ w_k^{Normalized} = \frac{w_k}{\sum w_k} \]
Es importante resaltar que, en ambos métodos, los pesos ajustados se obtienen mediante transformaciones multiplicativas directas de los pesos muestrales originales. Por esta razón, no deben emplearse para calcular totales poblacionales, ni modifican los coeficientes de variación asociados a ellos. Asimismo, estos procedimientos no afectan las estimaciones de razones, como medias o proporciones, ya que en dichos casos los pesos se cancelan en el cociente.
En la práctica, el uso de estos ajustes puede considerarse una solución pragmática en contextos donde no se dispone de software especializado para encuestas. Sin embargo, cuando se cuenta con programas estadísticos que permiten la incorporación directa de pesos y del diseño muestral —como los descritos en la Subsección 9.3.4— el reescalamiento de los pesos deja de ser necesario, ya que dichos programas realizan el tratamiento adecuado para preservar tanto la representatividad como las propiedades inferenciales del modelo.