DataSimulation.Rmd

---
title: "Data Simulation"
output: 
  prettydoc::html_pretty:
    theme: architect
    highlight: github
    toc: TRUE
---


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

The purpose of this training is to enable the participants to simulate dummy data which they can use for their analyses or visualizations. 

## Helper functions
We will need some functions from the dplyr package. If `library(dplyr)` gives you an error, you might have to install the package first, by typing into your console: `install.packages("dplyr")`.

```{r,message=FALSE}
set.seed(64) 

library(dplyr)
```


## Simulating one variable

Let's start with the most simple but most time-consuming way. Type everything manually and save it in a vector:

```{r}
client_gen <- c("Millenial","Gen X","Millenial",
                "Baby Boomer","Gen X","Millenial","Gen X")

data.frame(id=1:7,client_gen)
```


### Categorical variables
For categorical variables, we can save some time using the `sample` function. You specify first the possible values and then how many of these values you would like to pick. If you want to allow values to be picked more than once, make sure to set `replace=TRUE`.

```{r}
client_gen <- sample(c("Millenial","Gen X","Baby Boomer"),7,replace=T)

data.frame(id=1:7,client_gen)
```


### Numerical variables

```{r}
client_age <- sample(1:100,size=7,replace=T)

data.frame(id=1:7,client_age)
```


### Distributions

```{r}
client_age <- runif(7,min=1,max=100)

data.frame(id=1:7,client_age)

```

We can use the `round()` function to round each value to their next integer. But uniformly distributed variables are not always what we want. Imagine that we simulate 10000 clients and distributes their ages uniformly. Then there are as many 99 year old clients as there are 50 year old clients:

```{r}
runif(10000,1,100) %>% hist()
```

But we can easily access a whole list of other distribution functions, like the famous Normal distribution (with mean and standard deviation as parameters).

```{r}
rnorm(10000,mean=50,sd=20) %>% hist()
```

If we want to limit the values to not be smaller than 0 or larger than 100, we can use pmin and pmax.

```{r}
rnorm(10000,mean=50,sd=20) %>% pmax(0) %>% pmin(100) %>% hist()
```

For many applications (like balance distribution or any data that contains outliers) I like to use the Exponential distribution (with parameter rate and expectation 1/rate).

```{r}
rexp(10000,rate=0.01) %>% hist()
```

