#Data Colada 103 - Mediation Analysis is Counterintuitively Invalid (http://datacolada.org/103)
#
# 2022 09 21
# Uri Simonsohn (urisohn@gmail.com)
#
#-------------------------------------------------------
#This script simulates a single dataset where X is randomly assigned and
#has an effect on Y which is NOT mediated by M, but mediation analysis
#concludes that it is 100% mediated by it.
#1 Clear
rm(list = ls())
#2 per-cell sample size
n=500
#3 Seed
set.seed(1955)
#4 True parameters
a=1 # Effect of manipulation on mediator
b=0 # Effect of mediator on the DV (change to see what happes when there is mediation)
c.prime = .5 # Effect of manipulation on DV apart from the mediator ('direct' effect)
#X: Quiz randomly assigned
x = rep(c(0,1),n)
#Omitted variable: love of math
Z = rnorm(2*n)
#M: mediator, how much each person studies
m = a*x + Z + rnorm(2*n) #People treated to quizzes, and who love math, study more
#Y: Final test
y = b*m+c.prime*x+ Z +rnorm(n*2) #b: effect caused by m (default b=0)
#c.prime: direct effect
#Z is the ommited variable
#Mediation analysis
lm1=lm(m~x) #[1] #unbiased
lm2=lm(y~x) #[2] #unbiased
lm3=lm(y~x+m) #[3] #BIASED
#See results
summary(lm1)
summary(lm2)
summary(lm3)
#Extract key coefficients to print results on console
a_est = summary(lm1)$coefficients[2,1]
c_est = summary(lm2)$coefficients[2,1]
c.prime_est = summary(lm3)$coefficients[2,1]
b_est = summary(lm3)$coefficients[3,1]
#Show results on console
cat("True a=",a,", estimated a=",a_est)
cat("True b=",b,", estimated b=",b_est)
cat("Total mediation=",a*b,", estimated a*b=",a_est*b_est)
cat("Estimated share mediated: ", round(1 - (c_est - (a_est*b_est))/c_est,0)*100,"%")