Measurement Simulations: Test-Retest Reliability

Setup

Learning goals

  1. Simulate data for two experiments and compute test-retest reliability
  2. Practice some tidyverse (pivot_longer, mutate, select, and add onto existing base ggplot skills (geom_point, geom_jitter, facet_wrap, geom_line)
  3. Run a basic correlation (cor.test and interpret differences in observed reliability based on differences in the simulated data)

Import the libraries we need

library(tidyverse)
Warning: package 'ggplot2' was built under R version 4.2.3
Warning: package 'tidyr' was built under R version 4.2.3
Warning: package 'readr' was built under R version 4.2.3
Warning: package 'dplyr' was built under R version 4.2.3
Warning: package 'stringr' was built under R version 4.2.3
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2) # plotting

Define the simulation function - same as before

This makes “tea data”, a tibble (dataframe) where there are a certain number of people in each condition (default = 48, i.e., n_total, with n_total/2 in each half)

The averages of the two conditions are separated by a known effect (“delta”) with some variance (“sigma”). You can change these around since we’re simulating data!

set.seed(123) # good practice to set a random seed, or else different runs get you different results
make_tea_data <- function(n_total = 48,
                          sigma = 1.25,
                          delta = 1.5) {
  n_half <- n_total / 2
  tibble(condition = c(rep("milk first", n_half), rep("tea first", n_half)),
         rating = c(round(
           rnorm(n_half, mean = 3.5 + delta, sd = sigma)
         ),
         round(rnorm(
           n_half, mean = 3.5, sd = sigma
         )))) |>
    mutate(rating = if_else(rating > 10, 10, rating),
           # truncate if greater than max/min of rating scale
           rating = if_else(rating < 1, 1, rating))
}

1. Make a dataframe with our simulated data

  1. Input more participants (50 per condition) with a bigger average difference between conditions (2 points), with variance between participants at 2 points (sigma)
tea_data_largern <- make_tea_data(n_total = 50, sigma=2, delta=2)

2. Creating the second experiment

Now create a new column in your tibble for the second experiment.

Goals: practicing tidyverse and simulating data

Here, the rating of the simulated second experiment data is each participants first rating with some variance (people are likely to not say exactly the same thing, but their scores are likely to be similar)

TIPS:

Recommend running rowwise() in your pipe before creating the new condition to force tidyverse to sample a new random value for each row

Make your next dataframe A NEW NAME so that you’re not rewriting old dataframes with new ones and getting confused

Hint: you can use to_sample = -3:3 with the sample function and specifies the possible values you want to sample from

# YOUR CODE HERE
new_tea_data <- tea_data_largern |> rowwise() |> 
  rowwise() |>
  mutate(rate_exp1=rating) |>
  mutate(diff=sample(-3:3, 1)) |>
  mutate(rate_exp2=rate_exp1+diff) |>
  mutate(
    rate_exp2=if_else(rate_exp2>10, 10, rate_exp2),
    rate_exp2=if_else(rate_exp2<1, 1, rate_exp2))

3. Make a plot and compute the correlation

Examine how the ratings are correlated across these simulations, by making a plot

# install.packages("ggpubr")
# library(ggpubr)

ggplot(new_tea_data, aes(x = rate_exp1, y = rate_exp2)) +
  geom_point(aes(col=condition)) +
  labs(title="Scatter Plot", x="rating from experiment 1", y="rating from experiment 2") +
  geom_jitter(alpha=0.5) +
  geom_smooth(method='lm') +
  # ggpubr::stat_cor() + # what is this
  theme_minimal()
`geom_smooth()` using formula = 'y ~ x'

Hint: 1. Try out geom_point which shows you the exact values 2. Then try out geom_jitter which shows you the same data with some jitter around height / width

Bonus: 3. Use alpha to make your dots transparent 4. Use ylab and xlab to make nice axes labels 5. Use geom_smooth() to look at the trend_line 6. Try making each dot different by condition

Now examine – how correlated are your responses? What is your test-retest reliability?

cor.test(new_tea_data$rate_exp1, new_tea_data$rate_exp2)

    Pearson's product-moment correlation

data:  new_tea_data$rate_exp1 and new_tea_data$rate_exp2
t = 6.1302, df = 48, p-value = 1.585e-07
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.4712324 0.7944689
sample estimates:
      cor 
0.6626612

Hint: You’ll need a correlation function - cor.test Bonus: You can also find a ggplot function to overlay the correlation on top of your plot

4. Make another plot – like the one from the chapter, where each line should connect an individual subject

  1. First, Use pivot_longer to make the dataframe longer
  2. Use facet_wrap(~condition) to make two plots, one for milk_first and one for tea_first
  3. Use geom_line – with grouping specification by sub_id – to connect individual datapoints for each condition

Hint aes(group = sub_id))

# Make new, longer dataframe 
longer_df <- new_tea_data |> 
  ungroup() |>
  mutate(sub_id=row_number()) |>
  select(-c('rating', 'diff')) |>
  pivot_longer(
    cols = c(starts_with('rate_exp')),
    names_to = 'exp',
    values_to = 'rate_exp'
  )
ggplot(longer_df, aes(x = exp, y = rate_exp)) +
  geom_point(aes(col=condition)) +
  labs(title="Scatter Plot", x="experiment id", y="rating") +
  geom_jitter(alpha=0.5) +
  facet_wrap(~condition) +
  geom_line(aes(group=sub_id)) +
  # ggpubr::stat_cor() + # what is this
  theme_minimal()

5. OK, now go back and change things and test your intuition about how this works.

How does reliability change if you increase the variance between participants (sigma) in the first experiment simulated data?

How does reliability change if you change how much variation you allow between the first and second experiment?

How does reliability change if you increase sample size, holding those things constant?

Hint: copy the code from above where you made your new dataframe with experiment number 2, copy the correlation computation code, and just run this block over an over, modifying the code (command-shift-enter on Mac runs the block)

# YOUR CODE HERE
sigma_val <- 1 # change the sigma value, originally 2
sample_size <- 20 # change the sample size, originally 5
individual_diff <- 2 # change how much individual vary across experiments

tea_data_different_settings <- make_tea_data(n_total = sample_size, sigma=sigma_val, delta=2)

new_tea_data_different_settings <- tea_data_different_settings |> rowwise() |> 
  rowwise() |>
  mutate(rate_exp1=rating) |>
  mutate(diff=sample(-individual_diff:individual_diff, 1)) |>
  mutate(rate_exp2=rate_exp1+diff) |>
  mutate(
    rate_exp2=if_else(rate_exp2>10, 10, rate_exp2),
    rate_exp2=if_else(rate_exp2<1, 1, rate_exp2))

cor.test(new_tea_data_different_settings$rate_exp1, new_tea_data_different_settings$rate_exp2)

    Pearson's product-moment correlation

data:  new_tea_data_different_settings$rate_exp1 and new_tea_data_different_settings$rate_exp2
t = 5.9622, df = 18, p-value = 1.218e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5821419 0.9241029
sample estimates:
      cor 
0.8147699

Ans:

  1. Larger sigma, greater correlation
  2. Larger sampler size, greater correlation
  3. Larger individual difference, smaller correlation