assignment.Rmd

---
title: "Analysis"
author: "Yue Xiong"
date: '2022-09-15'
output:
  pdf_document: default
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## Step1: Data Preparation

In step1, we aim to read in the original datasets and the synthetic dataset while aligning the variables especially set in the syn dataset. Attempting to compare the two datasets, we should also concatenate the original datasets vertically with variables all set as the same. Also, we are filtering out all the non-GPDR countries.

```{r echo=FALSE, include=FALSE}
# install the required R packages
# install.packages("readr")
# install.packages("vroom")
# install.packages("tidyverse")
# install.packages("arsenal")
# install.packages("reshape2")
# install.packages("synthpop")
library(readr)
library(vroom)
library(tidyverse)
library(arsenal)
library(reshape2)
library(synthpop)
library(ggplot2)
library(dbplyr)
```

```{r}
# set the working directory
wd <- getwd()
setwd(wd)
# read in the synthetic dataset
syn_data <- read.csv(file = "./syn_2020-08-02_2020-08-08.csv")
# and on macos 
head(syn_data)
# syn_data <- read_csv("./syn_2020-08-02_2020-08-08.csv", show_col_types = FALSE)
# rename "sample weight" to "weight" to avoid conflicts
colnames(syn_data)[colnames(syn_data) == "sample_weight"] <- "weight"

# columns to be included from the original dataset
cols_list <- colnames(syn_data)

# initialize an empty dataset list
ori_dataset <- list()

for (i in 1:7){
ori_dataset[[i]] <- vroom(list.files(pattern = "*_full.csv$")[i],
                    show_col_types = FALSE) %>%
                    select(all_of(cols_list))
}

dim(ori_dataset[[2]])[2] == dim(syn_data)[2]
# check whether 2 dimensions coincide with each other

# bind the original datasets from 0802 to 0808 vertically
bindori_dataset <- bind_rows(ori_dataset)
bindori_dataset <- as.data.frame(bindori_dataset)
dim(bindori_dataset)
dim(syn_data)
```
Now we need to check all the `gpdr` countries and only do the alignment based on `gpdr` countries.
```{r}
gpdr_countries_data <- NA
gpdr_countries_data <- read.csv(file = "./gpdr.csv", sep = ",")
head(gpdr_countries_data$Country_GID, n = 10L)
```
We then filter the binded original datasets and syn dataset with `Region_GID` specified
```{r}
country_name <- unique(as.character(gpdr_countries_data$Country_GID))
length(country_name)
length(unique(bindori_dataset$GID_0))
bindori_dataset_filtered <- bindori_dataset %>%
                                        filter(GID_0 %in% country_name)


syn_data_filtered <- syn_data %>%
                          filter(GID_0 %in% country_name)
length(unique(bindori_dataset_filtered$GID_0)) == length(unique(syn_data_filtered$GID_0))  # check whether we have the same number of countries for both synthetic and original datasets.
```
## Step2: make all the vars aligned with each other
These data that need such alignment are the first filtered ones, `syn_data_filtered` and `bindori_dataset_filtered`. They are filtered only with those gpdr countries filtered. First of all, we need to output the types of all the variables included in the two datasets and compare the difference at the first place.

```{r echo=FALSE, include=FALSE}
bindori_dataset_filtered %>%
                    summarise_all(class) %>%
                    gather
```
```{r echo=FALSE, include=FALSE}
syn_data_filtered %>%
              summarise_all(class) %>%
              gather
```
Apparently shown above, one thing to find is that all the previously listed `double`-type variables, such as the `E4` variable. We should firstly collect all these variables and change the type of the original variables.

### 2.1 numeric to integer
```{r}
vars_list <- NULL
int_vars_list <- c()

for (i in 1:ncol(syn_data_filtered)) {
  if(class(syn_data_filtered[[i]])=="integer") {
    # print(i)
    int_vars_list <- append(int_vars_list, colnames(syn_data_filtered[i]))
  }
}

length(int_vars_list)
```
So let us convert all the numeric variables to integer.
```{r}
bindori_dataset_filtered[int_vars_list] <- sapply(bindori_dataset_filtered[int_vars_list], as.integer)

# check whether all the previous numeric types are converted to integer
sapply(bindori_dataset_filtered[int_vars_list], class)
```
### 2.2 var B2
Firstly, we change the remaining "B2", "B4", "E5", "E6" to character type in the original dataframe.

