PES Content Error

recodeST = function(st){
  st = st %>% 
    mutate(
      device = ifelse(sumber == "PAPI-Desktop", "D", ifelse(sumber == "PAPI-WebEntry", "W", NA))
    )  %>% 
    mutate(
      id_ruta_st = ifelse(sumber == "CAPI", id_ruta, paste0(id, device)),
      
      b2r202 = as.numeric(gsub(".*#(. *)", "\\1", id_art)),
      id_art_st = ifelse(sumber == "CAPI", id_art, as.character(b2r202))
    ) 
  st
}

expandDataset= function(pair,st,pes,group, varnames = "var", dropna = T,num_to_bin = 0,numeric = FALSE, matching_only = F){
  colnames(st) = c("id_1","id_2","var_st")
  colnames(pes) = c("id_1","id_2","var_pes")
  contentset = pair %>% left_join(st,by= c("id_ruta_st" = "id_1" ,"id_art_st"= "id_2" ))
  contentset = contentset %>% left_join(pes,by= c("id_ruta_pes" = "id_1" ,"id_art_pes"= "id_2" ))
  
  if(matching_only) return ( contentset)
  #Na treatment
  if (dropna ) {
    contentset = contentset  %>% filter(!is.na(var_st)| !is.na(var_pes))
    contentset[is.na(contentset)] = 0
    
  }else{
    contentset[is.na(contentset)] = 0
  }
  
  if(numeric) {return (contentset %>% select_at(c(group,"var_st","var_pes")))}
  
  #convert numeric to binary
  if(num_to_bin != 0 ){
    
    contentset = contentset %>% mutate(growth = 100*(var_pes-var_st)/var_st,
                                       growth = if_else(is.na(growth),0,growth),
                                       code = "1#1",
                                       code = ifelse(var_st>var_pes  & growth< -num_to_bin,"1#0",code),
                                       code = ifelse(var_st<var_pes & growth > num_to_bin,"0#1",code),
                                       code = ifelse(growth <0 & growth > -num_to_bin,"0#0",code)
    )%>% select(-var_st,-var_pes,-growth) 
    
    contentset = contentset  %>% separate_wider_delim(code,"#",names = c("var_st","var_pes")) 
  } 
  
  ctb = contentset %>% group_by_at(c(group,"var_st","var_pes")) %>% summarise(cnt = n()) %>% ungroup()
  rlevel = unique(c(ctb$var_st,ctb$var_pes))
  grid = expand.grid(unique(ctb$kdprov),rlevel,rlevel)
  colnames(grid) = colnames(ctb)[-NCOL(ctb)]
  mytab  = merge(grid,ctb,all.x = TRUE) 
  colnames(mytab) = c(group,paste0(varnames,"_st"),paste0(varnames,"_pes"),"cnt")
  mytab
}


find_mode <- function(x) {
  x = as.numeric(x)
  u <- sort(unique(x))
  tab <- tabulate(match(x, u))
  u[tab == max(tab)][1]
}




computeContent = function(dataset,freq,cencus_var,pes_var, group = "None" ){
  
  dataset = dataset %>% rename("Freq" = freq)
  if(group == "None"){
    dataset$group = "Nasional"
    group = "group"
  }
  #Produce Cross Tab
  ctab = dataset %>% group_by_at(c(group,cencus_var,pes_var)) %>% summarise_at("Freq",sum,na.rm=TRUE) %>% ungroup()  %>% spread(key = pes_var,value = Freq)
  ctab$total = rowSums(ctab[,-c(1,2)])
  idx = colnames(ctab)[-c(1:2)]
  
  coltot = ctab %>% group_by_at(c(group)) %>% summarise_at(idx,sum,na.rm=T)
  ctab = bind_rows(ctab, coltot) %>% arrange_at(c(group,cencus_var)) %>% mutate_at(cencus_var ,as.character)
  ctab[is.na(ctab)] = "total"
  
  #Prepare variables
  vrt = dataset %>% group_by_at(c(group,cencus_var)) %>% summarise_at("Freq",sum,na.rm=TRUE) %>% ungroup() %>% rename(vrt = Freq)
  hrz = dataset %>% group_by_at(c(group,pes_var)) %>% summarise_at("Freq",sum,na.rm=TRUE)  %>% ungroup()%>% .$Freq
  dig = dataset %>% filter(eval(parse(text = paste0(cencus_var,"==",pes_var)))) %>%group_by_at(c(group,pes_var)) %>% summarise_at("Freq",sum,na.rm=TRUE) %>%  ungroup()%>% .$Freq
  
  #produce basic indicators
  dataset = vrt %>%  mutate(vrt = as.numeric(vrt),hrz =as.numeric(hrz),dig = as.numeric(dig),vh = vrt*hrz  ) %>%
    group_by_at(group) %>% mutate(n = sum(vrt)) %>% ungroup()
  dataset = dataset %>% mutate(hrz =hrz ) %>% group_by_at(group) %>% mutate(n = sum(vrt)) %>% ungroup()
  dataset = dataset %>% mutate(NDR = ((vrt-hrz)/n)*100,
                               IOI = (vrt+hrz-2*dig)/((1/n)*(vrt*(n-hrz)+hrz*(n-vrt)))*100
  )
  
  #produce aggregate indicators
  datasetAgg = dataset %>% group_by_at(c(group,"n")) %>% summarise(across(vrt:vh,sum))
  datasetAgg = datasetAgg %>% mutate(AggIOI = (n-dig)/(n-(1/n)*vh),
                                     GDR = (n-dig)/(n-(1/n)*dig),
                                     ROA = dig/n
  )
  return(list(ctab = ctab,basic = dataset, agg =datasetAgg))
}

