#' ---
#' title: "DOL-ILAB SDC - UC Berkeley_Nepal_Awareness-Raising for Police_Public Data_Final"
#' author: "Pablo Diego-Rosell, PhD"
#' email: "pablo_diego_rosell@yahoo.com"
#' output:
#'    html_document:
#'      toc: true
#' theme: united
#' ---

rm(list=ls(all=t))

#' #Setup and crate dictionary

filename <- "UC Berkeley_Nepal_Awareness-Raising for Police_Public Data_Final" # !!!Update filename
source ("functions_1.7.R")

#' Visually inspect variables in "dictionary.csv" and flag for risk, using the following flags: 
# Direct PII: Respondent Names, Addresses, Identification Numbers, Phone Numbers
# Direct PII-team: Interviewer Names, other field team names 
# Indirect PII-ordinal: Date of birth, Age, income, education, household composition. 
# Indirect PII-categorical: Gender, education, ethnicity, nationality,
# occupation, employer, head of household, marital status
# GPS: Longitude, Latitude
# Small Location: Location (<100,000) 
# Large Location (>100,000)
# Weight: weightVar
# Household ID:  hhId, 
# Open-ends: Review responses for any sensitive information, redact as necessary 

#' #Direct PII: variables to be removed
# !!!No Direct PII variables
#'
#' #Direct PII-team: Encode interviewer names, which may be useful for analysis of interviewer effects
#' !!!Replace vector in "variables" field below with relevant variable names

mydata <- encode_direct_PII_team (variables=c("Srvyr"))

#' #Small locations: Encode locations  with pop <100,000 using random large numbers
#'  !!!Include relevant variables, but check their population size first to confirm they are <100,000

locvars <- c("vdc_id", "police_id") 
mydata <- encode_location (variables= locvars, missing=999999)

#' #Indirect PII - Ordinal: Global recode or Top/bottom coding for extreme values
# Recode birth year into 5 year groups.

break_year <- c(999,2020,2025,2030,2035,2040,2045,2050, 2055)
labels_year <- c("Missing"  =1,
                "2020-2024" =2, 
                "2025-2029" =3, 
                "2030-2034" =4, 
                "2035-2039" =5,
                "2040-2044" =6,
                "2045-2049" =7,
                "2050-2054" =8)
mydata <- ordinal_recode (variable="dem4", break_points=break_year, missing=999999, value_labels=labels_year)

# Recode education attainment into standard Nepal categories

break_edu <- c(0,7,10,12,13,14,16,17,18, 777, 888, 999)
labels_edu <- c("Primary or less (0-5)" = 1,
                "Lower secondary (6-8)" = 2,
                "Secondary (9-10)" = 3,
                "SLC (11)" = 4,
                "CLASS 12/Intermediate level (12)" = 5,
                "Bachelor/Postgraduate level" = 6,
                "Literate, but never attended school" = 7,
                "Illiterate, and never attended school"= 8,
                "Refused"= 9, 
                "Does not apply" = 10, 
                "Don't Know" = 11)
mydata <- ordinal_recode (variable="dem6", break_points=break_edu, missing=999, value_labels=labels_edu)

mydata <- top_recode ("dem9", break_point=2, missing=999999) # Topcode cases with 2 or more children, as there are only 16 cases with 3 or more. 
mydata <- top_recode ("dem10", break_point=2, missing=999999) # Topcode cases with 2 or more children, as there are only 36 cases with 3 or more. 

# Top code high income to the 99.5 percentile

percentile_99.5 <- floor(quantile(mydata$dem11[mydata$dem11>999], probs = c(0.995)))
mydata <- top_recode (variable="dem11", break_point=percentile_99.5, missing=999)

# Group year of employment into 5-year groups. 

break_year <- c(0, seq(2050, 2075, 5))
labels_year <- c("Before 2050"= 1,
                 "2050-2054" = 2, 
                 "2055-2059" = 3, 
                 "2060-2064" = 4, 
                 "2065-2069" = 5, 
                 "2070-2074" = 6)
mydata <- ordinal_recode (variable="emp1y", break_points=break_year, missing=999999, value_labels=labels_year)

# !!!Include relevant variables in list below

indirect_PII <- c("list_id", 
                  "police_rank", 
                  "dem3", 
                  "dem5", 
                  "dem7", 
                  "dem8", 
                  "emp20_1", 
                  "emp20_2", 
                  "emp20_5", 
                  "govlaw0", 
                  "enumimp1_a", 
                  "enumimp1_b", 
                  "enumimp1_c")

capture_tables (indirect_PII)

# Recode those with very specific values where more than half of the sample have actual data. 

# Top code high police ranks, as only a few "Superintendent of Police (SP)", "Deputy Superintendent of Police (DSP)", or "Inspector\t"

break_rank <- c(0, 4, 5, 6, 7,8)
labels_rank <- c("Superintendent, Deputy Superintendent, or Inspector"= 1,
                 "Sub Inspector" = 2, 
                 "Assistant Sub Inspector" = 3, 
                 "Head Constable" = 4, 
                 "Constable" = 5)
