---
title: "Simulating Data Stub File"
author: "Lisa DeBruine"
---

```{r lib, message=FALSE}
library(tidyverse)
# devtools::install_github("debruine/faux", build_vignettes = TRUE)
library(faux) 
library(afex) # for anova and lmer
library(broom.mixed) # to make tidy tables of lmer output
library(kableExtra) # to make nicer graphs
set.seed(8675309) # this makes sure your script uses the same set of random numbers each time you run the full script (never set this inside a function or loop)
```

## Independent samples

Let's start with a simple independent-samples design where the variables are from a normal distribution. Each subject produces one score (in condition A or B). What we need to know about these scores is:

* How many subjects are in each condition?
* What are the score means?
* What are the score variances (or SDs)?

### Parameters

I start simulation scripts by setting parameters for these values.

```{r ind-vars}

A_sub_n <- 50
B_sub_n <- 50
A_mean  <- 10
B_mean  <- 11
A_sd    <- 2.5
B_sd    <- 2.5

```

### Scores

We can the generate the scores using the `rnorm()` function.

```{r ind-dat}
A_scores <- rnorm(A_sub_n, A_mean, A_sd)
B_scores <- rnorm(B_sub_n, B_mean, B_sd)
```

You can stop here and just analyse your simulated data with `t.test(A_scores, B_scores)`, but usualy you want to get your simulated data into a data table that looks like what you might eventually import from a CSV file with your actual experimental data.

```{r ind-tibble}
dat <- tibble(
  sub_condition = rep( c("A", "B"), c(A_sub_n, B_sub_n) ),
  score = c(A_scores, B_scores)
)
```

### Check your data

Always examine your simulated data after you generate it to make sure it looks like you want.

```{r ind-check}
summary_dat <- dat %>%
  group_by(sub_condition) %>%
  summarise(n = n() ,
            mean = mean(score),
            sd = sd(score))

summary_dat
```


### Analysis

```{r ind-test}
t.test(score~sub_condition, dat)
```

### Function

You can wrap all this in a function so you can run it many times to do a power calculation. Put all your parameters as arguments to the function.

```{r ind-func}

ind_sim <- function(A_sub_n, B_sub_n, 
                    A_mean, B_mean, 
                    A_sd, B_sd) {
  A_scores <- rnorm(A_sub_n, A_mean, A_sd)
  B_scores <- rnorm(B_sub_n, B_mean, B_sd)
  
  dat <- tibble(
    sub_condition = rep( c("A", "B"), c(A_sub_n, B_sub_n) ),
    score = c(A_scores, B_scores)
  )
  t <- t.test(score~sub_condition, dat)
  
  # return just the values you care about
  list(
    t = t$statistic[[1]],
    ci_lower = t$conf.int[[1]],
    ci_upper = t$conf.int[[2]],
    p = t$p.value[[1]],
    estimate = t$estimate[[1]] - t$estimate[[2]]
  )
}

```

Now run your new function with the values you used above.

```{r}
# str() prints the resulting list in a shorter format
ind_sim(50, 50, 10, 11, 2.5, 2.5) %>% str()
```

Now you can use this function to run many simulations. 

```{r}
mysim <- map_df(1:1000, ~ind_sim(50, 50, 10, 11, 2.5, 2.5))
```

Now you can graph the data from your simulations.

```{r paired-sim-fig, fig.cap = "Distribution of results from simulated independent samples data"}
mysim %>%
  gather(stat, value, t:estimate) %>%
  ggplot() + 
  geom_density(aes(value, color = stat), show.legend = FALSE) +
  facet_wrap(~stat, scales = "free")
```

You can calculate power as the proportion of simulations on which the p-value was less than your alpha.

```{r}
alpha <- 0.05
power <- mean(mysim$p < alpha)
power
```



## Paired samples

Now let's try a paired-samples design where the variables are from a normal distribution. Each subject produces two scores (in conditions A and B). What we need to know about these two scores is:

* How many subjects?
* What are the score means?
* What are the score variances (or SDs)?
* What is the correlation between the scores?

### Parameters {#paired-params}

```{r paired-vars}

sub_n <- 100
A_mean <- 10
B_mean <- 11
A_sd <- 2.5
B_sd <- 2.5
AB_r <- 0.5

```


### Correlated Scores

You can then use `rnorm_multi()` to generate a data table with simulated values for correlated scores:

```{r sim-design}
dat <- faux::rnorm_multi(
    n = sub_n, 
    vars = 2, 
    r = AB_r, 
    mu = c(A_mean, B_mean), 
    sd = c(A_sd, B_sd), 
    varnames = c("A", "B")
  )
```

You can also do this using the `MASS::mvrnorm` function, but `faux::rnorm_multi` is easier when you have more variables to simulate.

