---
title: "Calorie Placement Re-Simulation"
output: 
  html_document:
    df_print: kable
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(error = TRUE)
library(tidyverse)
library(faux)
library(afex)
library(emmeans)
faux_options(plot = FALSE)

set.seed(8675309)
```

## Data Source

We will be replicating some of the re-analyses in Francis & Thunell's (2020) Meta-Psychology paper: Excess success in "Don't count calorie labeling out: Calorie counts on the left side of menu items lead to lower calorie food choices".

They ran power analyses for all 6 studies in Dallas, Liu, and Ubel's (2019) study showing that people order food with significantly fewer calories when the calorie count was placed to the left of the item than to the right (or having no calorie label). They then used these power estimates to calculate the probability of all 6 out of 6 studies being significant, given the observed power of each study.

* [Re-analysis](https://doi.org/10.15626/MP.2019.2266)
* [Re-analysis code](https://osf.io/xrdhj/)
* [Original paper](https://doi.org/10.1002/jcpy.1053)

Table 1 of the re-analysis paper provides all of the parameters we will need. 

## Reanalyses

### Study 2

We'll start with S2 because the analysis is very straightforward. It's a between-subjects design, where 143 subjects saw calorie placement on the left and their mean calories ordered were 1249.83 (SD = 449.07), while 132 subjects saw calorie placement on the right and their mean calories ordered were 1362.31 (SD = 447.35).

Let's first simulate a single data table with these parameters and set up our analysis.

```{r}
data <- NULL
```

Wrap the analysis in a function using the `tidy()` function from {broom} to get the results in a tidy table. Check that it works by running it on the single data set above.

```{r}
s2_analyse <- function(data) {
}

s2_analyse(data)
```


Now, simulate the data 500 times.

```{r}
s2 <- NULL
```

Run the analysis on each data set.

```{r}
s2_sim <- NULL

head(s2_sim)
```

Summarise the `p.value` column to get power.

```{r}
s2_power <- NULL
```

Compare this value (`r s2_power`) with the value in the paper (0.5426).

### Study 1

Study 1 is a little more complicated because the design includes a "no label" condition, so the decision rule for supporting the hypothesis is more complicated.

The data simulation is relatively straightforward, though.

```{r}
mu = c(left = 654.53, right = 865.41, none = 914.34)
sd = c(left = 390.45, right = 517.26, none = 560.94)
n =  c(left =  45,    right =  54,    none =  50)

data <- NULL
```

Set up the analysis. Here, we really just care about three p-values, so we'll just return those. We can use a function from the {emmeans} package to check the two pairwise comparisons.

```{r}
afex::set_sum_contrasts() # avoids annoying afex message on each run
afex_options(include_aov = TRUE) # we need aov for emmeans

s1_analyse <- function(data) {
  # main effect of placement
  a <- afex::aov_ez(
    id = "id",
    dv = "calories",
    between = "placement",
    data = data
  )
  
  # contrasts
  e <- emmeans(a, "placement")
  c1 <- list(lr = c(-0.5, 0.5, 0),
             ln = c(-0.5, 0, 0.5))
  b <- contrast(e, c1, adjust = "holm") |>
    broom::tidy()
  
  data.frame(
    p_all = a$anova_table$`Pr(>F)`[[1]],
    p_1 = b$adj.p.value[[1]],
    p_2 = b$adj.p.value[[2]]
  )
}

s1_analyse(data)
```

Let's just replicate this 100 times so the simulation doesn't take too long to run at first. We can always increase it later after we've run some sense checks.

```{r}
s1 <- NULL
```

Run the analysis on each data set.

```{r}
s1_sim <- NULL
```

Calculating power is a little trickier here, as all three p-values need to be significant here to support the hypothesis.

```{r}
s1_power <- NULL
```

Compare this value (`r s1_power`) with the value in the paper (0.4582).

### Study 3

Now you can use the pattern from Study 1 to analyse the data for Study 3. We'll start with the repeated data set.

```{r}
mu = c(left = 1428.24, right = 1308.66, none = 1436.79)
sd = c(left =  377.02, right =  420.14, none =  378.47)
n =  c(left =   85,    right =   86,    none =   81)

s3 <- NULL
```

These data were collected in the Hebrew language, which reads right to left, so the paired contrasts will be different.

```{r}
s3_analyse <- function(data) {

}
```

Run the analysis on each data set.

```{r}
s3_sim <- NULL
```

```{r}
s3_power <- NULL
```

Compare this value (`r s3_power`) with the value in the paper (0.3626).


### Study S1

Now you can use the pattern from Study 2 to analyse the data for Study S1. You can even reuse the analysis function `s2_analyse()`!

```{r}
mu = c(left = 185.94, right = 215.73)
sd = c(left =  93.92, right =  95.33)
n =  c(left =  99,    right =  77)

ss1 <- NULL
```

```{r}
ss1_sim <- NULL
```


```{r}
ss1_power <- NULL
```


### Study S2

Now you can use the pattern from Study 1 to analyse the data for Study S2. You can even reuse the analysis function `s1_analyse()`!

```{r}
mu = c(left = 1182.15, right = 1302.23, none = 1373.74)
sd = c(left =  477.60, right =  434.41, none =  475.77)
n =  c(left =  139,    right =  141,    none = 151)

ss2 <- NULL
```

```{r}
ss2_sim <- NULL
```

```{r}
ss2_power <- NULL
```

### Study S3

Now you can use the pattern from Study 1 to analyse the data for Study S3.

```{r}
mu = c(left = 1302.03, right = 1373.15, none = 1404.35)
sd = c(left =  480.02, right =  442.49, none =  422.03)
n  = c(left =  336,    right =  337,    none =  333)

ss3 <- NULL
```

```{r}
ss3_sim <- NULL
```

```{r}
ss3_power <- NULL
```

## Conclusion

Now that you've calculated power for each of the 6 studies, just multiply the 6 power values together to get the probability that all 6 studies will be significant.


```{r}
power_table <- tribble(
  ~study, ~power_ft, ~ power_my,
  "1",       0.4582,   s1_power,
  "2",       0.5426,   s2_power,
  "3",       0.3626,   s3_power,
  "S1",      0.5358,  ss1_power,
  "S2",      0.5667,  ss2_power,
  "S3",      0.4953,  ss3_power
)

power_table
```

The `reduce()` function from {purrr} applies a function sequentially over a vector, so can give up the product of all the values in the power columns.

```{r}
prob_ft <- purrr::reduce(power_table$power_ft, `*`)
prob_my <- purrr::reduce(power_table$power_my, `*`)
```

The Francis & Thunell paper showed a `r prob_ft` probability of getting 6 of 6 studies significant. Our re-simulation showed a `r prob_my` probability.