For variable B2, let's take a look at the table of the occurrences of indicators.
```{r echo=FALSE, include=FALSE}
table(bindori_dataset_filtered$B2)
```
And for the synthetic dataset, we observed that it is the "character" type.
```{r echo=FALSE, include=FALSE}
table(syn_data_filtered$B2)
```
In order to align with the synthetic dataset types, we make all the cell values corresponded to the the value/threshold shown above.
```{r echo=FALSE, include=FALSE}
table(bindori_dataset_filtered["B2"])
```

```{r}
bindori_dataset_filtered["B2"][bindori_dataset_filtered["B2"] == -99] <- "-99"
```

```{r}
bindori_dataset_filtered["B2"][bindori_dataset_filtered["B2"] >= 1000] <- "1000"
```
```{r}
bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 0 & bindori_dataset_filtered$B2 < 1] <- "[0, 1)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 1 & bindori_dataset_filtered$B2 < 3] <- "[1, 3)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 3 & bindori_dataset_filtered$B2 < 8] <- "[3, 8)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 8 & bindori_dataset_filtered$B2 < 15] <- "[8, 15)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 15 & bindori_dataset_filtered$B2 < 28] <- "[15, 28)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 28 & bindori_dataset_filtered$B2 < 90] <- "[28, 90)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 90 & bindori_dataset_filtered$B2 < 180] <- "[90, 180)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 180 & bindori_dataset_filtered$B2 < 366] <- "[180, 366)"

bindori_dataset_filtered$B2[bindori_dataset_filtered$B2 >= 366 & bindori_dataset_filtered$B2 < 1000] <- "[366, 1000)"

```
```{r}
# quickly check the class
class(bindori_dataset_filtered$B2)
```
### 3.3 var B4
```{r echo=FALSE, include=FALSE}
table(syn_data_filtered$B4)
```
```{r}
bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 < 0] <- "-99"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 1000] <- "1000"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 0 & bindori_dataset_filtered$B4 < 1] <- "[0, 1)"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 1 & bindori_dataset_filtered$B4 < 5] <- "[1, 5)"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 5 & bindori_dataset_filtered$B4 < 10] <- "[5, 10)"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 10 & bindori_dataset_filtered$B4 < 1000] <- "[10, 1000)"
```
```{r}
# quickly check the class
class(bindori_dataset_filtered$B4)
```
### 3.4 var E5
```{r echo=FALSE, include=FALSE}
table(syn_data_filtered$E5)
```
```{r}
bindori_dataset_filtered$E5[bindori_dataset_filtered$E5 < 0] <- "-99"

bindori_dataset_filtered$E5[bindori_dataset_filtered$E5 >= 1000] <- "1000"

bindori_dataset_filtered$E5[bindori_dataset_filtered$E5 >= 0 & bindori_dataset_filtered$E5 < 1] <- "[0, 1)"

bindori_dataset_filtered$E5[bindori_dataset_filtered$E5 >= 1 & bindori_dataset_filtered$E5 < 2] <- "[1, 2)"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 2 & bindori_dataset_filtered$B4 < 4] <- "[2, 4)"

bindori_dataset_filtered$B4[bindori_dataset_filtered$B4 >= 4 & bindori_dataset_filtered$B4 < 6] <- "[4, 6)"

bindori_dataset_filtered$E5[bindori_dataset_filtered$E5 >= 6 & bindori_dataset_filtered$E5 < 1000] <- "[6, 1000)"
```
```{r}
# quickly check the class
class(bindori_dataset_filtered$E5)
```
### 3.5 var E6
```{r echo=FALSE, include=FALSE}
table(syn_data_filtered$E6)
```
```{r}
bindori_dataset_filtered$E6[bindori_dataset_filtered$E6 < 0] <- "-99"

bindori_dataset_filtered$E6[bindori_dataset_filtered$E6 >= 26] <- "26"

bindori_dataset_filtered$E6[bindori_dataset_filtered$E6 >= 0 & bindori_dataset_filtered$E6 < 9] <- "[0, 9)"

bindori_dataset_filtered$E6[bindori_dataset_filtered$E6 >= 9 & bindori_dataset_filtered$E6 < 26] <- "[9, 26)"
```
```{r}
# quickly check the class
class(bindori_dataset_filtered$E6)
```
### 3.6 check for observations for var D1 and D2(ori)
For D1, we have
```{r}
table(bindori_dataset_filtered$D1)
table(syn_data_filtered$D1)
```
For D2, we have
```{r}
table(bindori_dataset_filtered$D2)
table(syn_data_filtered$D2)
```

