11.7 Introducción a la imputación múltiple.

La imputación múltiple, también conocida como “multiple imputations” en inglés, es una técnica estadística utilizada para tratar los datos faltantes o incompletos en un conjunto de datos. La imputación múltiple consiste en crear múltiples copias del conjunto de datos, donde los valores faltantes en cada copia son imputados utilizando modelos estadísticos. Estos modelos se basan en las relaciones entre las variables en el conjunto de datos y se utilizan para estimar los valores faltantes de manera plausible.

Luego de que se han creado múltiples copias completas del conjunto de datos, se realizan análisis separados en cada copia para generar resultados. Los resultados de cada análisis se combinan para obtener un único resultado final que refleje la incertidumbre causada por la imputación de los valores faltantes.

La imputación múltiple es una técnica poderosa para lidiar con los datos faltantes, ya que proporciona resultados más precisos y menos sesgados en comparación con otros métodos que simplemente eliminan las observaciones con valores faltantes. Sin embargo, la imputación múltiple es un proceso computacionalmente intensivo y requiere un conocimiento sólido de la teoría estadística para su implementación efectiva.

En este sentido, suponga entonces que existe un conjunto de \(n\) datos que relaciona dos variables \(X\), \(Y\), a través del siguiente modelo de regresión simple:

\[y_i = \beta x_i + \varepsilon_i\] Para todo individuo \(i = 1, \ldots, n.\), de tal manera que los errores tienen distribución normal con \(E(\varepsilon) = 0\) y \(Var(\varepsilon) = \sigma ^2\).

Ahora bien, se debe tener en cuenta, para utilizar la metodología, los siguientes atributos:

Sea \(Y_{Obs}\) los valores observados para un conjunto de individuos de tamaño \(n_1\).
Sea \(Y_{NoObs}\) los valores NO observados de la variable \(Y\) de tamaño \(n_0\), es decir, \(n_1 + n_0 = n\).
Suponga que sí fue posible observar los valores de la covariable \(X\) para todos los individuos en la muestra.

A modo de ejemplificar la teoría, realicemos un ejercicio de simulación de la siguiente manera. Simular un conjunto de \(n = 500\) datos con una pendiente \(\beta = 10\) y con una dispersión de \(\sigma = 2\). A su vez, el conjunto de datos tendrá \(n_0 = 200\) valores faltantes en la variable respuesta.

El algoritmo de simulación es el siguiente:

generar <- function(n = 500, n_0 = 200, 
                    beta = 10, sigma = 2){
  x <- runif(n)
  mu <- beta * x
  y <- mu + rnorm(n, mean = 0, sd = sigma)
  datos <- data.frame(x = x, y = y)
  faltantes <- sample(n, n_0)
  datos$faltantes <- "No"
  datos$faltantes[faltantes] <- "Si"
  datos$y.per <- y
  datos$y.per[faltantes] <- NA
  return(datos)
}

El código anterior realiza lo siguiente:

