upax_library: Librería de análisis

#' función de muestreo
#'
#' Función para calcular muestras complejas
#'
#' @param datos dataframe - marco muestral
#' @param var_estratos char vector - nombres de las variables con las que se construyen los estratos
#' @param var_tamano character - nombre de la variable con los tamanos
#' @param casos numeric - número de casos en la muestra
#' @param casos_x_unidad numeric - número de casos por unidad seleccionada
#' @param tolerancia_tamano numeric - multiplicador de tolerancia. Default = 1.5
#' @param unidad character - nombre de la variable que contiene el ID de las unidades seleccionadas
#' @param casos_maximos_seccion numeric - número máximo de casos por unidad
#' @param repeticiones_muestreo_estratos numeric - número de muestras generadas, siempre se seleccionará la muestra con mayor verosimilitud con base en ese número. Default = 1
#' @param semilla numeric - semilla aleatoria utilizada para mantener la replicabilidad. Default = 1
#' @param minimo_duplicar numeric - número mayor a 0 y menor a 1 que decir del punto de corte para n-plicar. Default = .8
#' @param repeticiones_muestreo_unidades numeric - número de asignaciones probabilísticas generadas, la asignación con mayor verosimilitud se selecciona. Default = 1
#' @param forzar_proporcionalidad numeric - número entre 0 y 1 que define el porcentaje a partir del cuál se forzará una asignación proporcional preferencial cuando se hace la asignación estratificada
#' @param porcentaje_aceptacion_estrato_unico numeric - número entre 0 y 1 que determina el tamaño mínimo que tendrá algo para ser considerado como estrato. Default = .8
#' @param fusiona_estratos boolean - valor que indica si el último estrato creado se fusionará con el penúltimo. Default = T
#' 
#' @return dataframe - el marco muestral original con 2 nuevas columnas correspondientes a la selección de la muestra (0,1) y a las probabilidades de inclusión
#' 
#' @details La estratificación para el caso electoral considera generalmente los municipios y los distritos federales, en el caso no electoral tomamos una combinación de municipios, colonias y códigos postales. 
#' 
#' Cuando una unidad tiene un tamaño igual o mayor a N / n se obtiene una probabilidad de inclusión 1 (Salvo ajustes cuando existen múltiples instancias de estas unidades).
#' 
#' La verosimilitud se refiere a la probabilidad de haber visto una ejecución específica de un evento aleatorio. Al quedarnos con muestras más verosímiles disminuímos la probabilidad de quedarnos con una muestra desafortunada
#' 
#' El muestreo se realiza con una estratificación seguida de un muestreo sistemático de probabilidades proporcionales al tamaño.