computeContentReg = function(dataset,cencus_var,pes_var, group = "None" ){
  
  
  if(group == "None"){
    dataset$group = "Nasional"
    group = "group"
  }
  
  res = dataset %>%  group_by_at(group) %>% summarise(Jumlah_observasi = n(),
                                                      Mean_ST =  mean(var_st,na.rm = T) ,
                                                      Mean_PES = mean(var_pes,na.rm = T) ,
                                                      MSE = mse(var_st,var_pes),
                                                      RMSE = rmse(var_st,var_pes),
                                                      R2 = R2_Score(var_st,var_pes),
                                                      MAE = mae(var_st,var_pes),
                                                      SMAPE = smape(var_st,var_pes)
                                                      )
  res
}

showres = function(df){
  df %>% kbl() %>% kable_classic_2(full_width=F)
}


exportRegression = function(nasional,provinsi, plot_reg ,var,folder=""){
  level = "Nasional"
  level2 = "Provinsi"
  title = paste0("Tabel Hasil Penghitungan Content Error ",var, " -Regresi")
  title1 = paste0("Tabel Hasil Penghitungan Content Error ",var," tingkat ",level)
  sheetname = "Nasional"
  wb=createWorkbook()
  addWorksheet(wb, sheetName = sheetname)
  writeData(wb, sheet = sheetname, x = title1, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname, x = nasional, startRow = 3)
  s2 = createStyle(border = "TopBottomLeftRight")
  s1 = createStyle( textDecoration = c("BOLD"))
  addStyle(wb, sheet = sheetname, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname, style = s2, rows = 3:(3+NROW(nasional)), cols = 1:NCOL(nasional), gridExpand = TRUE)
  
  title2 = paste0("Tabel Hasil Penghitungan Content Error ",var," tingkat ",level2)
  sheetname2 = "Provinsi"
  addWorksheet(wb, sheetName = sheetname2)
  writeData(wb, sheet = sheetname2, x = title2, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname2, x = provinsi, startRow = 3)
  addStyle(wb, sheet = sheetname2, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname2, style = s2, rows = 3:(3+NROW(provinsi)), cols = 1:NCOL(provinsi), gridExpand = TRUE)
  
  
  sheetname3 = "Plot"
  title3 = paste0("Diagram Pencar ",var)
  addWorksheet(wb, sheetName = sheetname3)
  writeData(wb, sheet = sheetname3, x = title3, startRow = 1, startCol = 1)
  addStyle(wb, sheet = sheetname3, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  print(plot_reg +xlab(paste(word(var,1),"pada","ST")) +ylab(paste(word(var,1),"pada","PES"))) 
  insertPlot(wb, sheet = sheetname3, fileType = "png", units = "in",startRow = 3 ,startCol = 2)
  
  saveWorkbook(wb, paste0(folder,Sys.Date()," ",title,".xlsx"))
  print("Cie ngesavenya berhasil !")
}

exportCategoric = function(nasional, provinsi,var,folder=""){
  level = "Nasional"
  level2 = "Provinsi"
  title = paste0("Tabel Hasil Penghitungan Content Error ",var," -Kategorik")
  title1 = paste0("Tabel Tabulasi ",var," tingkat ",level)
  
  sheetname = "Nasional Tabulasi"
  wb=createWorkbook()
  addWorksheet(wb, sheetName = sheetname)
  writeData(wb, sheet = sheetname, x = title1, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname, x = nasional$ctab, startRow = 3)
  s2 = createStyle(border = "TopBottomLeftRight")
  s1 = createStyle( textDecoration = c("BOLD"))
  addStyle(wb, sheet = sheetname, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname, style = s2, rows = 3:(3+NROW(nasional$ctab)), cols = 1:NCOL(nasional$ctab), gridExpand = TRUE)
  
  
  sheetname2 = "Nasional Kelas"
  title2 = paste0("Tabel Indikator NDR dan IOI untuk ",var," tingkat ",level)
  addWorksheet(wb, sheetName = sheetname2)
  writeData(wb, sheet = sheetname2, x = title2, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname2, x = nasional$basic, startRow = 3)
  addStyle(wb, sheet = sheetname2, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname2, style = s2, rows = 3:(3+NROW(nasional$basic)), cols = 1:NCOL(nasional$basic), gridExpand = TRUE)
  
  sheetname3 = "Nasional Agregat"
  title3 = paste0("Tabel Indikator Agregat IOI,GDR, dan ROA untuk ",var," tingkat ",level)
  addWorksheet(wb, sheetName = sheetname3)
  writeData(wb, sheet = sheetname3, x = title3, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname3, x = nasional$agg, startRow = 3)
  addStyle(wb, sheet = sheetname3, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname3, style = s2, rows = 3:(3+NROW(nasional$agg)), cols = 1:NCOL(nasional$agg), gridExpand = TRUE)
  
  level = level2
  sheetname = "Provinsi Tabulasi"
  title1 = paste0("Tabel Tabulasi ",var," tingkat ",level)
  addWorksheet(wb, sheetName = sheetname)
  writeData(wb, sheet = sheetname, x = title1, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname, x = provinsi$ctab, startRow = 3)
  addStyle(wb, sheet = sheetname, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname, style = s2, rows = 3:(3+NROW(provinsi$ctab)), cols = 1:NCOL(provinsi$ctab), gridExpand = TRUE)
  
  sheetname2 = "Provinsi Kelas"
  title2 = paste0("Tabel Indikator NDR dan IOI untuk ",var," tingkat ",level)
  addWorksheet(wb, sheetName = sheetname2)
  writeData(wb, sheet = sheetname2, x = title2, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname2, x = provinsi$basic, startRow = 3)
  addStyle(wb, sheet = sheetname2, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname2, style = s2, rows = 3:(3+NROW(provinsi$basic)), cols = 1:NCOL(provinsi$basic), gridExpand = TRUE)
  
  sheetname3 = "Provinsi Agregat"
  title3 = paste0("Tabel Indikator Agregat IOI,GDR, dan ROA untuk ",var," tingkat ",level)
  addWorksheet(wb, sheetName = sheetname3)
  writeData(wb, sheet = sheetname3, x = title3, startRow = 1, startCol = 1)
  writeData(wb, sheet = sheetname3, x = provinsi$agg, startRow = 3)
  addStyle(wb, sheet = sheetname3, style = s1, rows = 1, cols = 1, gridExpand = TRUE)
  addStyle(wb, sheet = sheetname3, style = s2, rows = 3:(3+NROW(provinsi$agg)), cols = 1:NCOL(provinsi$agg), gridExpand = TRUE)
  
  saveWorkbook(wb, paste0(folder,Sys.Date()," ",title,".xlsx"))
  print("Cie ngesavenya berhasil !")
}

library("rio")
library("dplyr")
library("tidyr")
library("Metrics")
library("MLmetrics")
library("kableExtra")
library(ggplot2)
library(openxlsx)
library(stringr)
#Load data pairs dan add kode provinsi
pairs= import("2024-01-15 Data Pairs Perbaikan Data Terbaru.xlsx")
pairs = pairs %>% select(idsubsls,id_ruta_pes,id_art_pes,id_ruta_st,id_art_st)
pairs$kdprov = substr(pairs$idsubsls,1,2)

pes_btiga = import("Data/SIS/PES/Data3_1_202312121016.xlsx")
#pes_btiga_2 = import("Data/SIS/PES/Data3_202312121016.xlsx")


# Load data ST dan recode ST ID ke ID Pairs
st_btiga  = import("ST 2023-01-14 FINAL/l2_blok3.Rds")
st_btiga  = recodeST(st_btiga )

tabulasi_numeric = expandDataset(pairs,st,pes,"kdprov",numeric = T)

resnas_310 = computeContentReg(tabulasi_numeric,"var_st","var_pes")
resnas_310  %>% kbl()%>%  kable_classic_2(full_width=F)

resprov_310 = computeContentReg(tabulasi_numeric,"var_st","var_pes",group = "kdprov")
resprov_310  %>%head() %>% kbl()%>%  kable_classic_2(full_width=F)

plot_310 = ggplot(tabulasi_numeric, aes(x=var_st, y= var_pes)) + 
  geom_point(shape=18, color="blue")+
  geom_smooth(method=lm,
             color="darkred", fill="blue") 
plot_310

exportRegression(resnas_310,resprov_310,plot_310,var = "B3R10 (Jumlah Bidang)")

tabulasi_recat = expandDataset(pairs,st,pes,"kdprov",num_to_bin = 30)
tabulasi_recat %>%head() %>% kbl()%>%  kable_classic_2(full_width=F)

resnas_310 = computeContent(tabulasi_recat,"cnt","var_st","var_pes")
resprov_310 = computeContent(tabulasi_recat,"cnt","var_st","var_pes",group = "kdprov")

exportCategoric(resnas_310,resprov_310,var = "B3R10 (Jumlah Bidang)")