6.1 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 siguientes es el proceso de estimación de los parámetros. A modo ilustrativo e introductorio al proceso de estimación de los parámetros con información provenientes de muestras complejas, si en lugar de observar una muestra de tamaño \(n\) de los \(N\) elementos de 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 (Tellez, 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}} \]
Ahora bien, cuando se desea estimar los parámetros de un modelo de regresión lineal, pero considerando que la información muestral proviene de encuestas con muestras complejas se altera el enfoque estándar que se le da a la estimación de coeficientes de regresión y sus errores estándar.
La principal razón por la que los métodos de estimación de parámetros de los coeficientes de regresión cambian es que la información recolectada por medio de una encuesta con muestra compleja generalmente no está distribuida de manera idéntica dado que el diseño muestral así es planeado. En este contexto, al ajustar modelos de regresión con este conjunto de datos, dado que los diseños complejos en su mayoría contienen estratificación, conglomerados, probabilidades de selección desiguales, etc, impiden el uso de estimadores de varianza convencionales que se pueden derivar por máxima verosimilitud puesto que, con esta metodología se asumen que los datos son independientes e idénticamente distribuidos y que provienen de alguna distribución de probabilidad (binomial, Poisson, exponencial, normal, etc.). En su lugar, según Wolter, (2007) se emplean métodos no paramétricos robustos basados en linealización de Taylor o métodos de estimación de la varianza usando replicación (Jackknife, bootstrapping, etc).
Con el fin de ilustrar la forma cómo se estiman los parámetros de regresión del modelo en el contexto de encuestas complejas se realizará estimando el parámetro \(\beta_{1}\) y su varianza para una regresión lineal simple. La extensión a la estimación de los parámetros de un modelo de regresión múltiple, algebraicamente es compleja y se sale del contexto de este libro. A continuación, se presenta la estimación del intercepto y su varianza en un modelo de regresión lineal simple:
\[ \hat{\beta_{1}} = \frac{{\sum_{h}^{H}\sum_{\alpha}^{a_{h}}\sum_{i=1}^{n_{h\alpha}}\omega_{h\alpha i}\left(y_{h\alpha i}-\bar{y}_{\omega}\right)\left(x_{h\alpha i}-\bar{x}_{\omega}\right)}}{{ \sum_{h}^{H}\sum_{\alpha}^{a_{h}}\sum_{i=1}^{n_{h\alpha}}\omega_{h\alpha i}\left(x_{h\alpha i}-\bar{x}_{\omega}\right)^{2}}} = \frac{t_{xy}}{t_{x^{2}}} \]
Como se puede observar en la ecuación anterior, el estimador del parámetro es un cociente de totales, por ende, su varianza estimada está dada por:
\[ var\left(\hat{\beta_{1}}\right) = \frac{var\left(t_{xy}\right)+\hat{\beta}_{1}^{2}var\left(t_{x^{2}}\right)-2\hat{\beta}_{1}cov\left(t_{xy},t_{x^{2}}\right)}{\left(t_{x^{2}}\right)^{2}} \]
A modo de generalización según Kish y Frankel, (1974) para la estimación de la varianza en un modelo de regresión lineal múltiple los métodos de aproximación requieren totales de muestra ponderados para los cuadrados y productos cruzados de todas las combinaciones \(y\) y \(x = {1 x_{1} … x_{p}}\). A continuación, se presenta la estimación:
\[\begin{eqnarray*} var\left(\hat{\beta}\right)=\hat{\Sigma}\left(\hat{\beta}\right) & = & \left[\begin{array}{cccc} var\left(\hat{\beta}_{0}\right) & cov\left(\hat{\beta}_{0},\hat{\beta}_{1}\right) & \cdots & cov\left(\hat{\beta}_{0},\hat{\beta}_{p}\right)\\ cov\left(\hat{\beta}_{0},\hat{\beta}_{1}\right) & var\left(\hat{\beta}_{1}\right) & \cdots & cov\left(\hat{\beta}_{1},\hat{\beta}_{p}\right)\\ \vdots & \vdots & \ddots & \vdots\\ cov\left(\hat{\beta}_{0},\hat{\beta}_{p}\right) & cov\left(\hat{\beta}_{1},\hat{\beta}_{p}\right) & \cdots & var\left(\hat{\beta}_{p}\right) \end{array}\right] \end{eqnarray*}\]
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:
fit <- lm(Income ~ Expenditure, data = BigCity)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:
| Pob | |
|---|---|
| (Intercept) | 123.337 |
| Expenditure | 1.229 |
| Num.Obs. | 150266 |
| R2 | 0.359 |
| R2 Adj. | 0.359 |
| RMSE | 461.74 |
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)
fit_sinP##
## Call:
## lm(formula = Income ~ Expenditure, data = encuesta)
##
## Coefficients:
## (Intercept) Expenditure
## 121.52 1.22stargazer(fit_sinP, header = FALSE,
title = "Modelo sin factores de expansion",
style = "ajps")##
## \begin{table}[!htbp] \centering
## \caption{Modelo sin factores de expansion}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lc}
## \\[-1.8ex]\hline \\[-1.8ex]
## \\[-1.8ex] & \textbf{Income} \\
## \hline \\[-1.8ex]
## Expenditure & 1.220$^{***}$ \\
## & (0.025) \\
## Constant & 121.500$^{***}$ \\
## & (11.410) \\
## N & 2605 \\
## R-squared & 0.487 \\
## Adj. R-squared & 0.487 \\
## Residual Std. Error & 345.000 (df = 2603) \\
## F Statistic & 2473.000$^{***}$ (df = 1; 2603) \\
## \hline \\[-1.8ex]
## \multicolumn{2}{l}{$^{***}$p $<$ .01; $^{**}$p $<$ .05; $^{*}$p $<$ .1} \\
## \end{tabular}
## \end{table}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)
fit_Ponde##
## Call:
## lm(formula = Income ~ Expenditure, data = encuesta, weights = wk)
##
## Coefficients:
## (Intercept) Expenditure
## 103.14 1.26stargazer(fit_Ponde, header = FALSE,
title = "Modelo encuesta ponderada",
style = "ajps")##
## \begin{table}[!htbp] \centering
## \caption{Modelo encuesta ponderada}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lc}
## \\[-1.8ex]\hline \\[-1.8ex]
## \\[-1.8ex] & \textbf{Income} \\
## \hline \\[-1.8ex]
## Expenditure & 1.263$^{***}$ \\
## & (0.024) \\
## Constant & 103.100$^{***}$ \\
## & (11.260) \\
## N & 2605 \\
## R-squared & 0.509 \\
## Adj. R-squared & 0.509 \\
## Residual Std. Error & 2627.000 (df = 2603) \\
## F Statistic & 2703.000$^{***}$ (df = 1; 2603) \\
## \hline \\[-1.8ex]
## \multicolumn{2}{l}{$^{***}$p $<$ .01; $^{**}$p $<$ .05; $^{*}$p $<$ .1} \\
## \end{tabular}
## \end{table}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:
fit_svy <- svyglm(Income ~ Expenditure,
design = diseno, family=stats::gaussian())
fit_svy## 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
)