```{r}
# make the correlation matrix
cormat <- matrix(c(   1, AB_r,
                   AB_r,    1), 
             nrow = 2, byrow = TRUE)

# make a corresponding matrix of the variance 
# (multiply the SDs for each cell)
varmat <- matrix(c(A_sd * A_sd, A_sd * B_sd,
                   A_sd * B_sd, B_sd * B_sd), 
             nrow = 2, byrow = TRUE) 

# create correlated variables with the specified parameters
S <- MASS::mvrnorm(n = sub_n, 
                   mu = c(A_mean, B_mean), 
                   Sigma = cormat * varmat)
dat <- data.frame(
  A = S[, 1],
  B = S[, 2]
)

```


### Check your data

Now check your data; `faux` has a function `get_params()` that gives you the correlation table, means, and SDs for each numeric column in a data table.

```{r paired-check}
faux::get_params(dat)
```

### Analysis

```{r paired-test}
t.test(dat$A, dat$B, paired = TRUE)
```

### Function

```{r paired-func}

paired_sim <- function(sub_n, A_mean, B_mean, A_sd, B_sd, AB_r) {

  dat <- faux::rnorm_multi(
    n = sub_n, 
    vars = 2, 
    r = AB_r, 
    mu = c(A_mean, B_mean), 
    sd = c(A_sd, B_sd), 
    varnames = c("A", "B")
  )
  t <- t.test(dat$A, dat$B, paired = TRUE)
  
  # return just the values you care about
  list(
    t = t$statistic[[1]],
    ci_lower = t$conf.int[[1]],
    ci_upper = t$conf.int[[2]],
    p = t$p.value[[1]],
    estimate = t$estimate[[1]]
  )
}

```

Run 1000 simulations and graph the results.

```{r}
mysim_p <- map_df(1:1000, ~paired_sim(100, 10, 11, 2.5, 2.5, .5))
```

```{r ind-sim-fig, fig.cap = "Distribution of results from simulated paired samples data"}
mysim_p %>%
  gather(stat, value, t:estimate) %>%
  ggplot() + 
  geom_density(aes(value, color = stat), show.legend = FALSE) +
  facet_wrap(~stat, scales = "free")
```

```{r}
alpha <- 0.05
power <- mean(mysim_p$p < alpha)
power
```


## Intercept model

Now I'm going to show you a different way to simulate the same design. This might seem excessively complicated, but you will need this pattern when you start simulating data for mixed effects models.

### Parameters

Remember, we used the following parameters to set up our simulation above:

```{r paired-vars2}
sub_n  <- 100
A_mean <- 10
B_mean <- 11
A_sd   <- 2.5
B_sd   <- 2.5
AB_r   <- 0.5
```

From these, we can calculate the grand intercept (the overall mean regardless of condition), and the effect of condition (the mean of B minus A).

```{r}
grand_i   <- (A_mean + B_mean)/2
AB_effect <- B_mean - A_mean
```

We also need to think about variance a little differently. First, calculate the pooled variance as the mean of the variances for A and B (remember, variance is SD squared).

```{r}
pooled_var <- (A_sd^2 + B_sd^2)/2
```

The variance of the subject intercepts is `r` times this pooled variance and the error variance is what is left over. We take the square root (`sqrt()`) to set the subject intercept and error SDs for simulation later.

```{r}
sub_sd   <- sqrt(pooled_var * AB_r)
error_sd <- sqrt(pooled_var * (1-AB_r))
```


### Subject intercepts

Now we use these variables to create a data table for our subjects. Each subject gets an ID and a **random intercept** (`sub_i`). The intercept is simulated from a random normal distribution with a mean of 0 and an SD of `sub_sd`. This represents how much higher or lower than the average score each subject tends to be (regardless of condition).

```{r}
sub <- tibble(
  sub_id = 1:sub_n,
  sub_i = rnorm(sub_n, 0, sub_sd)
)
```

### Observations

Next, set up a table where each row represents one observation. We'll use one of my favourite functions for simulation: `crossing()`. This creates every possible combination of the listed factors (it works the same as `expand.grid()`, but the results are in a more intuitive order). Here, we're using it to create a row for each subject in each condition, since this is a fully within-subjects design.

```{r}
obs <- crossing(
  sub_id = 1:sub_n,
  condition = c("A", "B")
)
```

### Calculate the score

Next, we join the subject table so each row has the information about the subject's random intercept and then calculate the score. I've done it in a few steps below for clarity. The score is just the sum of:

* the overall mean (`grand_i`)
* the subject-specific intercept (`sub_i`)
* the effect (`effect`): the numeric code for condition (`condition.e`) multiplied by the effect of condition (`AB_effect`)
* the error term (simulated from a normal distribution with mean of 0 and SD of `error_sd`)