generar <- function(n = 500, n_0 = 200, beta = 10, sigma = 2){:Esta línea de código define la función “generar” con cuatro argumentos opcionales: n, n_0, beta y sigma.
x <- runif(n): Esta línea genera una secuencia de n números aleatorios uniformemente distribuidos en el intervalo [0,1] y los asigna a la variable “x”.
mu <- beta * x: Esta línea de código calcula la media condicional de “y” (la variable de respuesta) dada “x” utilizando el parámetro de pendiente “beta” y lo asigna a la variable “mu”.
y <- mu + rnorm(n, mean = 0, sd = sigma): Esta línea genera una variable de respuesta “y” a partir de la media condicional “mu” y agrega un error aleatorio generado a partir de una distribución normal con media cero y desviación estándar “sigma”.
datos <- data.frame(x = x, y = y): Esta línea combina las variables “x” e “y” en un data frame llamado “datos”.
faltantes <- sample(n, n_0): Esta línea genera una muestra aleatoria de “n_0” valores únicos desde un rango de 1 hasta “n” (la longitud de los datos) y los asigna a la variable “faltantes”.
datos\(faltantes <- "No"*, *datos\)faltantes[faltantes] <- “Si”: Esta línea crea una nueva columna en el data frame “datos” llamada “faltantes” y la inicializa con valores “No” para todas las observaciones. Luego, establece los valores de esta columna en “Si” para las filas seleccionadas por la muestra aleatoria anterior.
datos\(y.per <- y*, *datos\)y.per[faltantes] <- NA: Esta línea crea una nueva columna en el data frame “datos” llamada “y.per” y la inicializa con los mismos valores que “y”. Luego, establece los valores de esta columna en “NA” para las filas seleccionadas por la muestra aleatoria anterior, lo que simula valores faltantes en los datos de “y”.

Una vez creada la función anterior, se genera la población usando una semilla para poder reproducirla:

set.seed(1234)
datos <- generar()
head(datos,12)

x	y	faltantes	y.per
0.1137034	2.0108953	No	2.010895
0.6222994	8.3432419	No	8.343242
0.6092747	6.9971281	No	6.997128
0.6233794	7.5601916	Si	NA
0.8609154	6.3364067	No	6.336407
0.6403106	5.6621110	No	5.662111
0.0094958	3.0488967	No	3.048897
0.2325505	-0.1223024	Si	NA
0.6660838	7.1769744	Si	NA
0.5142511	5.9525170	No	5.952517
0.6935913	8.8875196	No	8.887520
0.5449748	4.7519949	No	4.751995

A continuación, se grafican la relación entre “x” y “y” y “x” y “y.per” notando que, la distribución entre las dos relaciones es muy similar.

library(patchwork)
p1 <- ggplot(data = datos, aes(x = x, y = y)) +
  geom_point() + geom_smooth(formula = y~x , method = "lm")

p2 <- ggplot(data = datos, aes(x = x, y = y.per)) +
  geom_point() + geom_smooth(formula = y~x , method = "lm")
  
p1 | p1

Ahora, dado el 40% de valores faltantes, es necesario imputar dichos valores. Para esto, utilizaremos la técnica de imputación múltiple propuesta por Rubin (1987). La idea consiste en generar \(M > 1\) conjuntos de valores para los datos faltantes. Al final, el valor imputado corresponderá al promedio de esos \(M\) valores.

Hay varias maneras de realizar la imputación, a continuación, se hace un listado:

Ingenua: Esta clase de imputación carece de aleatoriedad y por tanto, la varianza de \(\beta\) va a ser subestimada.
Bootstrap: Se seleccionan \(m\) muestras bootstrap, y para cada una se estiman los parámetros \(\beta\) y \(\sigma\) para generar \(\dot{y}_i\). Al final se promedian los \(m\) valores y se imputa el valor faltante.
Bayesiana: Se definen las distribuciones posteriores de \(\beta\) y \(\sigma\) para generar \(M\) valores de estos parámetros y por tanto \(M\) valores de \(\dot{y}_i\). Al final se promedian los \(M\) valores y se imputa el valor faltante.

Dado que el interés es la estimación de la pendiente de la regresión lineal simple \(\beta\), entonces la esperanza estimada al utilizar la metodología de imputación múltiple está dada por:

\[E(\hat{\beta} | Y_{obs}) = E(E(\hat{\beta} | Y_{obs}, Y_{mis}) | Y_{obs})\]

Esta expresión es estimada por el promedio de las \(M\) estimaciones puntuales de \(\hat{\beta}\) sobre las \(M\) imputaciones, dado por:

\[\bar{\hat{\beta}} = \frac{1}{M} \sum_{m = 1} ^ M \hat{\beta}_m\]

La varianza estimada al utilizar la metodología de imputación múltiple está dada por la siguiente expresión:

\[ V(\hat{\beta} | Y_{obs}) = E(V(\hat{\beta} | Y_{obs}, Y_{mis}) | Y_{obs}) + V(E(\hat{\beta} | Y_{obs}, Y_{mis}) | Y_{obs}) \] La primera parte de la anterior expresión se estima como el promedio de las varianzas muestrales de \(\hat{\beta}\) sobre las \(M\) imputaciones, dado por:

\[\bar{U} = \frac{1}{M} \sum_{m = 1} ^ M Var(\beta)\]

El segundo término se estima como la varianza muestral de las \(M\) estimaciones puntuales de \(\hat{\beta}\) sobre las \(M\) imputaciones, dada por:

\[B = \frac{1}{M-1} \sum_{m = 1} ^ M (\hat{\beta}_m - \bar{\hat{\beta}})\]

Es necesario tener en cuenta un factor de corrección (puesto que \(M\) es finito). Por tanto, la estimación del segundo término viene dada por la siguiente expresión: \[ (1 + \frac{1}{M}) B \]

Por tanto, la varianza estimada es igual a: \[\hat{V}(\hat{\beta} | Y_{obs}) = \bar{U} + (1 + \frac{1}{M}) B\] A modo de ejemplo, a continuación, se crea una función que haga la estimación de los valores faltantes usando bootstrap:

im.bootstrap <- function(datos, M = 15){
  library(dplyr)
  n <- nrow(datos)
  datos1 <- na.omit(datos)
  n1 <- nrow(datos1)
  n0 <- n - n1
  Ind <- is.na(datos$y.per)
  faltantes.boot <- NULL
  beta1 <- NULL
  sigma1 <- NULL
  
  for (m in 1:M){
    datos.m <- dplyr::sample_n(datos1, n1, replace = TRUE)
    model1 <- lm(y ~ 0 + x, data = datos.m)
    beta <- model1$coeff
    sigma <- sqrt(anova(model1)[["Mean Sq"]][2])
    faltantes.boot <- rnorm(n0, datos$x[Ind] * beta, sd = sigma)
    datos$y.per[Ind] <-  faltantes.boot
    model.input <- lm(y.per ~ 0 + x, data = datos)
    beta1[m] <- model.input$coeff
    sigma1[m] <- summary(model.input)$coeff[2]
  }
  beta.input <- mean(beta1)
  u.bar <- mean(sigma1 ^ 2)
  B <- var(beta1)
  beta.sd <- sqrt(u.bar + B + B/M)
  result <- list(new = datos, beta = beta.input, sd = beta.sd)
}

La función anterior realiza internamente lo siguiente:

En la primera sección del código, desde library(dplyr) hasta sigma1 <- NULL se calcula el número total de observaciones “n”, se eliminan las filas con valores faltantes y se calculan las longitudes de los datos “datos1” sin valores faltantes, “n1”, y los datos faltantes “n0”. También se crea una variable de índice “Ind” que indica qué observaciones tienen valores faltantes, y se inicializan varias variables que se usarán más adelante en la función.
Luego, en la sección del bucle for se realiza el proceso de imputación múltiple. Se itera “M” veces, donde en cada iteración se muestran “n1” filas de los datos sin valores faltantes “datos1” con reemplazo, y se ajusta un modelo de regresión lineal con “y” como respuesta y “x” como variable explicativa. Luego se calculan los coeficientes de la regresión y la desviación estándar de los residuos. A continuación, se generan valores de imputación aleatorios para las filas faltantes, utilizando la media condicional y la desviación estándar estimadas en la muestra. Estos valores de imputación se asignan a las filas faltantes en el data frame “datos”. Finalmente, se ajusta un modelo de regresión lineal con la variable de respuesta imputada “y.per” como respuesta y “x” como variable explicativa. Se guardan los coeficientes y las desviaciones estándar de los residuos para cada iteración en las variables “beta1” y “sigma1”, respectivamente.
En la última sección del código se calcula el coeficiente de regresión promedio para las iteraciones, “beta.input”, y se calcula la desviación estándar de “beta1” y por último, la estimación de la desviación estándar del beta.

Al aplicar la función sobre el conjunto de datos creado, se obtienen las siguientes salidas:

datos <- generar()
im.bootstrap(datos)$beta

## [1] 10.29741

im.bootstrap(datos)$sd

## [1] 0.2222246

head(im.bootstrap(datos)$new)

x	y	faltantes	y.per
0.2173207	0.2872457	Si	0.4423490
0.2953090	1.5860882	Si	1.8694104
0.9609104	11.2309918	Si	11.8596299
0.3119669	5.0508837	No	5.0508837
0.0520595	-0.7714404	Si	-0.5380997
0.1597426	3.9523667	No	3.9523667

Nótese que existe una buena dispersión en los valores imputados.

nuevos <- im.bootstrap(datos)$new
ggplot(data = nuevos, aes(x = x, y = y.per, color = faltantes)) + geom_point()

Por otro lado, se ejemplificará la técnica de imputación múltiple para los datos de la encuesta que se utiliza de ejemplo en este texto:

encuesta$Income_imp <- encuesta$Income_missin
encuesta$Employment_imp <- encuesta$Employment_missin
encuesta_obs <- filter(encuesta, !is.na(Income_missin))
encuesta_no_obs <- filter(encuesta, is.na(Income_missin))
n0 <- nrow(encuesta_no_obs)
n1 <- nrow(encuesta_obs)

Inicialmente, se extraen los datos a imputar y se calculan los tamaños de los datos observados y no observados.

M = 10 
set.seed(1234)
for (ii in 1:M) {
vp <- paste0("Income_vp_",ii)
vp2 <- paste0("Employment_vp_",ii)

encuesta_temp <- encuesta_obs %>% sample_n(size = n1, replace = TRUE)

mod <- lm(Income~ Zone + Sex + Expenditure, data = encuesta_temp)
mod.mult <- multinom(Employment~Zone + Sex +Expenditure,data = encuesta_temp)

encuesta_no_obs[[vp]] <- predict(mod, encuesta_no_obs)
encuesta_obs[[vp]] <- encuesta_obs$Income

encuesta_no_obs[[vp2]] <- predict(mod.mult, encuesta_no_obs,type = "class")
encuesta_obs[[vp2]] <- encuesta_obs$Employment
}

## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1147.241384
## final  value 1127.485787 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1177.374646
## final  value 1161.652154 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1104.676334
## final  value 1090.350507 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1157.667055
## final  value 1138.773189 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1145.090973
## final  value 1094.863541 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1139.350146
## final  value 1127.134266 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1169.810267
## final  value 1141.798357 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1165.750155
## final  value 1125.933563 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1117.370321
## final  value 1093.345603 
## converged
## # weights:  15 (8 variable)
## initial  value 1651.214270 
## iter  10 value 1143.115259
## final  value 1129.053716 
## converged

El código anterior realiza lo siguiente:

Se crea una cadena de caracteres “vp” y “vp2” mediante la función paste0() concatenando la cadena “Income_vp_” y “Employment_vp_” respectivamente con el número de la iteración actual “ii”.
Se crea una muestra aleatoria con reemplazo de tamaño n1 a partir de la base de datos “encuesta_obs” utilizando la función sample_n() del paquete dplyr.
Se ajusta un modelo de regresión lineal con la variable “Income” como respuesta y las variables “Zone”, “Sex” y “Expenditure” como covariables utilizando la función lm().
Se ajusta un modelo de regresión multinomial con la variable “Employment” como respuesta y las variables “Zone”, “Sex” y “Expenditure” como covariables utilizando la función multinom() del paquete nnet.
Se utiliza el modelo ajustado en el punto 3 para predecir los valores de la variable “Income” en la base de datos “encuesta_no_obs” mediante la función predict().
Se guarda en la base de datos “encuesta_obs” los valores verdaderos de la variable “Income”.
Se utiliza el modelo ajustado en el punto 4 para predecir los valores de la variable “Employment” en la base de datos “encuesta_no_obs” mediante la función predict().
Se guarda en la base de datos “encuesta_obs” los valores verdaderos de la variable “Employment”.
Se repite el proceso para las siguientes iteraciones.
Al finalizar el bucle se obtendrán 20 nuevas variables en las bases de datos “encuesta_no_obs” y “encuesta_obs” correspondientes a las predicciones de “Income” y “Employment” respectivamente, para cada una de las 10 iteraciones del bucle.

Una vez corrido el código anterior, se seleccionan las variables de ingresos y sus 10 valores plausibles como se muestra a continuación:

dplyr::select(encuesta_no_obs, Income, matches("Income_vp_"))[1:10,1:4]

Income	Income_vp_1	Income_vp_2	Income_vp_3
409.87	550.2195	566.0432	567.8296
409.87	561.1194	529.3457	541.8310
90.92	210.5682	225.7715	163.9913
90.92	221.4682	189.0739	137.9927
90.92	210.5682	225.7715	163.9913
135.33	222.6937	237.9191	178.4083
135.33	222.6937	237.9191	178.4083
1539.75	784.8579	801.1103	846.8098
336.00	507.8955	472.7913	439.7467
685.48	593.0111	558.0623	540.9473

A continuación, se grafica la distribución de los ingresos y los 10 valores plausibles observándose que, las distribuciones son muy similares:

encuesta <- bind_rows(encuesta_obs, encuesta_no_obs)

dat_plot10 <- tidyr::gather(
  encuesta %>% dplyr::select(Zone,Sex,matches("Income_vp_")),
  key = "Caso", value = "Income2", -Zone,-Sex)

p1 <- ggplot(dat_plot10, aes(x = Income2, col = Caso)) + 
  geom_density(alpha = 0.2) + theme_bw() +
   theme(legend.position = "bottom") + geom_density(data = encuesta ,aes(x = Income), col =  "black", size = 1.2) 

p1

En el siguiente gráfico se presentan las frecuencias de los valores plausibles para la variable empleado. Se observa también que son muy aproximadas a la variable completa:

dat_plot11 <- tidyr::gather(
  encuesta %>% 
  dplyr::select(Zone,Sex, Employment,matches("Employment_vp_")),
  key = "Caso", value = "Employment2", -Zone,-Sex) %>%
  group_by(Caso,Employment2) %>% tally() %>% 
  group_by(Caso) %>% mutate(prop = n/sum(n))

p1 <- ggplot(dat_plot11, 
        aes(x = Employment2, y = prop,
            fill = Caso, color="red")) + 
       geom_bar(stat="identity",
          position = position_dodge(width = 0.5))  +
   theme_bw() +
   theme(legend.position = "bottom") +
  scale_fill_manual(values = c("Employment" = "black"))
p1

Con los valores plausibles enocntrados anteriormente, se procede a definir el diseño muestral utilizado en este ejemplo y así poder hacer la estimación de los parámetros. A continuación, se define el diseño muestral:

library(srvyr)

diseno <- encuesta %>%
  as_survey_design(
    strata = Stratum,
    ids = PSU,
    weights = wk,
    nest = T
  )

Con el diseño anterior, se estiman los ingresos medios para cada valor plausible junto ocn su varianza, como se muestra a continuación:

estimacion_vp <-  diseno %>% 
 summarise(
   vp1 = survey_mean(Income_vp_1, vartype = c("var")),
   vp2 = survey_mean(Income_vp_2, vartype = c("var")),
   vp3 = survey_mean(Income_vp_3, vartype = c("var")),
   vp4 = survey_mean(Income_vp_4, vartype = c("var")),
   vp5 = survey_mean(Income_vp_5, vartype = c("var")),
   vp6 = survey_mean(Income_vp_6, vartype = c("var")),
   vp7 = survey_mean(Income_vp_7, vartype = c("var")),
   vp8 = survey_mean(Income_vp_8, vartype = c("var")),
   vp9 = survey_mean(Income_vp_9, vartype = c("var")),
   vp10 =survey_mean(Income_vp_10, vartype = c("var")))
estimacion_vp

vp1	vp1_var	vp2	vp2_var	vp3	vp3_var	vp4	vp4_var	vp5	vp5_var	vp6	vp6_var	vp7	vp7_var	vp8	vp8_var	vp9	vp9_var	vp10	vp10_var
619.3525	845.706	617.6058	844.5521	617.4322	867.5133	617.748	856.7732	619.9729	857.805	617.1098	852.6605	618.3378	860.7981	619.2779	867.9485	616.9401	850.131	615.7853	870.0434

A continuación se presentan los datos anteriores discriminado por promedio y varianza:

require(tidyr)
(estimacion_vp %<>% tidyr::gather() %>%  separate(key, c("vp", "estimacion")) %>% 
mutate(estimacion = ifelse(is.na(estimacion), "promedio","var" )) %>% spread(estimacion,value) %>% mutate(vp = 1:10))

vp	promedio	var
1	619.3525	845.7060
2	615.7853	870.0434
3	617.6058	844.5521
4	617.4322	867.5133
5	617.7480	856.7732
6	619.9729	857.8050
7	617.1098	852.6605
8	618.3378	860.7981
9	619.2779	867.9485
10	616.9401	850.1310

Por último, para obtener la estimación de la media y su varianza utilizando la imputación múltiple, se realizan los siguientes cálculos que se derivan de las expresiones matemáticas antes mostradas:

Media_vp = mean(estimacion_vp$promedio)
(Ubar = mean(estimacion_vp$var))

## [1] 857.3931

(B = var(estimacion_vp$promedio))

## [1] 1.645758

var_vp = Ubar + (1 + 1/M) 
(resultado <- data.frame(Media_vp, 
                        Media_vp_se = sqrt(var_vp)))

Media_vp	Media_vp_se
617.9562	29.30005

estimacion_var_vp <-  diseno %>% 
  summarise_at(vars(matches("Income_vp")), 
               survey_var,  vartype = "var" )

Por otro lado, otro parámetro de interés es la varianza de los ingresos. Este parámetro permite medir la variabilidad de los ingresos de los ciudadanos de la base de datos de ejemplo. La forma de estimarla es la misma que para el promedio de los ingresos y se utilizarán los mismo códigos mostrados anteriormente, cambiando el parámetro a estimar:

(estimacion_var_vp %<>% tidyr::gather() %>%
   separate(key, c("A", "B","vp", "estimacion")) %>% 
mutate(estimacion = ifelse(is.na(estimacion), "promedio","var" ),
       A = NULL, B = NULL, vp = as.numeric(vp)) %>% 
  spread(estimacion,value))

vp	promedio	var
1	262689.8	3074674460
2	263092.0	3079895581
3	274370.1	3237991671
4	269127.4	3164842634
5	270450.0	3165285304
6	264992.6	3106582542
7	270916.0	3175865575
8	276069.5	3251836316
9	265060.8	3111426636
10	275462.0	3258214990

Por último, se utilizan las ecuaciones mostradas anteriormente:

Media_var_vp = mean(estimacion_var_vp$promedio)
Ubar = mean(estimacion_var_vp$var)
B = var(estimacion_var_vp$promedio)
var_var_vp = Ubar + (1 + 1/M)*B 
resultado$var_vp <- Media_var_vp
resultado$var_vp_se <- sqrt(var_var_vp)
cbind(Media_var_vp, var_var_vp)

Media_var_vp	var_var_vp
269223	3191037330

Otro parámetro de interés a estimar es la proporción. A continuación, se realizará la estimación de la proporción utilizando valores plausibles. Para ello se estimará la variable empleado, como se muestra a continuación:

estimacion_prop_vp <- lapply(paste0("Employment_vp_",1:10),
       function(vp){diseno %>% group_by_at(vars(Employment = vp)) %>% 
  summarise(prop = survey_mean(vartype = "var"),.groups = "drop") %>%
         mutate(vp = vp)}) %>% bind_rows()

Se presenta la estimación de la proporción para cada uno de los 10 valores plausibles en cada categoría de la variable:

(estimacion_prop_vp %<>% separate(vp, c("A", "B","vp")) %>% 
mutate(A = NULL, B = NULL, vp = as.numeric(vp)) %>%
  dplyr::select(vp,Employment:prop_var)) %>% slice(1:12L)

vp	Employment	prop	prop_var
1	Unemployed	0.0355675	0.0000375
1	Inactive	0.3753767	0.0002041
1	Employed	0.5890558	0.0001669
2	Unemployed	0.0355675	0.0000375
2	Inactive	0.3771091	0.0002035
2	Employed	0.5873234	0.0001661
3	Unemployed	0.0355675	0.0000375
3	Inactive	0.3867965	0.0002056
3	Employed	0.5776360	0.0001790
4	Unemployed	0.0355675	0.0000375
4	Inactive	0.3549487	0.0001915
4	Employed	0.6094838	0.0001524

Por último, utilizando las ecuaciones de Rubin se obtiene la varianza estimada:

resultado = estimacion_prop_vp %>% 
  group_by(Employment) %>% 
  summarise(prop_pv = mean(prop),
            Ubar = mean(prop_var),
            B = var(prop)) %>% 
  mutate(prop_pv_var = Ubar + (1 + 1/M)*B) |> 
  dplyr::select(Employment, prop_pv, prop_pv_var)
resultado

Employment	prop_pv	prop_pv_var
Unemployed	0.0355675	0.0000375
Inactive	0.3868270	0.0004665
Employed	0.5776055	0.0004328