Según la imputación en tres pasos ver presentación imputación.ppt Este es el ejercicio basado en los datos ENIGH 2012.

Cargar los datos desde la pagina Web de INEGI

url.zip<-"https://www.inegi.org.mx/contenidos/programas/enigh/tradicional/2012/microdatos/Tra_Concentrado_2012_concil_2010_csv.zip"
## para el caso de archivos más recientes por ejemplo, 2012 usar la siguiente ruta ## url.zip<-"https://www.inegi.org.mx/contenidos/programas/enigh/nc/2022/microdatos/enigh2022_ns_concentradohogar_csv.zip"

nombre_zip<-"Tra_Concentrado_2012_concil_2010_csv.zip" 

# para el caso 2022 el nombre del archivo es concentradohogar.csv

# crea un directorio "enigh_2012_datos" en la ruta actual de trabajo para que se #pueda guardar los datos. 
dir_unzip <- "enigh_2012_datos"
if (!dir.exists(dir_unzip)) {
  dir.create(dir_unzip)
}
 #descarga y descomprime  el archivo zip 
download.file(url.zip, destfile = nombre_zip, mode = "wb")
unzip(nombre_zip, exdir = dir_unzip)
#guarda la ruta donde se encuentra el archivo para despues leerlo
ruta_csv <- file.path(dir_unzip, "tra_concentrado_2012_concil_2010.csv") 
# caso 2022 concentradohogar.csv


#los datos con los que se trabajaran
data12<- read.csv(ruta_csv, encoding = "Latin-1", na.strings = c("", " "))
names(data12) # todas las variables que contiene el archivo 
file.remove(nombre_zip)
#los datos con los que se trabajaran
dir_unzip <- "enigh_2012_datos"
ruta_csv <- file.path(dir_unzip, "tra_concentrado_2012_concil_2010.csv")
data12<- read.csv(ruta_csv, encoding = "Latin-1", na.strings = c("", " "))
names(data12) # todas las variables que contiene el archivo 
##   [1] "folioviv"   "foliohog"   "ubica_geo"  "ageb"       "tam_loc"   
##   [6] "est_socio"  "est_dis"    "upm"        "factor_hog" "clase_hog" 
##  [11] "sexo_jefe"  "edad_jefe"  "educa_jefe" "tot_integ"  "hombres"   
##  [16] "mujeres"    "mayores"    "menores"    "p12_64"     "p65mas"    
##  [21] "ocupados"   "percep_ing" "perc_ocupa" "ing_total"  "ing_cor"   
##  [26] "ing_mon"    "trabajo"    "sueldos"    "horas_extr" "comisiones"
##  [31] "otra_rem"   "negocio"    "noagrop"    "industria"  "comercio"  
##  [36] "servicios"  "agrope"     "agricolas"  "pecuarios"  "reproducc" 
##  [41] "pesca"      "otros_trab" "rentas"     "utilidad"   "arrenda"   
##  [46] "transfer"   "jubilacion" "becas"      "donativos"  "remesas"   
##  [51] "bene_gob"   "otros_ing"  "gasto_nom"  "autoconsum" "remu_espec"
##  [56] "transf_esp" "transf_hog" "trans_inst" "estim_alqu" "percep_tot"
##  [61] "percep_mon" "retiro_inv" "prestamos"  "otras_perc" "erogac_nom"
##  [66] "gasto_tot"  "gasto_cor"  "gasto_mon"  "alimentos"  "ali_dentro"
##  [71] "cereales"   "carnes"     "pescado"    "leche"      "huevo"     
##  [76] "aceites"    "tuberculo"  "verduras"   "frutas"     "azucar"    
##  [81] "cafe"       "especias"   "otros_alim" "bebidas"    "ali_fuera" 
##  [86] "tabaco"     "vesti_calz" "vestido"    "calzado"    "vivienda"  
##  [91] "alquiler"   "pred_cons"  "agua"       "energia"    "limpieza"  
##  [96] "cuidados"   "utensilios" "enseres"    "salud"      "atenc_ambu"
## [101] "hospital"   "medicinas"  "transporte" "publico"    "foraneo"   
## [106] "adqui_vehi" "mantenim"   "refaccion"  "combus"     "comunica"  
## [111] "educa_espa" "educacion"  "esparci"    "paq_turist" "personales"
## [116] "cuida_pers" "acces_pers" "otros_gas"  "transf_gas" "erogac_tot"
## [121] "erogac_mon" "cuota_viv"  "mater_serv" "material"   "servicio"  
## [126] "deposito"   "prest_terc" "pago_tarje" "deudas"     "balance"   
## [131] "otras_erog" "smg"
#file.remove(nombre_zip)

