############# MUESTREO BIETAPICO piPT-MAS

rm(list=ls(all=TRUE))

# La base contiene el VIHOPE junto con 15 variables que se simularon altamente
# correlacionadas con el total de personas por Unidad geográfica
# 90, 70, 50, 30, 10 y variando el intercepto (porcentaje respecto a la media)
# Documentar

vihope <- read.csv2("D:/vihope5.csv", sep=";")
#View(vihope)
colnames(vihope)
library(sqldf)
Basei=sqldf("select Mpios_n,AREA_GEOGRAFICA as ID,Tot_PerHog,S2_4 from vihope")
#View(Basei)
colnames(Basei)
Ni=nrow(Basei)

# Agrupación por UPM (Municipio)
BaseI=sqldf('select Mpios_n,sum(Tot_PerHog)as Sum_Tot_PerHog,
            sum(S2_4) as SS2_4 from Basei group by Mpios_n')
#View(BaseI)
colnames(BaseI)
NI=nrow(BaseI)

library(stratification)

###### Tamaño de muestras primera etapa por estratificación
set.seed(12345)
# Se construyen 5 estratos y se calcula el tamaño de muestra por estrato de municipio
LH=strata.LH(BaseI$Sum_Tot_PerHog,CV=0.03,Ls=5)

# Calcular el tamaño de muestra para un MAS  áreas geográficas (USM)
library(samplesize4surveys)
# Se coloca un deff de tal forma que el deff del bietápico pipt - mas satisfaga el 3%
# Para el 5% serían aproximadamente 20 000 encuestas.
Muestra_ag <- ss4m(N = Ni, mu = mean(Basei$Tot_PerHog),sigma = mean(Basei$Tot_PerHog),
                   DEFF = 20, conf = 0.95, cve = 0.02, rme = 0.03, plot = FALSE)
#Cerca de 43000 encuestas para que satisfaga al final un cve del 3% en el diseño final
tamai <- Muestra_ag$n.cve
nI=sum(LH$nh)

#####' ACA EMPIEZA EL FOR MARA 5000 MUESTRAS'
# Número de simulaciones
m=5200 # Se hace un exceso de simulaciones por aquellas manzanas que tienen un indviduo.
# Inicializo una matriz para guardar el total, la varianza
resultado=matrix(c(numeric(m),numeric(m),numeric(m)),ncol=3)

for(i in 1:m){
  
  
  #### ESCOGER LOS MUNICIPIOS (MUESTRA PRIMERA ETAPA)
  
  library(TeachingSampling)
  #Se van a seleccionar municipios tomando como variable auxiliar los totales por munic (UPM)
  set.seed(1000 + i)
  MuestraIa <- S.piPS(n = nI, x = BaseI$Sum_Tot_PerHog)
  
  MuestraI <- BaseI[MuestraIa[,1],  ] 
  MuestraI$prob_inclu <- MuestraIa[,2]
  
  # SELECCIONAR LAS MUESTRAS POR ÁREA GEOGRÁFICAS (USM)
  # Basei: base de todos los registros del VIHOPE
  marcomuestrai <- merge(Basei, MuestraI[,c("Mpios_n","prob_inclu")], 
                         by = "Mpios_n", all.y = T)
  
  # Quedan sólo la tabla VIHOPE de los municipios seleccionadas
  
  # APMS es la cantidad de Àreas Geográficas (USM) por municipios seleccionados 
  APMS <-  sqldf("select Mpios_n, count(*) as 'num_mzXmpio' 
                 from marcomuestrai
                 group by  Mpios_n")
  
  APMS$ni <- tamai/sum(LH$nh) # El número de manzanas (USM) que se deben tomar por mpio
  # Si en algún municipio resulta un mayor tamaño de muestra que tamaño del mpio, 
  # Hacer censo
  APMS$ni <-  round(ifelse(APMS$num_mzXmpio > APMS$ni,
                           APMS$ni, APMS$num_mzXmpio))
  
  
  library(sampling)
  set.seed(1000 + i)
  Muestrai <- try(strata(data = marcomuestrai, stratanames="Mpios_n", size = APMS$ni,
                         description=F, method = "srswor"),silent=TRUE)
  # Extrae de la base del marco de USM de los municipios seleccionados lo resultante
  # con la función strata
  muefinal_final <- try(getdata(marcomuestrai, Muestrai),silent=TRUE)
  
  
  
  
  
  library(survey)
  diseno <-try(svydesign(id=~Mpios_n + ID, fpc=~prob_inclu +Prob, pps = "brewer",
                         data = muefinal_final),silent = TRUE)
  
  try(svymean(~S2_4, diseno, deff=TRUE),silent = TRUE)
  estimacion <- try(svytotal(~S2_4, diseno, deff=TRUE),silent=TRUE)
  # Notaciónd decimal
  options(sicipen = 1111)
  result <-try(matrix(c(ty = as.numeric(estimacion), var = (as.numeric(SE(estimacion)))^2, 
                        cve = 100 * as.numeric(SE(estimacion)) / as.numeric(estimacion)),ncol=3),silent = TRUE)
  
  
  resultado[i,1]=try(result[1,1],silent=T)
  resultado[i,2]=try(result[1,2],silent=T)
  resultado[i,3]=try(result[1,3],silent=T)
  
  print(paste("Muestreo",i,sep=" "))
  
}
write.csv2(resultado, file = 'D:/respiS2_4.csv2')