A demonstration of how we can get asymmetric funnel plots with no publication bias

funnel_no_pub_bias.R

Nmin = 20
Nmax = 100
Mmax = 100
Mtotal = 2*Nmin*Mmax

Num_sim_per_N = 10

N = rep(seq(Nmin,Nmax,10), each = Num_sim_per_N)
M = round( Mtotal / ( 2 * N ) )

eff.size=.2
sig.sub=1
sig.err=1

# For reproducibility of result
set.seed(1234)

sim.data = function(N,M,eff.size,sig.sub,sig.err){
  # Simulate data
  sub.eff = rnorm(N, 0, sig.sub)
  x = outer(sub.eff, c(eff.size, 0), "+") + rnorm(2*N, 0, sig.err/sqrt(M))
  
  # t test
  tstat = t.test(x[,1],x[,2],paired=TRUE)$statistic
  
  # effect size
  raw = mean(x[,1] - x[,2])
  sd_diff = sd(x[,1] - x[,2])
  serr_diff = sd_diff / sqrt(N)
  corr_scores = cor(x[,1],x[,2])
  # Cohen's d, no within correction
  d1 = raw/sd_diff
  # Cohen's d, within correction
  sd_within = sd_diff / sqrt( 2 * ( 1 - corr_scores ) )
  d2 = raw/sd_within
  J = 1 - 3 / (4 * (N-1) - 1)
  # Hedge's g, no within correction
  g1 = J*d1
  # Hedge's g, within correction
  g2 = J*d2
  return(c(t=tstat, raw=raw, serr_diff=serr_diff, d1=d1, d2=d2, g1=g1, g2=g2))
}

results = mapply(sim.data, N=N, M=M, 
                 MoreArgs = list(eff.size=eff.size,
                                 sig.sub=sig.sub,
                                 sig.err=sig.err)
                 )


sim_results = data.frame(N=N,M=M,t(results))

plot(x = sim_results$g1,
     y = sim_results$N, 
     pch=19,col=rgb(1,0,0,.5),ylab="Sample size",xlab="Hedge's g")

plot(x = sim_results$raw,
     y = sim_results$N, 
     pch=19,col=rgb(1,0,0,.5),ylab="Sample size",xlab="Raw effect")