El siguiente código es para obtener el vector de valores de ingresos y factores de expansión

 muestra<-data12$ing_cor
 fac<-data12$factor # aveces el factor de expansión tiene el nombre "factor_hog"
 total_ingreso <-sum(muestra*fac) # 
 tot_hog<-sum(fac) 
 print(paste("total del ingreso de la muestra: ",total_ingreso,sep=""))
## [1] "total del ingreso de la muestra: 1203202597623.71"
 print(paste("total de hogares en la población: ",tot_hog,sep=""))
## [1] "total de hogares en la población: 31559379"

Se quitan los valores menores a 2000 y se redistribuye la población

w<-c()
w<-which(muestra<2000)
if(length(w)==0)
{muestra<-muestra} else 
{muestra<-muestra[-w]
 fac<-fac[-w]
}
fquitados<-tot_hog-sum(fac)  ## re-austar los valores de expansión 
coci<-fquitados/length(fac)
nfac<-(fac+coci)
fac<-nfac  ## nuevos factores ajustados 
w<-order(muestra)
muestra<-muestra[w]
fac<-fac[w]
head(fac) # los nuevos factores de expansión 
## [1] 1967.884 6202.884 2484.884 2854.884 1730.884 3189.884

Algunos datos obtenidos del ingreso del año 2012

El valor estimado del vector de parámetros de la función GG sin restricciones ( obtenido antes)

library(alabama)
## Cargando paquete requerido: numDeriv
library(gamlss)
## Warning: package 'gamlss' was built under R version 4.4.3
## Cargando paquete requerido: splines
## Cargando paquete requerido: gamlss.data
## 
## Adjuntando el paquete: 'gamlss.data'
## The following object is masked from 'package:datasets':
## 
##     sleep
## Cargando paquete requerido: gamlss.dist
## Warning: package 'gamlss.dist' was built under R version 4.4.3
## Cargando paquete requerido: nlme
## Cargando paquete requerido: parallel
##  **********   GAMLSS Version 5.4-22  **********
## For more on GAMLSS look at https://www.gamlss.com/
## Type gamlssNews() to see new features/changes/bug fixes.
xGGsin<-c(25084.35679,0.805540247,-0.199996903)
xGGsin
## [1] 25084.3567900     0.8055402    -0.1999969

El valor del vector estimado del parámetro de la función GG con restricciones ( obtenido antes)

xGGcon<-c(29580.06446, 1.068264447, -0.47753398)
xGGcon
## [1] 29580.064460     1.068264    -0.477534

Algunos datos importantes del 2012

a=min(muestra)
b<-max(muestra)
tot_hog<-sum(fac) # total de hogares 
p=0.999999 #Para el caso del ajuste a nivel nacional para algún año de interés nos interesa: el promedio de ingresos declarado al SAT por el millonesimo percentil(.999999) que corresponde los 32 más ricos, y el ingreso del treintaidosavo más rico.
mean_Sat_anual<- 1624174001.00 # promedio del ingreso de los 32 más ricos anual                   
mean_sat<-mean_Sat_anual/4 # trimestralizado 2012
qi<-177965502 # ingreso del 32 del top-32 más ricos ya trimestralizado
prom<- mean_sat # promedio del ingreeso de los 32 más ricos según SAT 2012 trimestralizado
prom<-prom*0.000001
idan<-  11635413   # valor de IDAN del 2012 
tot_cn<-idan *1000000 # se multiplica por 1 millón  
ci<-c(tot_cn/4/tot_hog)  # promedio de ingreso por hogar según CN

funciones truncadas

