Why p-values are not measures of evidence

p_not_evidence.R

library(metafor)

g <- escalc(
  measure = "SMD",
  n1i = 43, # sample size in group 1 is 50
  m1i = 1.3, # observed mean in group 1 is 5.6
  sd1i = 1, # observed standard deviation in group 1 is 1.2
  n2i = 43, # sample size in group 2 is 53
  m2i = 1, # observed mean in group 1 is 4.9
  sd2i = 1
) # observed standard deviation in group 2 is 1.3

# print results
meta_data <- rbind(g, g, g) 

# Perform the meta-analysis
res <- rma(yi, vi, data = meta_data)
res

forest(res)


jpeg(file="plot2.jpg",width=2000,height=1400, units = "px", res = 300)
forest(res)
dev.off()

# Lindley plot

N <- 150
p <- 0.05
p_upper <- 0.05 + 0.00000001
p_lower <- 0.00000001
ymax <- 50 # Maximum value y-scale (only for p-curve)

# Calculations

# p-value function
pdf2_t <- function(p) 0.5 * dt(qt(p / 2, 2 * N - 2, 0), 2 * N - 2, ncp) / dt(qt(p / 2, 2 * N - 2, 0), 2 * N - 2, 0) + dt(qt(1 - p / 2, 2 * N - 2, 0), 2 * N - 2, ncp) / dt(qt(1 - p / 2, 2 * N - 2, 0), 2 * N - 2, 0)

jpeg(file="plot1.jpg",width=3400,height=2000, units = "px", res = 300)

plot(-10,
     xlab = "P-value", ylab = "", axes = FALSE,
     main = "P-value distribution for d = 0, 50% power, and 99% power", xlim = c(0, 1), ylim = c(0, ymax), cex.lab = 1.5, cex.main = 1.5, cex.sub = 1
)
axis(side = 1, at = seq(0, 1, 0.05), labels = formatC(seq(0, 1, 0.05),format="f", digits=2), cex.axis = 1)
#Draw null line
ncp <- (0 * sqrt(N / 2)) # Calculate non-centrality parameter d
curve(pdf2_t, 0, 1, n = 1000, col = "grey", lty = 1, lwd = 2, add = TRUE)
#Draw 50% low power line
N <- 146
d <- 0.23
se <- sqrt(2 / N) # standard error
ncp <- (d * sqrt(N / 2)) # Calculate non-centrality parameter d
curve(pdf2_t, 0, 1, n = 1000, col = "black", lwd = 3, add = TRUE)
#Draw 99% power line
N <- 150
d <- 0.5
se <- sqrt(2 / N) # standard error
ncp <- (d * sqrt(N / 2)) # Calculate non-centrality parameter d
curve(pdf2_t, 0, 1, n = 1000, col = "black", lwd = 3, lty = 3, add = TRUE)

dev.off()