## Step3: Evaluating the utility of the syn data
In step3, we try to evaluate the utility of the synthetic dataset with one-way marginal and two-way marginal measures. 

- As for the one-way utility, the syn dataset is measured with the compare plots and pMSE/S_pMSE.

- And for the two-way utility, the synthetic dataset is evaluated with the utility tables which takes up a heatmap fashion/manner.

```{r}
# # subset some example columns and try plotting the compare histograms
# # these are symptoms variables
# symptoms <- c("B3","B4","B1_1","B1_2","B1_3","B1_4","B1_5","B1_6","B1_7",
#               "B1_8","B1_9","B1_10","B1_12","B1_13",
#               "B1b_x1","B1b_x2","B1b_x3","B1b_x4","B1b_x5","B1b_x6","B1b_x7",
#               "B1b_x8","B1b_x9","B1b_x10","B1b_x12","B1b_x13","B2")
# 
# # these are testing variables
# testing <- c("B7")
# 
# # Columns `B0`, `B8a`, `B15_1`, `B15_2`, `B15_3`, etc. don't exist
# 
# ori_dataset_symptoms <- as.data.frame(bindori_dataset_filtered[, symptoms])
# syn_dataset_symptoms <- as.data.frame(syn_data_filtered[, symptoms])
# 
# ori_dataset_testing <- as.data.frame(bindori_dataset_filtered$B7)
# syn_dataset_testing <- as.data.frame(syn_data_filtered$B7)
```

### (1). one-way marginals using compare()
First, we need to filter out variables `GID_0` and `GID_1` as they simply stand for the country/region. We do the operations both for the synthetic and ordinary datasets.
```{r}
syn_data_filtered2 <- subset(syn_data_filtered, select = -c(GID_0, GID_1))

bindori_dataset_filtered2 <- subset(bindori_dataset_filtered, select = -c(GID_0, GID_1))
```
Now we can see there are 90 variables remained in the datasets.
```{r}
ncol(syn_data_filtered2)
ncol(bindori_dataset_filtered2)
```
We need to filter out B2,B4,E5,E6.
```{r}
# try with 86 vars firstly
syn_select_vars <- subset(syn_data_filtered2, select = -c(B2, B4, E5, E6, D6_2))

bindori_select_vars <- subset(bindori_dataset_filtered2, select = -c(B2, B4, E5, E6, D6_2))
```

```{r}
compare_plots<- c()

for (i in 1:85) {
  cat(colnames(bindori_select_vars[i]), "\n")  # print the var string under analysis
  
  compare_plots[[i]] <- compare(object = data.frame(Pdata = syn_select_vars[i]),
          data = data.frame(Pdata = bindori_select_vars[i]),
          vars = c(colnames(bindori_select_vars[i])), cont.na = NULL,
          msel = NULL, stat = "percents", breaks = 10,
          nrow = 2, ncol = 2, rel.size.x = 1,
          utility.stats = c("pMSE", "S_pMSE"),
          cols = c("#1A3C5A","#4187BF"),
          plot = TRUE, table = TRUE)
  
}
```
Now we need to make all of these compare plots into a pdf file.
```{r}
# specify the file path to store the pdf
destination_path <- "F:/Master-Thesis-DifferentialPrivacy/Figures/one_way_compare.pdf"

# Print plots to a pdf file
pdf(destination_path)

for (i in 1:85) {
  print(compare_plots[[i]]$plots)  # Plot 1 --> in the first page of PDF
}

dev.off()
```

