Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment