################################################################################# ## "Can the Behavioral Sciences Self-Correct? A Social Epistemic Study". ## Studies in History and Philosophy of Science, 2016 ## Simulations Code ## ## Feel free to contact me if you have questions or troubles making this code work. ## ## Felipe Romero ## f.romero@uvt.nl ## ################################################################################# # Loads the required packages and setup empty variables. loadPackages<-function() { # Uncomment to Install required packages #install.packages("MBESS") #install.packages("metafor") #install.packages("pwr") #install.packages("plotrix") #install.packages("tizkDevice") #install.packages("ggplot2") library(lattice) library(pwr) library(MBESS) library(metafor) library(plotrix) library(tikzDevice) library(ggplot2) } # intializes variables for an iteration run setup<-function() { # Number of simulations per scenario # 500 iterations per scenario takes quite some time. Reduce this value for testing purposes sims<<-500 #number of replications attempted nreplications<<-0 #variables to keep track of the evolution of meta-analysis Rd<<-NULL Rn<<-NULL Rlci<<-NULL Ruci<<-NULL met<<-NULL MOE<<-NULL mlci<<-NULL muci<<-NULL } # Simulates replications of an experiment until scientists achieve 100 publications # effectSize: real effect size # publicationRule # ALL : All experiments are published # SIG : Only statistically significant results p<0.05 are published. # sampleSize : replication sample size # metaSource: Sources for meta-analysis # ALL : All the published experiments # RECENT: Most recent five experiments # bias : Only results of positive magnitude are published simulateCommunity<-function(effectSize, publicationRule, sampleSize, metaSource, bias) { D0<<-effectSize publications<-0 while (publications<50) { #replication sample size rn<-sampleSize #degrees of freedom rdf<-2*rn-2 #replication non-centrality parameter rncp<-D0 * 1/sqrt(2/rn) #replication random deviate from t distribution rept<-rt(1,rdf,rncp) #replication observed Cohen's d rdb<- rept * sqrt(1/rn + 1/rn) #unbiased estimate of Cohen's d using Hedges (1981) formula rd<-(1-3/(4*rdf-1))*rdb #replication observed non-centrality parameter roncp<-rd * 1/sqrt(2/rn) #replication lower bound of the 95% confidence interval rlci<-ci.smd(smd=rd,n.1=rn,n.2=rn)$Lower.Conf.Limit.smd #replication upper bound of the 95% confidence interval ruci<-ci.smd(smd=rd,n.1=rn,n.2=rn)$Upper.Conf.Limit.smd #increase the replications count nreplications<<-nreplications+1 #initially the result is not published published<-FALSE # Decide whether the study is published or not according to bias and selected publication rule if (publicationRule=="ALL" && bias ==FALSE) { published<-TRUE publications<-publications+1 } else if (publicationRule=="SIG" && bias == FALSE) { # Compute p-value if (2*pt(-abs(rept),rdf)<0.025) { published<-TRUE publications<-publications+1 } } else if (publicationRule=="ALL" && bias == TRUE) { if (rd>0) { published<-TRUE publications<-publications+1 } } else if (publicationRule=="SIG" && bias == TRUE) { if (2*pt(-abs(rept),rdf)<0.025 && rd>0) { published<-TRUE publications<-publications+1 } } # Add published study to meta-analysis if (published==TRUE) { # Variables to compute cumulative meta-analytic effect #list of replication sample sizes of all the experiments done so far Rn<<-c(Rn,rn) #list of effect sizes of experiments done so far Rd<<-c(Rd,rd) #confidence intervals of experiments done so far Rlci<<-c(Rlci,rlci) Ruci<<-c(Ruci,ruci) #compute the variances of effect sizes of experiments so far v<-variancefromd(Rd,Rn) if (metaSource=='ALL') { #compute and store the meta-analytical result of the experiments done so far m<-rma(Rd,v,method='FE') met<<-c(met,m$b[1]) mlci<<-c(mlci, m$ci.lb) muci<<-c(muci, m$ci.ub) } else if (metaSource=='RECENT') { #compute meta-analytical result of recent publications only #window of publications window<-10 inbound<-length(Rd)-window upbound<-length(Rd) if (inbound<0) { inbound<-0 } m<-rma(Rd[inbound:upbound],v[inbound:upbound],method='FE') met<<-c(met,m$b[1]) mlci<<-c(mlci, m$ci.lb) muci<<-c(muci, m$ci.ub) } } } } #computes variance from Cohen's d and sample size n variancefromd<-function(d,n) { v<-(2/n) + (d^2/(2*n^2)) return(v) } ################# # PLOTS ################# # Plots a t-distribution # Corresponds to Figures 1, 2, and 3 in the paper. plotT<-function() { #tikz('fragilityfig3.tex', width = 3.25, height = 3.25) alpha=.05 #effect size d<-0.41 #sample size sample_size<-156 #degrees of freedom Df<-2*sample_size-2 #non-centrality parameter ncp<-d/sqrt(2/sample_size) #plot density of distribution using a smooth line over histogram xfit<-seq(fro=-4,to=8,by=.01) yfit<-dt(xfit,Df, ncp) plot(0,0, xlim=c(-2,6), ylim=c(0,0.45), main='', xlab='t', lwd=1, type='n', ylab='P(t)', cex=0.5) lines(xfit,yfit,col='black',lwd=1.5) lower_critical_t<-qt(.5*alpha,2*sample_size-2,lower.tail=T) upper_critical_t<-qt(1-.5*alpha,2*sample_size-2,lower.tail=T) #plot critical values abline(v=c(lower_critical_t,upper_critical_t),col = "black", lwd=1) # dev.off() } # Plots the simulation of experiments in one society. # Figures 4 plotFigure4<-function() { setup() simulateCommunity(0.41,"ALL",156,"ALL", FALSE) #uncomment for exporting plots in tex format #tikz('fragilityfig4.tex', width = 3.25, height = 3.25) fr<-data.frame(x=c(1:10), fit=Rd[1:10], lower=Rlci[1:10], upper=Ruci[1:10]) ggplot(fr, aes(x=factor(x), y=fit))+ theme_grey()+ theme(legend.key = element_rect(colour = 0),axis.title = element_text(size = rel(0.8)), axis.text = element_text(size = rel(0.8)),legend.position="bottom",legend.direction="horizontal")+ geom_point(size=1)+ geom_errorbar(aes(ymax=upper, ymin=lower), width=0.4)+ labs(y = "Effect Size (d)")+ labs(x = "Experiment Number")+ geom_hline(aes(yintercept=D0, linetype="Real effect size", color="Real effect size"))+ scale_colour_manual("",labels=c("Real effect size"),values=c("Real effect size"="red"))+ scale_linetype_manual("",labels=c("Real effect size"), values=c("Real effect size"=2))+ ylim(-0.1,0.9) #dev.off() } # Plots the simulation of meta-analysis in one society. plotFigure5<-function() { #tikz('fragilityfig5.tex',width = 3.25, height = 3.25) fr<-data.frame(x=c(1:10), fit=met[1:10], lower=mlci[1:10], upper=muci[1:10]) ggplot(fr, aes(x=factor(x), y=fit))+ theme_grey()+ theme(legend.key = element_rect(colour = 0),axis.title = element_text(size = rel(0.8)), axis.text = element_text(size = rel(0.8)),legend.position="bottom",legend.direction="horizontal")+ geom_point(size=1)+ geom_errorbar(aes(ymax=upper, ymin=lower), width=0.4)+ labs(y = "Effect Size (d)")+ labs(x = "Meta-Analysis Number")+ geom_hline(aes(yintercept=D0, linetype="Real effect size", color="Real effect size"))+ scale_colour_manual("",labels=c("Real effect size"),values=c("Real effect size"="red"))+ scale_linetype_manual("",labels=c("Real effect size"), values=c("Real effect size"=2))+ ylim(-0.1,0.9) #dev.off() } # Plots the simulation of experiments in one society when there is no real effect # Figures 10 plotFigure10<-function() { setup() #tikz('fragilityfig10.tex', width = 3.25, height = 3.25) simulateCommunity(0,"SIG",156,"ALL", FALSE) fr<-data.frame(x=c(1:50), fit=Rd[1:50], lower=Rlci[1:50], upper=Ruci[1:50]) ggplot(fr, aes(x=x, y=fit))+ theme_grey()+ theme(legend.key = element_rect(colour = 0),axis.title = element_text(size = rel(0.8)), axis.text = element_text(size = rel(0.8)),legend.position="bottom",legend.direction="horizontal")+ geom_point(size=1)+ geom_errorbar(aes(ymax=upper, ymin=lower), width=0.4)+ labs(y = "Effect Size (d)")+ labs(x = "Experiment Number")+ geom_hline(aes(yintercept=D0, linetype="Real effect size", color="Real effect size"))+ scale_colour_manual("",labels=c("Real effect size"),values=c("Real effect size"="red"))+ scale_linetype_manual("",labels=c("Real effect size"), values=c("Real effect size"=2)) #dev.off() } # Plots for scenarios S1-S4 # Figure 6 plotsAllPublishedFig6<-function() { S1<-robustnessWrapperPlot(0.41,"ALL",156,"ALL", FALSE,sims) S1$suff<-"Sufficient resources" S1$bias<-"No direction bias" S2<-robustnessWrapperPlot(0.41,"ALL",36,"ALL", FALSE,sims) S2$suff<-"Limited resources" S2$bias<-"No direction bias" S3<-robustnessWrapperPlot(0.41,"ALL",156,"ALL", TRUE,sims) S3$suff<-"Sufficient resources" S3$bias<-"Direction bias" S4<-robustnessWrapperPlot(0.41,"ALL",36,"ALL", TRUE,sims) S4$suff<-"Limited resources" S4$bias<-"Direction bias" D0<-0.41 fr<-rbind(S1,S2,S3,S4) fr$suff <- factor(fr$suff, levels = rev(sort(unique(fr$suff)))) fr$bias <- factor(fr$bias, levels = rev(sort(unique(fr$bias)))) #tikz('fragilityfig6.tex',width = 7.1, height = 5,pointsize=10) robustnessWrapperGrid("S1","S2","S3","S4", fr) #dev.off() } # Plots for scenarios S5-S8 # Figure 7 plotsAllPublishedFig7<-function() { S5<<-robustnessWrapperPlot(0,"ALL",156,"ALL", FALSE,sims) S5$suff<-"Sufficient resources" S5$bias<-"No direction bias" S6<<-robustnessWrapperPlot(0,"ALL",36,"ALL", FALSE,sims) S6$suff<-"Limited resources" S6$bias<-"No direction bias" S7<<-robustnessWrapperPlot(0,"ALL",156,"ALL", TRUE,sims) S7$suff<-"Sufficient resources" S7$bias<-"Direction bias" S8<<-robustnessWrapperPlot(0,"ALL",36,"ALL", TRUE,sims) S8$suff<-"Limited resources" S8$bias<-"Direction bias" D0<-0 fr<-rbind(S5,S6,S7,S8) fr$suff <- factor(fr$suff, levels = rev(sort(unique(fr$suff)))) fr$bias <- factor(fr$bias, levels = rev(sort(unique(fr$bias)))) #tikz('fragilityfig7.tex',width = 7.1, height = 5,pointsize=10) robustnessWrapperGrid("S5","S6","S7","S8", fr) #dev.off() } # Plots for scenarios S9-S12 # Figure 8 plotsAllPublishedFig8<-function() { S9<<-robustnessWrapperPlot(0.41,"SIG",156,"ALL", FALSE,sims) S10<<-robustnessWrapperPlot(0.41,"SIG",36,"ALL", FALSE,sims) S11<<-robustnessWrapperPlot(0.41,"SIG",156,"ALL", TRUE,sims) S12<<-robustnessWrapperPlot(0.41,"SIG",36,"ALL", TRUE,sims) S9$suff<-"Sufficient resources" S9$bias<-"No direction bias" S10$suff<-"Limited resources" S10$bias<-"No direction bias" S11$suff<-"Sufficient resources" S11$bias<-"Direction bias" S12$suff<-"Limited resources" S12$bias<-"Direction bias" D0<-0.41 fr<-rbind(S9,S10,S11,S12) fr$suff <- factor(fr$suff, levels = rev(sort(unique(fr$suff)))) fr$bias <- factor(fr$bias, levels = rev(sort(unique(fr$bias)))) #tikz('fragilityfig8.tex',width = 7.1, height = 5,pointsize=10) robustnessWrapperGrid("S9","S10","S11","S12", fr) #dev.off() } # Plots for scenarios S13-S16 # Figure 9 plotsAllPublishedFig9<-function() { S13<<-robustnessWrapperPlot(0,"SIG",156,"ALL", FALSE,sims) S14<<-robustnessWrapperPlot(0,"SIG",36,"ALL", FALSE,sims) S15<<-robustnessWrapperPlot(0,"SIG",156,"ALL", TRUE,sims) S16<<-robustnessWrapperPlot(0,"SIG",36,"ALL", TRUE,sims) S13$suff<-"Sufficient resources" S13$bias<-"No direction bias" S14$suff<-"Limited resources" S14$bias<-"No direction bias" S15$suff<-"Sufficient resources" S15$bias<-"Direction bias" S16$suff<-"Limited resources" S16$bias<-"Direction bias" D0<-0 fr<-rbind(S13,S14,S15,S16) fr$suff <- factor(fr$suff, levels = rev(sort(unique(fr$suff)))) fr$bias <- factor(fr$bias, levels = rev(sort(unique(fr$bias)))) #tikz('fragilityfig9.tex',width = 7.1, height = 5,pointsize=10) robustnessWrapperGrid("S13","S14","S15","S16", fr) #dev.off() } #constructs a grid with four plots based on a data frame that contains them all # n1-n4 names of the plots robustnessWrapperGrid<-function(n1,n2,n3,n4, fr) { ggplot(fr, aes(x, fit,lower,upper))+ theme_grey()+ theme(legend.key = element_rect(colour = 0),axis.title = element_text(size = rel(0.8)), axis.text = element_text(size = rel(0.8)),legend.position="bottom",legend.direction="horizontal")+ geom_ribbon(aes(ymin=lower, ymax=upper, color="Standard dev.", linetype="Standard dev.",fill="Standard dev."), alpha=0.2)+ geom_line(aes(y = fit,color='Mean observed effect', linetype="Mean observed effect", fill="Mean observed effect"))+ geom_hline(aes(yintercept=D0, color="Real effect size",linetype="Real effect size",fill="Real effect size") )+ scale_colour_manual("",labels=c("Mean observed effect", "Real effect size","Standard dev."),values=c("Mean observed effect"="blue", "Real effect size"="red","Standard dev."="transparent"))+ scale_linetype_manual("",labels=c("Mean observed effect", "Real effect size","Standard dev."),values=c("Mean observed effect"=1, "Real effect size"=2, "Standard dev."=1))+ scale_fill_manual("",labels=c("Mean observed effect", "Real effect size","Standard dev."),values=c("grey", "grey", "grey12"))+ # scale_fill_manual("", values=c("Standard dev."="grey12"))+ facet_grid(bias~suff) + labs(y = "Effect Size (d)")+ labs(x = "Meta-Analysis Number")+ geom_text(data=data.frame(x=Inf, y=Inf, label=c(n1,n2), suff=c("Sufficient resources","Limited resources"), bias=c("No direction bias", "No direction bias") ), aes(x,y,label=label,label.padding=0.5), inherit.aes=FALSE, vjust=1.8,hjust=1.8,size=3)+ geom_text(data=data.frame(x=Inf, y=Inf, label=c(n3,n4), suff=c("Sufficient resources","Limited resources"), bias=c("Direction bias", "Direction bias")), aes(x,y,label=label), inherit.aes=FALSE,vjust=1.8,hjust=1.8,size=3) } #runs the selected scenario multiple times to compute expected values. # numCommunities: number of communities to be simulated robustnessWrapperPlot<-function(effect, rule, sample, sourcem, bias, numCommunities) { nreplicationsall<<-NULL metall<<-NULL for (i in 1:numCommunities) { setup() if (i%%100==0) { print(i) } simulateCommunity(effect,rule,sample,sourcem, bias) metall<<-rbind(metall,met) nreplicationsall<<-c(nreplicationsall,nreplications) } means<-NULL sds<-NULL for (i in 1:50) { means<-c(means,mean(metall[,i])) sds<-c(sds,sd(metall[,i])) } low<-means-sds upp<-means+sds S<-data.frame(x=c(1:50), fit=means, lower=low, upper=upp) return(S) }