If you want to explore further probability distributions check out this [link](https://www.stat.umn.edu/geyer/old/5101/rlook.html).

### Putting variables together

To create our first simulated dataframe, we can start by simulating the variables separately and then putting them together.

```{r}
set.seed(61)

k <- 7

id <- 1:k
name <- c("Frank","Dorian","Eva","Elena","Andy","Barbara","Yvonne")
age <- rnorm(k,mean=30,sd=10) %>% pmax(18) %>% round()
ocupation <- sample(c("analyst","manager","sr analyst"),k,replace=T,prob=c(10,2,3))
balance <- rexp(k,rate=0.001) %>% round(2)
married <- sample(c("Yes","No"),k,replace=T,prob=c(0.6,0.4))

data <- data.frame(client_id=id,name,age,ocupation,balance,married_flg=married)
data
```

Great! We just simulated a dataset which we can use now for visualization or modelling purposes. 

## Simulating dependent variables
If we want to have a dataset that "makes sense" from a real world perspective, it would be great if managers in general had higher balances than analysts? Or if 18 years old clients are less likely to be married than 30 years olds?
In this section we are going to have a look at techniques to create dependence between variables.

### Rule based

We can use `ifelse()` and `case_when()` to create new variables that depend on others.

```{r}
k <- 7
married <- sample(c("Y","N"),k,replace=T)

data <- data.frame(id=1:k,married)

data %>% mutate(
  age=ifelse(married=="Y",rnorm(k,45,10),rnorm(k,30,10)) %>% pmax(18) %>% round()
)

```

In this small example we will not see the effect, but when we simulate 1000 clients and take a look at their average age, we can see that there is a difference between the two groups.

```{r}
k <- 1000
married <- sample(c("Y","N"),k,replace=T)

data <- data.frame(id=1:k,married)

data %>% mutate(
  age=ifelse(married=="Y",rnorm(k,45,10),rnorm(k,30,10)) %>% 
    pmax(18) %>% 
    round()
    ) %>% 
  group_by(married) %>% 
  summarise(avg_age=mean(age))
```

```{r}
k <- 1000

ocupation <- sample(c("analyst","manager","sr analyst"),k,replace=T,prob=c(10,2,3))

data <- data.frame(id=1:k,ocupation)

data <- data %>% mutate(balance=case_when(
  ocupation=="analyst" ~ rexp(k,0.01),
  ocupation=="sr analyst" ~ rexp(k,0.005),
  TRUE ~ rexp(k,0.001) #this is the else case
))

#Check the average balance per group
data %>% 
  group_by(ocupation) %>% 
  summarise(avg_balance=mean(balance))
```


### Correlation based

If we just deal with numeric variables and want to have a slightly more complex connection between the different variables, we can also try this approach, for which we specify a correlation matrix beforehand and reorder our variables afterwards so that they match the desired correlation.

Of course, we need to find reasonable correlation values, for example between age and number of kids (slightly positively correlated) or between savings and number of kids (slightly negatively correlated).

```{r}
set.seed(64)

k <- 2000

age <- rnorm(k,mean=35,sd=10) %>% pmax(18) %>% round()
balance <- rexp(k,rate=0.001) %>% round(2)
tenure <- rnorm(k,mean=15,sd=5) %>% pmax(1) %>% round()
kids_cnt <- sample(0:5,k,replace=T,prob=c(100,120,80,30,5,1))


data <- data.frame(age,balance,kids_cnt,tenure)
data %>% head(7)
```

We directly see that there are things that don't make sense, like the 22-years-old with a tenure of 22 years.

To improve this, we want to reshuffle the rows and get a distribution close to a desired one. First we simulate a helping dataset of same size, where every entry is random normal distributed.

```{r}
nvars <- ncol(data)
numobs <- nrow(data)

set.seed(3)
rnorm_helper <- matrix(rnorm(nvars*numobs,0,1),nrow=nvars)

```

The correlation of this matrix should be close to the identity matrix.

```{r}
cor(t(rnorm_helper))
```

Next, we specify our desired correlation matrix:

```{r}
Q <- matrix(c(1,0.3,0.4,0.2,  0.3,1,0,0.3,  0.4,0,1,-0.3,  0.2,0.3,-0.3,1),ncol=nvars)

Q
```

We can now multiply the `rnorm_helper` matrix with the Cholesky decomposition of our desired correlation matrix `Q`. Why this works, is explained in the following comment. If you are not interested in mathematical details, you can skip this part.

![](img/Cholesky.png)

(Explanation found [here](https://math.stackexchange.com/q/163472))

```{r}
L <- t(chol(Q))
Z <- L %*% rnorm_helper
```

Good, now we convert this new data to a data frame and give it the name of our original data. The correlation of this dataset is close to our desired outcome.

```{r}
raw <- as.data.frame(t(Z),row.names=NULL,optional=FALSE)
names(raw) <- names(data)
head(raw,7,addrownums=FALSE)
cor(raw)
```

However, this dataset `raw` does not have anything to do with our original data. It is still our transformed random normal data. But as we know that this dataset has the correct correlation, we can use this to reorder the rows of our other dataset.

And then we just replace the largest value of the random normal dataset with the largest value in our dataset, the second largest with the second largest etc. We go column by column and repeat this procedure.

```{r}

for(name in names(raw)) {
  raw <- raw[order(raw[,name]),]
  data <- data[order(data[,name]),]
  raw[,name] <- data[,name]
}

```

Let's check the correlation of this new dataset. It is not exactly what we wished, but close enough. The reason for this is that our variables take less values than a random normal distributed variable (e.g. kids count just takes values between 0 and 5).

```{r}
cor(raw)
```


We can also take a look at the reshuffled dataset.

```{r}
head(raw,7,addrownums=FALSE)
```



## Exporting 

To export this dataset, we can use the base R function.

```{r,eval=F}
write.csv(raw,"data.csv")
```

Alternatively, if this is not fast enough, we can also use the `fwrite` function from the data.table package which is much faster.

```{r,eval=F}
library(data.table)
fwrite(raw,"data.csv")
```