9.3 Preparando insumos para STAN | Estimación de áreas pequeñas con el modelo multinomial

9.3 Preparando insumos para `STAN`

Lectura y adecuación de covariables

satelital_media <- readRDS('01 Modelo de area/PER/2017/Data/satelital_media.rds')
statelevel_censo <- readRDS('01 Modelo de area/PER/2017/Data/statelevel_predictors_df_dam2.rds') %>% 
  mutate(dam2 = case_when(dam2 == "120606"~"120699",
                          TRUE ~ dam2))

statelevel_predictors_df <- left_join(satelital_media,statelevel_censo) 
## Estandarizando las variables para controlar el efecto de la escala. 
statelevel_predictors_df %<>%
  mutate_at(vars("luces_nocturnas", 
                 "cubrimiento_cultivo",
                 "cubrimiento_urbano",
                 "modificacion_humana",
                 "accesibilidad_hospitales",
                 "accesibilidad_hosp_caminado"),
            function(x)as.numeric(scale(x)))

Seleccionar las variables del modelo y crear matriz de covariables.

names_cov <- c(
  "dam2",
  "luces_nocturnas",
  "cubrimiento_cultivo",
  "cubrimiento_urbano",
  "modificacion_humana",
  "tasa_desocupacion",
  "discapacidad1",
  "rezago_escolar"
) 
X_pred <-
  anti_join(statelevel_predictors_df %>% select(all_of(names_cov)),
            indicador_dam1 %>% select(dam2))

En el bloque de código se identifican que dominios serán los predichos.

X_pred %>% select(dam2) %>% 
  saveRDS(file = "01 Modelo de area/PER/2017/Data/dam_pred.rds")

Creando la matriz de covariables para los dominios no observados (X_pred) y los observados (X_obs)

## Obteniendo la matrix 
X_pred %<>%
  data.frame() %>%
  select(-dam2)  %>%  as.matrix()

## Identificando los dominios para realizar estimación del modelo

X_obs <- inner_join(indicador_dam1 %>% select(dam2, id_orden),
                    statelevel_predictors_df %>% select(all_of(names_cov))) %>%
  arrange(id_orden) %>%
  data.frame() %>%
  select(-dam2, -id_orden)  %>%  as.matrix()

Calculando el n_efectivo y el \(\tilde{y}\)

D <- nrow(indicador_dam1)
P <- 3 # Ocupado, desocupado, inactivo.
Y_tilde <- matrix(NA, D, P)
n_tilde <- matrix(NA, D, P)


# n efectivos ocupado
n_tilde[,1] <- (indicador_dam1$Ocupado*(1 - indicador_dam1$Ocupado))/indicador_dam1$Ocupado_var
Y_tilde[,1] <- n_tilde[,1]* indicador_dam1$Ocupado


# n efectivos desocupado
n_tilde[,2] <- (indicador_dam1$Desocupado*(1 - indicador_dam1$Desocupado))/indicador_dam1$Desocupado_var
Y_tilde[,2] <- n_tilde[,2]* indicador_dam1$Desocupado

# n efectivos Inactivo
n_tilde[,3] <- (indicador_dam1$Inactivo*(1 - indicador_dam1$Inactivo))/indicador_dam1$Inactivo_var
Y_tilde[,3] <- n_tilde[,3]* indicador_dam1$Inactivo

Compilando el modelo

X1_obs <- cbind(matrix(1,nrow = D,ncol = 1),X_obs)
K = ncol(X1_obs)
D1 <- nrow(X_pred)
X1_pred <- cbind(matrix(1,nrow = D1,ncol = 1),X_pred)

sample_data <- list(D = D,
                    P = P,
                    K = K,
                    hat_y = Y_tilde,
                    X_obs = X1_obs,
                    X_pred = X1_pred,
                    D1 = D1)


library(rstan)
fit_mcmc2 <- stan(
  file = "01 Modelo de area/0funciones/01 Multinomial_simple_pred.stan",  # Stan program
  data = sample_data,    # named list of data
  verbose = TRUE,
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 4,              # number of cores (could use one per chain)
)

saveRDS(fit_mcmc2,
        "01 Modelo de area/PER/2017/Data/fit_rtanmultinomial_con_covariable_satelite.Rds")