With PNGs,
```{r}
# specify the file path to store the PNG files
destination_path <- "F:/Master-Thesis-DifferentialPrivacy/Figures/one_way_"

# Print plots to a pdf file


for (i in 1:85) {
  png(paste0(destination_path, compare_plots[[i]]$vars, ".png"))
  print(compare_plots[[i]]$plots)  # Plot 1 --> in the first page of PDF
  dev.off()
}


```
We then try exporting the __tab_utility__ as a csv file in convenience of comparing and choose vars which performed better in the synthesis.

First of all, we need to append each of the instances based on tab_utility to dataframes, respectively.
```{r}
pMSE_list <- c()
S_pMSE_list <- c()
for (i in 1:85) {
  pMSE_list <- append(pMSE_list, compare_plots[[i]]$tab.utility[1])
  S_pMSE_list <- append(S_pMSE_list, compare_plots[[i]]$tab.utility[2])
  }
```
```{r}
length(pMSE_list)
length(S_pMSE_list)
```

Then, we create a dataframe.
```{r}
#create data frame
df_utility <- data.frame(vars_list=colnames(syn_select_vars),
                 pMSE=pMSE_list,
                 S_pMSE=S_pMSE_list)
```
```{r}
write_csv_dest <- "F:/Master-Thesis-DifferentialPrivacy/Tab_Utility/Terrence_oneway_utility.csv"
write.csv(df_utility, write_csv_dest, row.names=FALSE)
```
We need to output those variables whose *S_pMSE* is smaller than 10.


```{r}
vars2show <- df_utility[df_utility[, "S_pMSE"]<10, ][1]
```
```{r}
nrow(vars2show)  # there are 37 in total
```
Apart from that, we noticed that for var `D9`, the corresponding statistic of 
interest has returned a NA value. We then need to redo the evaluation again
just for `D9`.
```{r}
# we try to get the distribution of D9
table(bindori_dataset_filtered2$D9)
table(syn_data_filtered2$D9)
```



```{r}
compare(object = data.frame(D9 = syn_select_vars["D9"]),
        data = data.frame(D9 = bindori_select_vars["D9"]),
        vars = c("D9"), cont.na = NULL,
        msel = NULL, stat = "percents", breaks = 5,
        nrow = 2, ncol = 2, rel.size.x = 1,
        utility.stats = c("pMSE", "S_pMSE"),
        cols = c("#1A3C5A","#4187BF"),
        plot = TRUE, table = TRUE)
```

### (2). two-way marginals with utility.tables()
In this part, we focus on the utility in the two-way manner/fashion.
However, with a large number of vars to be included, we need to select a few variables to have a first glance of what does the `utility.tables` look like.

```{r}
## filter out a few variables to run the evaluation
## pls make sure that the vars selected are in alignment 
## with the one-way maginal ones
selected_cols <- c("B3", "B1_1", "B7", "D1", "D2", "E2", "E3")
syn_select_vars <- syn_data_filtered[, selected_cols]
bindori_select_vars <- bindori_dataset_filtered[, selected_cols]
```

```{r}
## S3 method for class 'data.frame'
utility.twoway <- utility.tables(object = data.frame(syn_select_vars), 
                                 data = data.frame(bindori_select_vars),
                                 tables = "twoway",
                                 tab.stats = c("pMSE", "S_pMSE"), 
                                 plot.stat = "S_pMSE", plot = TRUE,  
                                 print.tabs = TRUE)
```
Now, we print the results out with the heatmap-like output plot.
```{r}
utility.twoway$utility.plot
```
```{r}
df_twoway_utility<- data.frame(utility.twoway$tabs)
```
```{r}
write_csv_dest <- "F:/Master-Thesis-DifferentialPrivacy/Tab_Utility/Terrence_twoway_utility.csv"
write.csv(df_utility, write_csv_dest, row.names=FALSE)
```


```{r}
twoway_pMSE_list <- c()
twoway_SpMSE_list <- c()
for (i in 1:85) {
  twoway_pMSE_list <- append(pMSE_list, compare_plots[[i]]$tab.utility[1])
  twoway_SpMSE_list <- append(S_pMSE_list, compare_plots[[i]]$tab.utility[2])
  }
```

```{r}
df_utility_twoway <- data.frame()
data.frame(vars_list=colnames(syn_select_vars),
                 pMSE=pMSE_list,
                 S_pMSE=S_pMSE_list)
```