a_muestra <- function (datos, var_estratos, var_tamano, casos, casos_x_unidad, 
          tolerancia_tamano = 1.5, unidad = "seccion", casos_maximos_seccion, 
          repeticiones_muestreo_estratos = 1, repeticiones_muestreo_unidades = 1, 
          semilla = 42, minimo_duplicar = 0.8, forzar_proporcionalidad = 0, 
          porcentaje_aceptacion_estrato_unico = 0.8, fusiona_estratos = T) 
{
  require(sampling)
  require(tidyr)
  round_preserve_sum <- function(x, digits = 0) {
    up <- 10^digits
    x <- x * up
    y <- floor(x)
    indices <- tail(order(x - y), round(sum(x)) - sum(y))
    y[indices] <- y[indices] + 1
    y/up
  }
  funcion_trim_estratos <- function(tabla, forzar_proporcionalidad = 0) {
    resultado <- tabla
    if (forzar_proporcionalidad > 0) {
      casos_totales <- sum(tabla$casos)
      resultado <- resultado %>% mutate(porcentaje = if_else(porcentaje < 
                                                               forzar_proporcionalidad, 0, porcentaje), porcentaje = porcentaje/sum(porcentaje), 
                                        casos = porcentaje * casos_totales) %>% filter()
    }
    return(resultado)
  }
  funcion_asignacion_proporcional <- function(datos, estratos, 
                                              tamano, casos) {
    resultado <- eval(parse(text = paste0("datos %>% \n      group_by(", 
                                          paste0(estratos, collapse = ","), ") %>% \n      summarise(tamano_estrato = sum(", 
                                          tamano, ")) %>%\n      ungroup %>%\n      mutate(porcentaje = tamano_estrato/sum(tamano_estrato),\n      casos = porcentaje * ", 
                                          casos, ")")))
    return(resultado)
  }
  funcion_creacion_estratos <- function(asignacion_proporcional, 
                                        casos_seccion, estratos, porcentaje_aceptacion_estrato_unico, 
                                        fusiona_estratos = T) {
    resultado <- asignacion_proporcional %>% arrange(-casos)
    resultado_normal <- resultado %>% filter(casos >= casos_seccion * 
                                               porcentaje_aceptacion_estrato_unico)
    resultado_normal$temporal = 0
    resultado_normal$estrato <- eval(parse(text = paste0("resultado_normal$", 
                                                         estratos)))
    resultado_resto <- resultado %>% filter(casos < casos_seccion * 
                                              porcentaje_aceptacion_estrato_unico) %>% mutate(temporal = casos) %>% 
      mutate(estrato = "estrato1")
    if (nrow(resultado_resto) > 0) {
      estrato = 1

      for(i in 1:nrow(resultado_resto)){
        resultado_resto$estrato[i] <- paste0('estrato',estrato)
        if(sum(resultado_resto$temporal[1:i]) > casos_seccion){

          resultado_resto$temporal[1:i] <- 0
          if (i != nrow(resultado_resto)) 
            estrato <- estrato + 1
        }
      }
      if (estrato > 1 & fusiona_estratos == T) {
        resultado_resto1 <- eval(parse(text = paste0("resultado_resto %>% \n        filter(estrato != 'estrato", 
                                                     estrato, "')")))
        resultado_resto2 <- eval(parse(text = paste0("resultado_resto %>% \n        filter(estrato == 'estrato", 
                                                     estrato, "')")))
        suma <- resultado_resto2 %>% summarise(suma = sum(casos)) %>% 
          unlist
        if (suma < casos_seccion * porcentaje_aceptacion_estrato_unico) {
          resultado_resto2$estrato <- paste0("estrato", 
                                             estrato - 1)
        }
      }
      else {
        resultado_resto1 <- resultado_resto
        resultado_resto2 <- resultado_resto[0, ]
      }
      resultado <- rbind(resultado_normal, resultado_resto1, 
                         resultado_resto2)
    }
    else {
      resultado <- resultado_normal
    }
    return(resultado)
  }
  funcion_seleccion_intraestrato <- function(asignacion_proporcional) {
    f_seleccion <- function(tabla) {
      if (nrow(tabla) == 1) {
        tabla$inclusion = 1
        tabla$probabilidades_estrato = 1
      }
      else {
        probabilidades <- inclusionprobabilities(tabla$tamano_estrato, 
                                                 1)
        tabla$inclusion <- 0
        tabla$inclusion <- UPmaxentropy(probabilidades)
        tabla$probabilidades_estrato <- probabilidades
      }
      return(tabla)
    }
    lista_estratos <- asignacion_proporcional %>% group_split(estrato) %>% 
      map(., f_seleccion) %>% reduce(., rbind)
    return(lista_estratos)
  }
  funcion_cambios_por_estratificacion <- function(datos, asignacion, 
                                                  estratos, casos_x_unidad) {
    require(stringr)
    asignacion_estratificada <- asignacion %>% group_by(estrato) %>% 
      summarise(tamano_estrato = sum(tamano_estrato), porcentaje = sum(porcentaje), 
                casos = sum(casos)) %>% mutate(casosb = if_else(str_detect(estrato, 
                                                                           "estrato"), casos_x_unidad, casos))
    diferencia <- sum(asignacion_estratificada$casos) - sum(asignacion_estratificada$casosb)
    asignacion_estratos <- asignacion_estratificada %>% filter(str_detect(estrato, 
                                                                          "estrato"))
    asignacion_municipios <- asignacion_estratificada %>% 
      filter(!str_detect(estrato, "estrato"))
    asignacion_municipios <- asignacion_municipios %>% mutate(t1 = porcentaje/sum(porcentaje), 
                                                              casos_extra = t1 * diferencia, casosb = casosb + 
                                                                casos_extra) %>% dplyr::select(-t1, -casos_extra)
    asignacion_estratificada <- rbind(asignacion_municipios, 
                                      asignacion_estratos)
    asignacion_estratificada$casos_r <- round_preserve_sum(asignacion_estratificada$casosb/10) * 
      10
    asignacion$casos_r <- w_transferencia_variable("casos_r", 
                                                   "estrato", "estrato", asignacion_estratificada, asignacion)
    asignacion <- asignacion %>% mutate(casos_r = casos_r * 
                                          inclusion) %>% data.frame
    eval(parse(text = paste0("asignacion <- asignacion %>%\n      arrange(", 
                             estratos, ")")))
    return(asignacion)
  }
  funcion_seleccion_muestra_estratos <- function(lista) {
    f_verosimilitud <- function(tabla) {
      resultado <- tabla %>% mutate(verosimilitud = inclusion * 
                                      probabilidades_estrato) %>% dplyr::select(verosimilitud) %>% 
        sum
      return(resultado)
    }
    valores <- map(lista_asignacion, f_verosimilitud) %>% 
      unlist
    mejor_muestra <- which(valores == max(valores))[1]
    resultado <- lista[[mejor_muestra]]
    return(resultado)
  }
  funcion_split_datos <- function(datos, estratos) {
    lista_marcos <- eval(parse(text = paste0("datos %>% group_split(", 
                                             paste0(estratos, collapse = ","), ")")))
    return(lista_marcos)
  }
  
  
  # tabla <- funcion_split_datos(datos = datos, 
  #                     estratos = var_estratos)[[1]]
  # 
  # casos <- asignacion$casos_r[1]
  
  
  seleccion_secciones <- function(tabla, casos, tamano, casos_seccion, 
                                  tolerancia_tamano, unidad, casos_maximos_seccion, minimo_duplicar) {
    funcion_atipicos <- function(tabla_atipicos, tamano_maximo, 
                                 casos_maximos_seccion, casos_seccion, tamano) {
      tamano_atipico <- tabla_atipicos[tamano] %>% unlist
      repeticiones <- ceiling(tamano_atipico/tamano_maximo - 
                                minimo_duplicar)
      resultado <- tabla_atipicos
      if (repeticiones >= 2) {
        for (i in 2:repeticiones) {
          resultado <- rbind(resultado, tabla_atipicos)
        }
      }
      repeticiones_maximas <- casos_maximos_seccion/casos_seccion
      resultado <- resultado[1:repeticiones_maximas, ] %>% 
        drop_na
      eval(parse(text = paste0("resultado <- resultado %>%\n    mutate(", 
                               tamano, " = tamano_atipico/nrow(resultado))")))
      return(resultado)
    }
    a <- "normal"
    if (casos == 0) {
      casos = casos_seccion
      a <- "anormal"
    }
    n <- round(casos/casos_seccion, 0)
    tamanos <- tabla[tamano] %>% unlist
    tamano_maximo <- sum(tamanos)/n
    normal <- eval(parse(text = paste0("tabla %>% filter(", 
                                       tamano, " <= tamano_maximo * ", tolerancia_tamano, 
                                       ")")))
    atipico <- eval(parse(text = paste0("tabla %>% filter(", 
                                        tamano, " > tamano_maximo * ", tolerancia_tamano, 
                                        ")")))
    lista_atipico <- eval(parse(text = paste0("atipico %>% group_split(", 
                                              unidad, ")")))
    if (length(lista_atipico) > 0) {
      lista_modificada <- lista_atipico %>% map(., funcion_atipicos, 
                                                tamano_maximo, casos_maximos_seccion, casos_seccion, 
                                                tamano) %>% reduce(., rbind) %>% rbind(., normal)
    }
    else {
      lista_modificada <- normal
    }
    tamanos <- lista_modificada[tamano] %>% unlist
    probabilidades <- inclusionprobabilities(tamanos, n)
    seleccion <- UPsystematic(probabilidades)
    if (a == "normal") {
      lista_modificada$seleccion <- seleccion
    }
    else {
      lista_modificada$seleccion <- 0
    }
    lista_modificada$probabilidades <- probabilidades
    return(lista_modificada)
  }
  resumen_muestra <- function(tabla, estrato, asignacion, unidad, 
                              tamano, casos_x_unidad) {
    resultado <- eval(parse(text = paste0("tabla %>%\n      dplyr::select(", 
                                          unidad, ",", tamano, ",seleccion,probabilidades) %>%\n      group_by(", 
                                          unidad, ") %>%\n      summarise(", tamano, " = sum(", 
                                          tamano, "),\n      casos = sum(seleccion * casos_x_unidad),\n      seleccion = max(seleccion),\n\n      probabilidades_unidad = min(sum(probabilidades),1))")))
    resultado$estrato <- estrato
    resultado <- resultado[, c(length(resultado), 1:(length(resultado) - 
                                                       1))]
    return(resultado)
  }
  pegar_probabilidades_estrato <- function(tabla, asignacion_estratos, 
                                           estratos) {
    resultado <- tabla
    eval(parse(text = paste0("resultado$probabilidades_estrato <- w_transferencia_variable(\"probabilidades_estrato\",\"", 
                             estratos, "\",\"estrato\",asignacion_estratos,tabla)")))
    names(resultado)[1] <- estratos
    eval(parse(text = paste0("resultado$estrato <- w_transferencia_variable(\"estrato\",\"", 
                             estratos, "\",\"", estratos, "\",asignacion_estratos,resultado)")))
    resultado <- resultado[, c(length(resultado), 1:(length(resultado) - 
                                                       1))]
    return(resultado)
  }
  funcion_seleccion_muestra_unidades <- function(lista) {
    f_verosimilitud <- function(tabla) {
      resultado <- tabla %>% mutate(verosimilitud = seleccion * 
                                      probabilidades_estrato * probabilidades_unidad) %>% 
        dplyr::select(verosimilitud) %>% sum
      return(resultado)
    }
    valores <- map(lista, f_verosimilitud) %>% unlist
    mejor_muestra <- which(valores == max(valores))[1]
    resultado <- lista[[mejor_muestra]]
    return(resultado)
  }
  set.seed(semilla)
  lista_asignacion <- list()
  for (i in 1:repeticiones_muestreo_estratos) {
    lista_asignacion[[i]] <- funcion_asignacion_proporcional(datos = datos, 
                                                             estratos = var_estratos, tamano = var_tamano, casos = casos) %>% 
      funcion_trim_estratos(., forzar_proporcionalidad = forzar_proporcionalidad) %>% 
      funcion_creacion_estratos(., casos_x_unidad, var_estratos, 
                                porcentaje_aceptacion_estrato_unico = porcentaje_aceptacion_estrato_unico, 
                                fusiona_estratos = fusiona_estratos) %>% funcion_seleccion_intraestrato(.) %>% 
      funcion_cambios_por_estratificacion(datos, ., estratos = var_estratos, 
                                          casos_x_unidad = casos_x_unidad)
  }
  asignacion <- lista_asignacion %>% funcion_seleccion_muestra_estratos
  lista_seleccion <- list()
  for (i in 1:repeticiones_muestreo_unidades) {
    lista_seleccion[[i]] <- funcion_split_datos(datos = datos, 
                                                estratos = var_estratos) %>% map2(., asignacion$casos_r, 
                                                                                  seleccion_secciones, tamano = var_tamano, casos_seccion = casos_x_unidad, 
                                                                                  tolerancia_tamano = tolerancia_tamano, unidad = unidad, 
                                                                                  casos_maximos_seccion = casos_maximos_seccion, minimo_duplicar = minimo_duplicar) %>% 
      map2(., asignacion[, 1], resumen_muestra, asignacion = asignacion, 
           unidad = unidad, tamano = var_tamano, casos_x_unidad = casos_x_unidad) %>% 
      reduce(., rbind) %>% pegar_probabilidades_estrato(., 
                                                        asignacion, estratos = var_estratos)
    
    
  }
  seleccion <- lista_seleccion %>% funcion_seleccion_muestra_unidades
  return(seleccion)
}