GGtrunc_max<-function(x, m=muestra,w=fac,a=min(muestra),b=max(muestra))
{
  fi <- c()
  mu<-x[1]
  sigma<-x[2]
  nu<-x[3]
  Fa<-pGG(a, mu, sigma, nu) ## devuelve probabilidad F(x)
  Fb<-pGG(b, mu, sigma, nu)
  fi <- (w*(dGG(m, mu, sigma, nu,log=TRUE)-log(Fb-Fa)))
 sumaf <- -1*sum( fi )
 #sumaf <- sum( fi )
  return( sumaf )
}

hin<-function(x)
{
  h <- rep(NA, 2)
  h[1]<-x[1]
  h[2]<-x[2]
  h
}

valor inicial para obtener

theta<-c(25084.35679,0.805540247,-0.199996903) # valor inicial para empezar la optimizacion 2012
ansGGsin <- constrOptim.nl(par=theta, fn=GGtrunc_max,hin=hin) 
## Min(hin):  0.8055402 
## par:  25084.36 0.8055402 -0.1999969 
## fval:  359727156 
## Min(hin):  0.8007333 
## par:  24802 0.8007333 -0.2434864 
## fval:  359723467 
## Min(hin):  0.8007593 
## par:  24802 0.8007593 -0.2434856 
## fval:  359723467 
## Min(hin):  0.8007593 
## par:  24802 0.8007593 -0.2434856 
## fval:  359723467
parsin<-ansGGsin$par  # parámetro sin restricciones para función truncada

Las probabilidades truncadas de la muestra

parsin # parámetro de la vesión truncada
## [1] 24802.0001781     0.8007593    -0.2434856
# las probabilidades de la muestra
pp<-pGG(muestra, mu=parsin[1],sigma=parsin[2], nu=parsin[3],lower.tail = TRUE, log.p = FALSE)
head(pp)
## [1] 0.0001848681 0.0002047245 0.0002221342 0.0002316695 0.0006235241
## [6] 0.0006291202

ingreso imputado

xGGcon<-c(29580.06446, 1.068264447, -0.47753398)  # con restricciones 2012
# anteriormente para el 2012
qq<-qGG(pp, mu=xGGcon[1],sigma=xGGcon[2], nu=xGGcon[3],lower.tail = TRUE, log.p = FALSE)
head(qq) #" valores imputados"
## [1] 1740.019 1770.855 1796.147 1809.395 2172.274 2176.060

Obtención de algunos valores importantes del truncamiento

#limites
par<-xGGcon
aprima<-min(qq) # valor mínimo ingreso  de la función sin restricciones truncada 
bprima<-max(qq) # valor max de ingreso  de la función sin restricciones truncada 
fa<-min(pp)   # probabilidad función sin restricciones truncada para aprima
fb<-max(pp)   # probabilidad función sin restricciones truncada para bprima
# integral x*f(x)
n<-length(qq)

promedio_ingreso_aprima<-integrate(function(x) x*dGG(x=x, mu=par[1], sigma=par[2],nu=par[3]), 
          lower=0, upper=aprima)$v  # note los limites 

promedio_ingreso_aprima_bprima<-integrate(function(x) x*dGG(x=x, mu=par[1], sigma=par[2],nu=par[3]), 
          lower=aprima, upper= bprima)$v  # note los limites

promedio_ingreso_bprima<-integrate(function(x) x*dGG(x=x, mu=par[1], sigma=par[2],nu=par[3]), 
          lower=bprima, upper= 1000000000000)$v   # note los limites
ingreso_muestra<-sum(muestra*fac)
ingreso_muestra
## [1] 1.204045e+12
hogares_0_aprima<-fa*tot_hog
hogares_0_aprima
## [1] 5834.323
hogares_aprima_bprima<-(fb-fa)*tot_hog
hogares_aprima_bprima
## [1] 31550479
hogares_bprima_<-(1-fb)*tot_hog
hogares_bprima_
## [1] 3065.88
total_ingreso_aprima<-(promedio_ingreso_aprima/fa)*hogares_0_aprima
total_ingreso_aprima_bprima<-(promedio_ingreso_aprima_bprima/(fb-fa))*hogares_aprima_bprima
total_ingreso_bprima<-(promedio_ingreso_bprima/(1-fb))*hogares_bprima_