```{r}
dat <- obs %>%
  left_join(sub, by = "sub_id") %>%
  mutate(
    condition.e = recode(condition, "A" = -0.5, "B" = 0.5),
    effect = AB_effect * condition.e,
    error = rnorm(nrow(.), 0, error_sd),
    score = grand_i + sub_i + effect + error
  )
```

You can use the following code to put the data table into a more familiar "wide" format.

```{r}
dat_wide <- dat %>%
  select(sub_id, condition, score) %>%
  spread(condition, score)
```


Then you can use the `get_params` function to check this looks correct (remove the subject ID to leave it out of the table, since it's numeric).

```{r}
dat_wide %>%
  select(-sub_id) %>%
  faux::get_params()
```

### Analyses

You can analyse the data with a paired-samples t-test from the wide format:

```{r}
t.test(dat_wide$A, dat_wide$B, paired = TRUE)
```

Or in the long format:

```{r}
t.test(score ~ condition, dat, paired = TRUE)
```

You can analyse the data with ANOVA using the `aov_4()` function from `afex`. (Notice how the F-value is the square of the t-value above.)

```{r}
afex::aov_4(score ~ (condition.e | sub_id), data = dat)
```


You can even analyse the data with a mixed effects model using the `lmer` function (the `afex` version gives you p-values, but the `lme4` version does not). 

```{r}
lmem <- afex::lmer(score ~ condition.e + (1 | sub_id), data = dat)

broom.mixed::tidy(lmem, effects = "fixed")%>% 
  # nicer formatting with p-values
  mutate_if(~max(.) < .001, ~as.character(signif(., 3)))
```

## Functions

We can put everything together into a function where you specify the subject number, means, SDs and correlation, it translates this into the intercept specification, and returns a data table.

```{r sim-paired-dat-func}

sim_paired_data <- function(sub_n = 100, 
                            A_mean, B_mean, 
                            A_sd, B_sd, AB_r) {
  
  grand_i <- (A_mean + B_mean)/2
  AB_effect <- B_mean - A_mean
  pooled_var <- (A_sd^2 + B_sd^2)/2
  sub_sd <- sqrt(pooled_var * AB_r)
  error_sd <- sqrt(pooled_var * (1-AB_r))
  
  sub <- tibble(
    sub_id = 1:sub_n,
    sub_i = rnorm(sub_n, 0, sub_sd)
  )
  
  crossing(
    sub_id = 1:sub_n,
    condition = c("A", "B")
  ) %>%
  left_join(sub, by = "sub_id") %>%
  mutate(
    condition.e = recode(condition,  "A" = -0.5, "B" = +0.5),
    effect = AB_effect * condition.e,
    error = rnorm(nrow(.), 0, error_sd),
    score = grand_i + sub_i + effect + error
  ) %>%
  select(sub_id, condition, score) %>%
  spread(condition, score)
}
```

```{r}
sim_paired_data(100, 10, 12, 2.5, 2.5, .5) %>%
  select(-sub_id) %>%
  faux::get_params()
```

Set the `sub_n` to a very high number to see that the means, SDs and correlations are what you specified.

```{r}
sim_paired_data(10000, 0, 0.5, 1, 1, .25) %>%
  select(-sub_id) %>%
  faux::get_params()
```

It might be more useful to create functions to translate back and forth from the distribution specification to the intercept specification.

### Distribution to intercept specification

```{r dist2int-func}
dist2int <- function(mu = 0, sd = 1, r = 0) {
  A_mean <- mu[1]
  # set B_mean to A_mean if mu has length 1
  B_mean <- ifelse(is.na(mu[2]), mu[1], mu[2])
  A_sd <- sd[1]
  # set B_sd to A_sd if sd has length 1
  B_sd <- ifelse(is.na(sd[2]), sd[1], sd[2])
  AB_r <- r
  pooled_var <- (A_sd^2 + B_sd^2)/2
  
  list(
    grand_i = (A_mean + B_mean)/2,
    AB_effect = B_mean - A_mean,
    sub_sd = sqrt(pooled_var * AB_r),
    error_sd = sqrt(pooled_var * (1-AB_r))
  )
}
```


```{r}
dist2int()
```



```{r}
dist2int(mu = c(100, 105), sd = c(10.5, 9.5), r = 0.5)
```


### Intercept to distribution specification

```{r int2dist-func}
int2dist <- function(grand_i = 0, 
                     AB_effect = 0, 
                     sub_sd = 0, 
                     error_sd = 1) {
  pooled_var <- sub_sd^2 + error_sd^2
   
  list(
    A_mean = grand_i - 0.5 * AB_effect,
    B_mean = grand_i + 0.5 * AB_effect,
    A_sd = sqrt(pooled_var),
    B_sd = sqrt(pooled_var),
    AB_r = sub_sd^2 / pooled_var
  )
}
```

```{r}
int2dist()
```


```{r}
int2dist(102.5, 5, 7.08, 7.08)
```