mydata <- ordinal_recode (variable="police_rank", break_points=break_rank, missing=999999, value_labels=labels_rank)
mydata <- ordinal_recode (variable="enumimp1_c", break_points=break_rank, missing=999999, value_labels=labels_rank)

# Encode caste and political party, as many small groups

mydata <- encode_location (variables= "dem5", missing=999999)
mydata <- encode_location (variables= "govlaw0", missing=999999)

# Recode religion

break_rel <- c(1,2,3, 777, 888, 999)
labels_rel <- c("Hindu" = 1,
                "Buddhist" = 2,
                "Other" = 3,
                "Refused" = 4,
                "Not applicable" = 5,
                "Don't know" = 6)
mydata <- ordinal_recode (variable="dem7", break_points=break_rel, missing=999, value_labels=labels_rel)

# Recode marital status

break_mar <- c(1,2,777, 888, 999)
labels_mar <- c("Married" = 1,
                "Other" = 2,
                "Refused" = 3,
                "Not applicable" = 4,
                "Don't know" = 5)
mydata <- ordinal_recode (variable="dem8", break_points=break_mar, missing=999, value_labels=labels_mar)

#' #Matching and crosstabulations: Run automated PII check 

# Based on dictionary inspection, select variables for creating sdcMicro object
# See: https://sdcpractice.readthedocs.io/en/latest/anon_methods.html
# All variable names should correspond to the names in the data file
# selected categorical key variables: gender, occupation/education and age
selectedKeyVars = c('dem3', 'police_rank', "dem4") ##!!! Replace with candidate categorical demo vars
# weight variable
selectedWeightVar = c('weight') ##!!! Replace with weight var
# creating the sdcMicro object with the assigned variables
sdcInitial <- createSdcObj(dat = mydata, keyVars = selectedKeyVars, weightVar=selectedWeightVar)
sdcInitial

#' Show values of key variable of records that violate k-anonymity
#mydata <- labelDataset(mydata)
notAnon <- sdcInitial@risk$individual[,2] < 2 # for 2-anonymity
mydata[notAnon,selectedKeyVars]
sdcFinal <- localSuppression(sdcInitial)

# Recombining anonymized variables

extractManipData(sdcFinal)[notAnon,selectedKeyVars] # manipulated variables HH
mydata [192,"police_rank"] <- NA
mydata [192,"enumimp1_c"] <- NA
mydata [notAnon,"dem4"] <- NA
mydata [192,"dem4"] <- 3

# Check that 2-anonimity has been achieved. 

sdcInitial <- createSdcObj(dat = mydata, keyVars = selectedKeyVars, weightVar=selectedWeightVar)
sdcInitial

#' #Open-ends: review responses for any sensitive information, redact as necessary
# !!! Identify open-end variables here: 

open_ends <- c("start_failure_reason_other",
               "start_failure_reason_other_eng",
               "mc1_other",
               "mc1_other_eng",
               "ht25o",
               "ht25o_eng",
               "ht26o",
               "ht26o_eng",
               "govlaw0_other",
               "govlaw0_other_eng",
               "enumimp9_other",
               "enumimp7o",
               "enumimp7o_eng",
               "enumimp10",
               "enumimp10_eng",
               "poster_read")

report_open (list_open_ends = open_ends)

# Review "verbatims.csv". Identify variables to be deleted or redacted and their row number 

mydata <- mydata[!names(mydata) %in% c("govlaw0_other","govlaw0_other_eng")] # Specific political party mentioned

#' #GPS data: Displace
# Setup map

countrymap <- map_data("world") %>% filter(region=="Nepal")  #!!! Select correct country
#admin <- raster::getData("GADM", country="NP", level=0) #!!! Select correct country map using standard 2-letter country codes: https://en.wikipedia.org/wiki/ISO_3166-1_alpha-2
admin <- readRDS(file="gadm36_NPL_0_sp.rds")

# Displace all pairs of GPS variables (Longitude, Latitude). Check summary statistics and maps before and after displacement. 

gps.vars <- c("Longitude", "Latitude") # !!!Include relevant variables, always longitude first, latitude second.
mydata <- displace(gps.vars, admin=admin, samp_num=1, other_num=100000) # May take a few minutes to process.

gps.vars <- c("GPSinitial_LO", "GPSinitial_LA") # !!!Include relevant variables, always longitude first, latitude second.
mydata <- displace(gps.vars, admin=admin, samp_num=1, other_num=100000) # May take a few minutes to process.

gps.vars <- c("gpsht13a_LO", "gpsht13a_LA") # !!!Include relevant variables, always longitude first, latitude second.
mydata <- displace(gps.vars, admin=admin, samp_num=1, other_num=100000) # May take a few minutes to process.

gps.vars <- c("gpsenumip_LO", "gpsenumip_LA") # !!!Include relevant variables, always longitude first, latitude second.
mydata <- displace(gps.vars, admin=admin, samp_num=1, other_num=100000) # May take a few minutes to process.

#' #Save processed data in Stata and SPSS format
#' Adds "_PU" (Public Use) to the end of the name 

haven::write_dta(mydata, paste0(filename, "_PU.dta"))
haven::write_sav(mydata, paste0(filename, "_PU.sav"))
