Instantly share code, notes, and snippets.

# richarddmorey/funnel_no_pub_bias.R Created Jan 9, 2016

A demonstration of how we can get asymmetric funnel plots with no publication bias
 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")