############# MUESTREO BIETAPICO ESTRATIFICADO-MAS (ESTMAS-MAS)
gc(reset=T) #Eliminar la basura
rm(list=ls(all=TRUE))

vihope <- read.csv2("F:/Una mas/vihope6.csv", sep=";")
#View(vihope)
colnames(vihope)
library(sqldf)
#Seleccionando las variables de interés
Basei=sqldf("select Mpios_n,AREA_GEOGRAFICA as ID,Tot_PerHog,
            S1_1 from vihope")
#View(Basei)

colnames(Basei)
# Número de USM
Ni=nrow(Basei)

#Sumando a nivel de UPM la variable de interés y la auxiliar
BaseI=sqldf('select Mpios_n,sum(Tot_PerHog) as Sum_Tot_PerHog,
            sum(S1_1) as SS1_1 from Basei group by Mpios_n')
#View(BaseI)
colnames(BaseI)
NI=nrow(BaseI)

library(stratification)

###### ESTRATIFICACION
set.seed(12345)
LH=strata.LH(BaseI$Sum_Tot_PerHog,CV=0.03,Ls=5)

cortes <- c(min(BaseI$Sum_Tot_PerHog),LH$bh,max(BaseI$Sum_Tot_PerHog))
BaseI$Estrato=cut(x = BaseI$Sum_Tot_PerHog, breaks = cortes, labels = NULL,
                  include.lowest = T, right = F)

BaseI$Estrato=cut(x=BaseI$Sum_Tot_PerHog, breaks = cortes, labels = c("Estrato 5","Estrato 4",
                                                                      "Estrato 3","Estrato 2","Estrato 1"),include.lowest = T, right = F)

table(BaseI$Estrato)

# Calcular el tamaño de muestra áreas geográficas

library(samplesize4surveys)
Muestra_ag <- ss4m(N = Ni, mu = mean(Basei$Tot_PerHog),sigma = mean(Basei$Tot_PerHog),
                   DEFF = 8, conf = 0.95, cve = 0.02, rme = 0.03,plot = FALSE)
tamai <- Muestra_ag$n.cve


#Seleccionar la muestra de municipios
tam_pobXest <- table(BaseI$Estrato)
tam_mueXest <- as.numeric(as.character(LH$nh))

# for repetir esto 5000 veces

m=5 
resultado=matrix(c(numeric(m),numeric(m),numeric(m),numeric(m)),ncol=4)
for(i in 1:m){
  set.seed(1000+i)
  library(sampling)
  MuestraI <- strata(data = BaseI, stratanames="Estrato", size = tam_mueXest,
                     description=F, method = "srswor")
  muestraI <- getdata(BaseI, MuestraI)
  
  
  marcomuestrai <- merge(Basei, muestraI, all.y = T, by = "Mpios_n")
  
  # SELECCIONAR LAS MUESTRAS POR ÁREA GEOGRÁFICAS
  
  # APMS es la cantidad de Àreas Geográficas por municipios seleccionados 
  APMS <-  sqldf("select Mpios_n, count(*) as 'Num_areas_geogr' 
                 from marcomuestrai
                 group by Mpios_n")
  
  
  APMS$ni <- tamai/sum(LH$nh)
  APMS$ni <- round(ifelse(APMS$ni > APMS$Num_areas_geogr,
                          APMS$Num_areas_geogr, APMS$ni ))
  
  set.seed(1000+i)
  Muestrai <- try(strata(data = marcomuestrai, stratanames="Mpios_n", size = APMS$ni,
                         description=F, method = "srswor"),silent=T)
  muesfinal <- try(getdata(marcomuestrai, Muestrai),silent=T)
  dim(muesfinal)
  
  
  
  #Verificar para garantizar que está bien: Calcular el número de municipios por estrato (Para tener los tamaños en survey)
  cons_tammueEstXMunic <- data.frame(Estrato = paste0("Estrato ", 1:5), NmunXEstrato = LH$Nh[c(5,4,3,2,1)])
  cons_tammue_mzXmpio <- APMS[c("Mpios_n", "ni")]
  muefinal1 <- try(merge(muesfinal, cons_tammueEstXMunic),silent=T) # Pegarle el tamaño de muestra por estrato
  muefinal <- try(merge(muefinal1, cons_tammue_mzXmpio),silent=T) # Pegarle el tamaño de muestra por municipio
  
  
  library(survey)
  diseno<-try(svydesign(id=~Mpios_n + ID, strata = ~  Estrato,  fpc=~NmunXEstrato +ni,
                        data = muefinal),silent=T)
  
  try(svymean(~S1_1, diseno, deff=TRUE),silent=T)
  estimacion <- try(svytotal(~S1_1, diseno, deff=TRUE),silent=T)
  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=T)
  
  
  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)
  resultado[i,4]=try(sum(APMS$ni),silent=T)
  print(paste("Muestreo",i,sep=" "))
  
  
}



write.csv2(resultado, file = 'F:/Una mas/resESTS1_1.csv2')