In debruine/data-sim-workshops: Data Simulation Workshops

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
• Re-analysis code
• Original paper

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.

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.

s2_analyse <- function(data) {
}

s2_analyse(data)

Now, simulate the data 500 times.

s2 <- NULL

Run the analysis on each data set.

s2_sim <- NULL

head(s2_sim)

Summarise the p.value column to get power.

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.

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.

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.

s1 <- NULL

Run the analysis on each data set.

s1_sim <- NULL

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

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.

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.

s3_analyse <- function(data) {

}

Run the analysis on each data set.

s3_sim <- NULL

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()!

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

ss1 <- NULL

ss1_sim <- NULL

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()!

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

ss2_sim <- NULL

ss2_power <- NULL

Study S3

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

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

ss3_sim <- NULL

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.

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.

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.

debruine/data-sim-workshops documentation built on April 27, 2024, 8:13 a.m.