← Back to Code Reviews
Ongchoco et al., 2023
Code Review · AI-assisted with ReproAI, complementary to the manuscript review.
# Effect Size Summaries from Ongchoco et al., 2023 ------------------------
# Paper: https://doi.org/10.1073/pnas.2303883120
# Dataset: https://www.pnas.org/doi/suppl/10.1073/pnas.2303883120/suppl_file/pnas.2303883120.sd01.xlsx
library(openxlsx) # to open the original dataset
library(MOTE) # for etasq
library(esc) # for cohen's d
library(ggplot2) # for plotting
library(metafor) # package needed for TES
# Open dataset
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E1 Data (PostExclusions)") # dataset available from the PNAS website
### Plot Eta-Squared with 90% CIs for interaction effects; https://doi.org/10.1073/pnas.2303883120 Studies 1-6
library(ggplot2)
library(MOTE)
etas <- data.frame("eta" = as.numeric()
, "etalow" = as.numeric()
, "etahigh" = as.numeric()
)
etas <- rbind(etas, unlist(MOTE::eta.F(dfm = 1, dfe = 929, Fvalue = 6.57, a = .1)[1:3])) # 1
etas <- rbind(etas, unlist(MOTE::eta.F(dfm = 1, dfe = 796, Fvalue = 7.01, a = .1)[1:3]))# 2 # reported eta = 0.094
etas <- rbind(etas, unlist(MOTE::eta.F(1, 796, 5.17)[1:3]))
etas <- rbind(etas, unlist(MOTE::eta.F(1, 796, 3.98)[1:3]))
etas <- rbind(etas, unlist(MOTE::eta.F(1, 796, 7.17)[1:3])) # reported eta = 0.005
etas <- rbind(etas, unlist(MOTE::eta.F(1, 996, 7.01)[1:3]))
names(etas) <- c("eta", "lower", "upper")
etas$study <- as.character(1:6)
ggplot(etas, aes(x = eta, y = study, xmin = lower, xmax = upper)) + geom_point() + geom_errorbar() + xlim(c(-.01, .05)) + geom_abline(slope = 9999999, intercept = 0, lty = 2, col = "dark red") +
theme_bw() + ylab("Study") + xlab("Effect Size (etasq)") + ggtitle("Effect Sizes and CIs for Ongchoco et al., 2023") + scale_y_discrete(limits=rev)
### Plot Cohen's d with 95% CIs for https://doi.org/10.1073/pnas.2303883120 Studies 1-6 (incidental anchoring effect in the no doorway condition)
library(ggplot2)
library(esc)
ds <- data.frame("d" = as.numeric()
, "dlow" = as.numeric()
, "dhigh" = as.numeric()
, "se" = as.numeric()
)
ds <- rbind(ds, unlist(esc::esc_mean_sd(grp1m = 103.46, grp1sd = 57.88, grp1n = 200, grp2m = 90.11, grp2sd = 53.79, grp2n = 200)[c(1, 4, 5, 2)])) #s1
ds <- rbind(ds, unlist(esc::esc_mean_sd(grp1m = 90.30, grp1sd = 44.18, grp1n = 200, grp2m = 81.39, grp2sd = 42.10, grp2n = 200) [c(1, 4, 5, 2)])) #s2
ds <- rbind(ds, unlist(esc::esc_mean_sd(grp1m = 76.11, grp1sd = 39.63, grp1n = 200, grp2m = 66.34, grp2sd = 40.60, grp2n = 200) [c(1, 4, 5, 2)])) #s
ds <- rbind(ds, unlist(esc::esc_mean_sd(grp1m = 12.62, grp1sd = 9.98, grp1n = 200, grp2m = 10.71, grp2sd = 8.88, grp2n = 200) [c(1, 4, 5, 2)]))
ds <- rbind(ds, unlist(esc::esc_mean_sd(grp1m = 12.96, grp1sd = 10.96, grp1n = 200, grp2m = 11.46, grp2sd = 10.02, grp2n = 200) [c(1, 4, 5, 2)]))
ds <- rbind(ds, unlist(esc::esc_mean_sd(grp1m = 79.37, grp1sd = 39.16, grp1n = 250, grp2m = 47.48, grp2sd = 28.31, grp2n = 250) [c(1, 4, 5, 2)]))
names(ds) <- c("d", "lower", "upper", "se")
ds$study <- c("1: Item Valuations"
, "2: Item Valuations (Controlling for Visual Complexity)"
, "3: Factual Questions"
, "4: Legal Judgments"
, "5: Legal Judgments (Larger Replication)"
, "6: Explicit Anchoring in Item Valuations")
cohens_ds <- ds
ggplot(ds, aes(x = d, y = study, xmin = lower, xmax = upper)) + geom_point() + geom_errorbar() +
# xlim(c(-.01, .05)) +
geom_abline(slope = 9999999, intercept = 0, lty = 2, col = "dark red") +
theme_bw() + ylab("Study") +
xlab("Effect Size (Cohen's d)") +
ggtitle("Effect Sizes and CIs for Ongchoco et al., 2023\nMain Effects for Anchoring in the No-Door Condition")+
scale_y_discrete(limits=rev)
### p-curve data for interaction effects
# F(1, 929) = 6.57, p = .01053
# F(1, 796) = 7.01, p = .00827
# F(1, 796) = 5.17, p = .02325
# F(1, 796) = 3.98, p = .04638
# F(1, 796) = 7.27, p = .00716
# F(1, 996) = 7.01, p = .00823
# Flexibility in data analysis --------------------------------------------
### load all datasets
## s1
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E1 Data (PreExclusions)")
ds_sub1 <- ds[ds$room_condition == "hall", ]
ds_sub1$study <- 1
head(ds_sub1)
## s2
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E2 Data (PreExclusions)")
ds_sub2 <- ds[ds$room_condition == "hall", ]
ds_sub2$study <- 2
## s3
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E3 Data (PreExclusions)")
ds_sub3 <- ds[ds$room_condition == "hall", ]
ds_sub3$study <- 3
## s4
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E4 Data (PreExclusions)")
ds_sub4 <- ds[ds$room_condition == "hall", ]
ds_sub4$study <- 4
## s5
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E5 Data (PreExclusions)")
ds_sub5 <- ds[ds$room_condition == "hall", ]
ds_sub5$study <- 5
## s6
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E6 Data (PreExclusions)")
ds_sub6 <- ds[ds$room_condition == "hall", ]
ds_sub6$study <- 6
## s7
ds <- openxlsx::read.xlsx("pnas.2303883120.sd01.xlsx", sheet = "E7 Data (PreExclusions)")
ds_sub7 <- ds[ds$expt_condition == "hall", ]
ds_sub7$study <- 7
ds_sub7 <- ds_sub7[, c("anchor", "response_value", "rt", "study")]
names(ds_sub7) <- c("anchor", "response", "rt", "study")
### combine all datasets into one big dataset
ds_hall <- rbind( ds_sub1[, c("anchor", "response", "rt", "study")]
, ds_sub2[, c("anchor", "response", "rt", "study")]
, ds_sub3[, c("anchor", "response", "rt", "study")]
, ds_sub4[, c("anchor", "response", "rt", "study")]
, ds_sub5[, c("anchor", "response", "rt", "study")]
, ds_sub6[, c("anchor", "response", "rt", "study")]
, ds_sub7[, c("anchor", "response", "rt", "study")]
)
# show anchors for all studies
table(ds_hall$anchor, ds_hall$study) # anchors differed between studies
### compute cohen's ds for all studies given different exclusion criteria (post data collection)
## description from the dataset ("Summary" sheet)
# 1. Compute the response time threshold of 2 standard deviations from the grand mean response time
# 2. Compute the response value threshold of x standard deviations from the grand mean response (where x varies per experiment)
# E1: 1.5SD, E2: 1.5SD, E3: 2SD, E4: 2SD, E5: 1.5SD, E6: 1.5SD, E7: 2SD
# (with these values based on pilot ranges of values depending on the prompt and platform [mTurk vs. Prolific], and pre-registered)
# 3. Remove subjects who go beyond the response time threshold computed from #1.
# 4. Remove subjects who go beyond the response value threshold computed from #2.
# 5. E1-E3: Remove item valuations below 10 and above 350.
# 6. E4: Remove reported number of days in the year above 365.
# 7. Take the first 200 (or 400 for E6; or 250 for E7) subjects for each of the four door-high, door-low, hall-high, hall-low conditions.
# define conditions and all of their combinations to be tried out
response_values <- c(NA, 1.5, 2)
response_times <- c(NA, 2)
additional <- c("no", "yes")
# prepare empty dataset
cors <- data.frame( "r" = as.numeric()
, "rlo" = as.numeric()
, "rhi" = as.numeric()
, "exc_value" = as.numeric()
, "exc_rt" = as.numeric()
, "study" = as.numeric()
# , "exc_additional" = as.numeric()
)
for (i in response_values) {
for (j in response_times) {
for (l in unique(ds_hall$study)) {
# filter by study
ds_temp <- ds_hall[ds_hall$study == l, ]
# exclude response value outliers
if (!is.na(i)) {
ds_temp[, "response"] <- ifelse(abs(scale(ds_temp[, "response"])) > i, NA, ds_temp[, "response"])
}
# exclude response time outliers
if (!is.na(j)) {
ds_temp[, "rt"] <- ifelse(abs(scale(ds_temp[ "rt"])) > j, NA, ds_temp[, "rt"])
}
# exclude additional outliers [not yet implemented]
# compute effect size
results <- c(as.numeric(unlist(cor.test(ds_temp[, "anchor"], ds_temp[, "response"]))[c(4, 9, 10)]), i, j, l) # , k
# attach results to dataset
cors <- rbind(cors
, results
)
}
}
}
names(cors) <- c("r"
, "rlo"
, "rhi"
, "exc_value"
, "exc_rt"
, "study")
cors
## recode studies
cors$study <- dplyr::recode(cors$study, "1" = "1: Item Valuations"
, "2" = "2: Item Valuations (Controlling for Visual Complexity)"
, "3" = "3: Factual Questions"
, "4" = "4: Legal Judgments"
, "5" = "5: Legal Judgments (Larger Replication)"
, "6" = "6: ???"
, "7" = "7: Explicit Anchoring in Item Valuations")
ggplot(cors, aes(x = r, xmin = rlo, xmax = rhi, y = as.character(study), col = as.factor(exc_value), fill = as.factor(exc_rt))) + geom_point() +
geom_errorbar() + facet_wrap(~exc_value + exc_rt) +
geom_vline(xintercept = 0) +
scale_y_discrete(limits=rev)
# Test of excess significance ---------------------------------------------
library(metafor) # package needed for TES
# dataset for tes (anchoring w. comparative question is not included)
dst <- cohens_ds[1:5, ]
# test if observed number of significant findings is larger than expected number of significant findings
metafor::tes(x = dst$d
, sei = dst$se
, data = dst
, alternative = "greater"
, test = "exact" # default an recommended if number of tests is less equal 10^6
)