6.3 Modelo para la varianza

El código ajusta un modelo de regresión lineal múltiple (utilizando la función lm()), donde ln_sigma2 es la variable respuesta y las variables predictoras son pobreza, nd, y varias transformaciones de éstas. El objetivo de este modelo es estimar la función generalizada de varianza (FGV) para los dominios observados.

library(gtsummary)
FGV1 <- lm(ln_sigma2 ~ pobreza + 
             I(sqrt(nd)) + 
             I(sqrt(pobreza)) +
             I(sqrt(nd^pobreza)) 
             ,  
     data = baseFGV)
tbl_regression(FGV1) %>% 
  add_glance_table(include = c(r.squared, adj.r.squared))
Characteristic Beta 95% CI1 p-value
pobreza -40 -64, -15 0.002
I(sqrt(nd)) -0.06 -0.08, -0.03 <0.001
I(sqrt(pobreza)) 31 15, 47 <0.001
I(sqrt(nd^pobreza)) 1.4 0.14, 2.7 0.030
R² 0.298
Adjusted R² 0.261
1 CI = Confidence Interval

Después de tener la estimación del modelo se debe obtener el valor de la constante \(\Delta\) para lo cual se usa el siguiente código.

delta.hat = sum(baseFGV$vardir) / 
  sum(exp(fitted.values(FGV1)))

De donde se obtiene que \(\Delta = 1.6825614\). Final es posible obtener la varianza suavizada ejecutando el siguiente comando.

hat.sigma <- 
  data.frame(dam2 = baseFGV$dam2,
             hat_var = delta.hat * exp(fitted.values(FGV1)))

baseFGV <- left_join(baseFGV, hat.sigma)
tba(head(baseFGV, 10))
dam2 pobreza nd vardir ln_sigma2 hat_var
00101 0.2225 6796 0.0004 -7.9505 0.0004
00201 0.1822 531 0.0004 -7.7493 0.0039
00206 0.3366 230 0.0031 -5.7661 0.0041
00301 0.4266 666 0.0043 -5.4432 0.0050
00302 0.4461 261 0.0014 -6.5697 0.0029
00303 0.5587 566 0.0142 -4.2552 0.0075
00304 0.5406 412 0.0116 -4.4550 0.0042
00401 0.3359 1219 0.0010 -6.9320 0.0042
00402 0.1496 172 0.0007 -7.2136 0.0047
00403 0.4644 309 0.0015 -6.4935 0.0031

Validación del modelo para la FGV

par(mfrow = c(2, 2))
plot(FGV1)

Comparación entre la varianza estimada versus la pronosticada por la FGV

ggplot(baseFGV , 
       aes(y = vardir, x = hat_var)) + 
  geom_point() +
  geom_smooth(method = "loess") + 
    labs(x = "FGV", y = "VarDirEst") +
  ylab("Varianza del Estimador Directo")

Predicción de la varianza suavizada

base_sae <- base_sae %>%  left_join(hat.sigma, by = "dam2")

Ahora, realizamos un gráfico de linea para ver la volatilidad es la estimaciones de las varianzas.

ggplot(base_sae %>%
         arrange(nd), aes(x = 1:nrow(base_sae))) +
  geom_line(aes(y = vardir, color = "VarDirEst")) +
  geom_line(aes(y = hat_var, color = "FGV")) +
  labs(y = "Varianzas", x = "Tamaño muestral", color = " ") +
  scale_x_continuous(breaks = seq(1, nrow(base_sae), by = 10),
                     labels = base_sae$nd[order(base_sae$nd)][seq(1, nrow(base_sae), by = 10)]) +
  scale_color_manual(values = c("FGV" = "Blue", "VarDirEst" = "Red"))

El siguiente código utiliza la función mutate() del paquete dplyr para crear nuevas variables de la base de datos base_sae y luego guarda el resultado en un archivo RDS llamado base_FH_2018.rds.

En concreto, el código realiza las siguientes operaciones:

  • La variable deff_dam2 se ajusta a 1 cuando es NaN.

  • La variable deff_FGV se calcula a partir de otras dos variables hat_var y vardir. Si vardir es 0, entonces deff_FGV se ajusta a 1. En caso contrario, se divide hat_var por vardir / deff_dam2 para obtener deff_FGV.

  • La variable deff_FGV se regulariza utilizando el criterio MDS: si deff_FGV es menor que 1, se ajusta a 1.

  • Finalmente, se calcula la variable n_eff_FGV dividiendo nd (el tamaño de la muestra) por deff_FGV.

base_FH <- base_sae %>%
  mutate(
    deff_dam2 = ifelse(is.nan(deff_dam2), 1,
                         deff_dam2),
    deff_FGV = ifelse(
      vardir == 0 ,
      1,
      hat_var / (vardir / deff_dam2)
    ),
    # Criterio MDS para regularizar el DeffFGV
    deff_FGV = ifelse(deff_FGV < 1, 1, deff_FGV),
    n_eff_FGV = nd / deff_FGV
  )

saveRDS(object = base_FH, "Recursos/Día2/Sesion2/Data/base_FH_2018.rds")