# http://DataColada.org/76
#
#Code to reproduce MTurk study on Maluma vs Takiti using Dutch Baby names from 2016
#Written by Uri Simonsohn (urisohn@gmail.com)
#Please contact directly if you see any errors or have any comments.
#
#
#**Note from 2019 04 29** --- After posting, and following requests by Blake McShane, #
# we added a dummy 1/0 to both datasets for whether the same MTurk ID
# was observed across waves, and we later added a masked MTurk id as well.
# When first posted, there was one observation that was was miscoded as being duplicate when it was not.
# (as could be verified by the non-matching masked MTurk id's).
# Blake McShane spotted the error whose origin we were unable to identify, as re-running the
# code we had used to generate the file, did not generate the error.
# The .csv files relied on the code below, which load the data from the Colada server, have that error corrected.
#
#Version: 2019 05 06 (typo fixed in dropping missing values)
#
###########################################################################################################
#Outline
# 0. Preliminaries
# 1. Functions (just two)
# 2. Get data ready to be analyzed
#3 Analysis.
#3.1 Function that computes d and n for every study, outputs data.frame with three columns: name, d, n
#3.2 Compute d & n for collected data
#3.3 Correlation in observed effect size across waves within study
#3.4 Meta analysis of the first wave, the I2=78% for Wave 1 reported in the post
#3.5 Figure 1
#3.6 Figure 2 - Scatterplot
#3.7 Bootstraping individual data rows (original, see 3.9 by bootstrapping by participant)
#3.8 Absolute mean error and Figure 3
#3.9 Bootstraping participants
#######################################################################
# 0. Preliminaries
rm(list = ls())
library(metafor) #to run meta-analysis
library(stringr) #to handle strings
library(rio) #to import any data format
#1) Functions
#1.1 Function 1: like vlookup in Excel for dealing with the names
vlookup=function(x,table,k=2) table[match(x,table[,1]),k]
#1.2 Function 2: format two-decimal numbers without a 0
round.fun <- function(val) { sub("^(-?)0.", "\\1.", sprintf("%.2f", val)) }
#2) Get data ready to be analyzed
#2.1) load directly from the Colada server
#These files are the raw Qualtrics .csv files, except the following modifications
#IP, latitude, longitude and MTURK id have been deleted for privacy reasons
#Load from server datafiles which include masked MTurk ids
a1=import("http://datacolada.org/appendix/76/Data%201%20-%20Colada%2076%20-%202019%2004%2026%20(masked%20mturk%20id).csv")
a2=import("http://datacolada.org/appendix/76/Data%202%20-%20Colada%2076%20-%202019%2004%2026%20(masked%20mturk%20id).csv")
#2.1.5) Add a subject id, just a counter as the true MTurk id is hidden for privacy reasons (this code was written before the masked id had been generated or used)
a1$id=1:nrow(a1)
a2$id=1:nrow(a2)
#Note: 2019 04 25
#Late (on European time) on the day of the post we received an email, from Blake McShane, asking some clarifying questions
#about the data, including if there were any repeat participant. We went back to the original files which ahd the
#MTUrk id that we had dropped, and created new columns in both datasets for when the same MTurk ID was used
#by un-commenting teh following two lines, you can re-run the entire code dropping repeat participants
#we added a footnote to the post with those results
#if writing code from scratch we would have made it easier to carry out the exclusions without relying on commenting out, but we didn't.
#a1=a1[a1$duplicate.mturk.id==0,] #drop participants from Wave 1 that also participated in Wave 2 (at least, same MTUrk id)
#a2=a2[a2$duplicate.mturk.id==0,] #drop participants from Wave 2 that also participated in Wave 1 (at least, same MTUrk id)
#2.2 Reshape the data to grab the multiple columns with values and put them in a single vector, plus the name
#They come in default Qualtrics format, where each question version is its own column, which are mostly empty
#see the variable names
names(a1) #Data starts in column 5, so below the loop will start grabbing ys there.
y1=y2=names1=names2=c()
for (k in 1:50)
{
#Grab the names being used
y1=c(y1,a1[,4+k]) #4+k becuase as noted above, starting in column 5 we have the data
y2=c(y2,a2[,4+k])
names1=c(names1,rep(names(a1)[4+k],nrow(a1)))
names2=c(names2,rep(names(a2)[4+k],nrow(a2)))
}
#2.2.1 Create participant id for every row (copy vertically all the ids 50 times), this will allow knowing which person we are drawing later on in the bootstrap
#Note 2019 04 29: This would no long ber relevant as we have masked MTurk id, but it is the way it was coded originally
id1=rep(a1$id,50)
id2=rep(a2$id+1000,50)
#2.3 Delete empty rows
keep1=ifelse(is.na(y1),F,T)
keep2=ifelse(is.na(y2),F,T)
names1=names1[keep1]
names2=names2[keep2]
y1=y1[keep1]
y2=y2[keep2]
id1=id1[keep1]
id2=id2[keep2]
#2.4 Turn to numeric
y1=as.numeric(y1)
y2=as.numeric(y2)
#y2=as.numeric(y2)+runif(length(y2)) #use this to make sure the process works, when run this way, adding noise, indeed expected variation is below what's observed
#2.5 Sort names so girl name appears first always
#Create table with names and whether it is reverse-code (i.e. girl name in the right)
reverse.raw=rep(c(0,1),25) #Names in the original stimuli are reversed girl-boy,boy-girl, etc.
names.raw=names(a1[5:54]) #grab the names
names_split=str_split(names.raw,"-") #split the first and second name
#Turn "list" int otwo arrays
first = unlist(names_split)[2*(1:50)-1] #the name appearing first
second = unlist(names_split)[2*(1:50) ] #the name appearing second
#put back together, always girl-boy
names.adj=ifelse(reverse.raw==0,paste(first," - ",second),paste(second, " - ",first))
names.adj
#Now table names
reverse.table=data.frame(names.raw, reverse.raw, names.adj)
#3 Analysis.
#3.1 Function that computes d and n for every study, outputs data.frame with three columns: name, d, n
run.once=function(y,names)
#y: is the dv
#nmes: is a vector with the name pair associated to each observation
{
#M and SD by name
m =aggregate(y,list(names),mean)
sd=aggregate(y,list(names),sd)
sum1=aggregate(y,list(names),FUN=sum)
n=sum1$x/m$x
#Lookup the names to know if they need to be reverse coded
reverse=vlookup(m$Group.1,reverse.table,2) #uses custom formula 1, like vlookup() in Excel
#Adjust
m.adjusted=ifelse(reverse==1,6-m$x,m$x) #Reverse coded when girl name is on the right
#Cohen d
d=(m.adjusted-3)/sd$x #standardize the effect size
#output the cohen d and the sample size for each name
data.frame(m$Group.1, d, n)
}
#3.2 Compute d & n for collected data
df1.obs=run.once(y1,names1) #this outpus a data.frame with the name pairs, cohen d, and sample size, in wave 1
df2.obs=run.once(y2,names2) #this is a data.frame with the name pairs, cohen d, and sample size, in wave 1
#3.3 Correlation in observed effect size across waves
r.obs=cor(df1.obs$d,df2.obs$d) #each row is a name-pair, we compute the correlation between them
r.obs
#3.4 Meta analysis of the first wave, the I2=78% for Wave 1 reported in the post
rma(df1.obs$d,1/df1.obs$n) #This runs a meta-analysis with d as the effect size and 1/N as the variance
#it's 1/N because d has Mean=0 and Var=1 by construction.
#The mean of that variable, then has 1/N as variance
#3.5 Figure 1
# setwd("c:/dropbox/DataColada/76 Heterogeneity") #Change to folder you want to use, if you will be saving with png()
#Effect in Wave 1
d1=df1.obs$d
#Names re-coded in Wave 1, using vlookup (Function 2 above)
n1=vlookup(df1.obs$m.Group.1[order(d1)],reverse.table,3)
# png("F1 - Barplot for Maluma Takiti.png", width=1200,height=800,res=125) #commented out, this saves the graph to a file
#Margin
par(mar=c(10.5,5.1,4.1,2.1))
#Plot
bp1=barplot(sort(d1),col='red3',las=1,ylim=c(-.1,.8))
#x-axis
axis(side=1,at=bp1,n1,las=2,cex=.9)
mtext(side=1,line=8.5,cex=1.5,font=2,"Girl - Boy name pair")
#y-axis
mtext(side=2,line=3,cex=1.5,font=2,"Cohen's d")
#mtext(side=2,line=3,cex=1.1,font=3,"(Rating i girl name to Maluma - 3)/SD")
#Value labels
dx1=df1.obs$d[order(df1.obs$d)]
text(bp1,dx1+.02*sign(dx1),round.fun(dx1),cex=.7)
mtext(side=3,font=2,cex=1.55,"Standardized rating assigning the girl name to the rounded shape ")
mtext(side=3,font=3,cex=1.15,line=-1.25,"(Standardized: SDs of 1-5 rating above midpoint)")
#dev.off() #Commented out, this completes the .png file and closes it
#3.6 Figure 2 - Scatterplot
#Simple vectors to work with
d1=df1.obs$d #effect size in wave 1
d2=df2.obs$d #effect size in wave 2
#png("F2 - Scatterplot Maluma Takiti.png", width=1200,height=800,res=125) #Commented out, this save figure to file
par(mar=c(5.1,6.1,4.1,2.1))
plot(d1, d2,pch=16,col='gray77', xlab="",ylab="",xlim=c(-.2,1.1),ylim=c(-.2,1.1))
abline(0,1,lty=3,lwd=3)
abline(lm(d2~d1),col='red3',lwd=3)
#legend
# Perfectly Replicable (identical estimate waves 1 & 2)
# Observed Replicability (OLS line)
legend(-.2,1.1,col=c("black","red3"),lty=c(3,1), lwd=c(3,3),
legend=c("Perfect replicability line (45 degrees: Wave 2 = Wave 1)",
"Observed replicability line (regression: Wave 2 = a + b * Wave 1)"),cex=1.15)
#x-axis
mtext(side=1,font=2,cex=1.4,"Effect Size in Wave 1",line=2.5)
mtext(side=1,font=3,cex=1.2,"(June 2017)",line=3.75)
#y-axis
mtext(side=2,font=2,cex=1.4,"Effect Size in Wave 2",line=3.75)
mtext(side=2,font=3,cex=1.2,"(April 2019)",line=2.55)
#Text for correlation
cor(d1,d2)
cor(d1,d2,method="spearman") #robustness, not reported
text(.6,0,"corr(wave 1, wave2) = .77",col='gray44',cex=1.25)
#Header
mtext(side=3,line=2,cex=2,font=2,"Observably Heterogenous Effects Do Replicate")
mtext(side=3,line=.2,cex=1.35,font=3,col='gray20',"(Name-pairs showing larger Maluma effect in 2017, show larger effect in 2019)")
#dev.off() #Commented out, this would save the .png and close the figure file
#3.7 Bootstraping individual data rows (original, see 3.8 by bootstrapping by participant)
#Preliminaries
set.seed(190)
simtot=5000 #Number of bootstrap iterations
#Combine y1 and y2, and names, to draw under null they all come from same population, no heterogeneity
y12=c(y1,y2)
names12=c(names1,names2)
#Start the loop, saving the correlations in rk1 (when drawing only from 1) and rk12 (when drawing from 1 and 2)
rk1=rk12=c()
d1.all=d2.all=matrix(nrow=simtot,ncol=50)
for (simk in 1:simtot)
{
N1=length(y1) #How many observations in Wave 1
N2=length(y2) #How many observations in Wave 2
#3.7.1 Draw from Wave 1 only (in footnote)
k1=sample(N1,size=N1,replace=T) #Draw Wave 1 in the size that it was, from Wave 1
k2=sample(N1,size=N2,replace=T) #Draw Wav2 2 in the size that it was, from Wave 1
y1.k1=y1[k1] #This is boostrapped Wave 1
y1.k2=y1[k2] #This is boostrapped Wave 2
names1.k1=names1[k1] #The corresponding names
names1.k2=names1[k2]
df1.bootk1=run.once(y1.k1,names1.k1) #Compute the cohen d's
df2.bootk2=run.once(y1.k2,names1.k2)
rk1[simk]=cor(df1.bootk1$d, df2.bootk2$d) #Store the correlation between d1 and d2 in this sample
#3.7.2 Draw from waves 1 and 2 combined (In main text)
k1=sample(N1+N2,size=N1, replace=T) #Draw these observations for Wave 1 from 1 and 2
k2=sample(N1+N2,size=N2, replace=T) #Draw these observations for Wave 2 from 1 and 2
#DV (here is the actual drawing, above just which observations are selected)
y1.k12=y12[k1]
y2.k12=y12[k2]
#Names
names1.k12=names12[k1]
names2.k12=names12[k2]
#Compute cohen ds
df1.bootk12=run.once(y1.k12,names1.k12)
df2.bootk12=run.once(y2.k12,names2.k12)
#Compute correlations
rk12[simk]=cor(df1.bootk12$d, df2.bootk12$d)
#Save all ds
d1.all[simk,]=df1.bootk12$d
d2.all[simk,]=df2.bootk12$d
#Counter every 500
if (simk%%500==0) cat("...",simk)
}
#3.7.3 Results of the bootstrap
#3.7.3.1 Mean correlation, and p-value, bootstrapping from merging from waves
#See the result
mean(rk12) #Average correlation under the null
#Make the histogram in the footnotes
hist(rk12,col='gray77',border=F,breaks=35,ylim=c(0,550),
xlab='Correlation: r(Wave 1, wave 2))',
main="Distribution of Correlations between Wave 1 and Wave 2 under null")
#Vertical lines in histogram
abline(v=r.obs,col='red4',lwd=4)
abline(v=mean(rk12),col='gray44',lwd=3,lty=2)
text(.76,550,adj=1,"Observed correlation r=.77",col='red4')
text(.80,550,adj=0,"Expected correlation r=.79",col='gray44')
#3.7.3.2 Mean correlation, and p-value, bootstrapping from only wave 1
mean(rk1)
hist(rk1)
abline(v=r.obs,col='red',lwd=4)
#p-values for both appraochs
mean(abs(rk12) <= abs(r.obs))
mean(abs(rk1) <= abs(r.obs))
#3.8 Absolute mean error
#3.8.1 Compute absolute error in bootstrapped samples
#Computer average absolute deviation between waves in the bootstrapped data
dif.abs.null=abs(d1.all-d2.all) #d1.all and d2.all are data.frames generated in the boostrapping
#loop above, they contain simtot rows and 50 columns, with the estimates for each simulation
#Verify visually it works as intended: the 3rd row is the difference of first 2
rbind(d1.all[,1],d2.all[,1],dif.abs.null[,1])[,1:5]
#Expected mean absolute difference is teh mean of all absolute differences in the simtot bootstraps
mean.dif.null=mean(dif.abs.null)
#3.8.2 Compute absolute mean error in the observed sample
mean.dif.obs=mean(abs(d1-d2))
#3.8.3 Compute ranked absolute errors
#Sort each row in the matrix with difference in effect so that we can analyze ranks
sorted.rows=matrix(nrow=simtot, ncol=50) #This matrix has simtot rows, with all studies in each row sorted by the absolute error
for (k in 1:simtot) sorted.rows[k,]=sort(dif.abs.null[k,])
null.sorted=colMeans(sorted.rows)
#3.8.4 Plot figure 3
#png("F3 - Expected vs observed absolute errors.png", width=1400,height=800,res=125) #commented out, this saves the graph to a file
par(mar=c(13.1,6.1,4.1,2.1))
plot(sort(abs(d1-d2)),xaxt="n",pch=16,col='red4',cex=1.25,
xlab="",ylab="",las=1,xlim=c(0,55))
points(null.sorted,col='blue',type='l',lwd=2)
legend(2,.4,cex=1,lty=c(NA,1),pch=c(16,NA), col=c("red4","blue"),
legend=c("Observed difference between waves for each name-pair",
"Expected difference between waves under ZERO heterogeneity"))
axis(side=1,at=1:50,paste0(n1[order(abs(d1-d2))]," [",1:50,"]"),las=2,cex=.9)
#Names for bars
axis(side=1,at=53,"Observed average ",las=2,cex=1,font=2)
axis(side=1,at=55,"Expected average ",las=2,cex=1,font=2)
#Bars at the end for the aveage
polygon(x=c(52,52,54,54),y=c(0,mean.dif.obs,mean.dif.obs,0),col='red4')
polygon(x=c(54,54,56,56),y=c(0,mean.dif.null,mean.dif.null,0),col='blue3')
#Value labels
text(52.5,mean.dif.obs+.01,round(mean.dif.obs,3),cex=.9,col='red4')
text(55.5,mean.dif.null+.01,round(mean.dif.null,3),cex=.9,col='blue3')
mtext(side=1,line=11,font=2,cex=1.75,"The 50 studies sorted by absolute difference in effects between waves")
#Y-axis
mtext(side=2,font=2,cex=1.5,line=4,"Difference in effect size across waves")
mtext(side=2,font=3,cex=1.35,line=2.55,"(abs(Wave 2 - Wave1))",col='gray60')
#Vertical line
abline(v=51,col='gray66')
#Header
mtext(side=3,font=2,cex=1.8,line=2,
"Effect sizes between waves differ as expected with ZERO heterogeneity")
#dev.off()
##########################################################################
#3.9 Bootstraping participants
#Note: In feedback provided upon seeing a draft of this post, Blake, Ulf and Karsten commented
#that the bootstrap above does not take into account the dependency across
#participants, to address that issue, here the bootstrap is done at the participant level
#Preliminaries
set.seed(191)
simtot=100 #Number of bootstrap iterations
#Combine y1 and y2, and names, to draw under null they all come from same population, no heterogeneity
y12=c(y1,y2)
names12=c(names1,names2)
id12=c(id1,id2)
#Combine ids for participants for both waves
N1.subjects=length(unique(id1)) #how many participants in wave 1
N2.subjects=length(unique(id2)) #how many participants in wave 2
id12.subjects=unique(c(id1,id2)) #unique identifier of each particiopant in both waves
#Start the loop, saving the correlations in rk1 (when drawing only from 1) and rk12 (when drawing from 1 and 2)
rk.subject=c()
simtot=250
for (simk in 1:simtot)
{
#3.9.1 Draw from both waves, entire subject
k1.subject=sample(id12.subjects,size=N1.subjects,replace=T) #Draw the # of participants seen in Wave 1, from Wave 1 & 2
k2.subject=sample(id12.subjects,size=N2.subjects,replace=T) #Draw the # of participants seen in Wave 2, from Wave 1 & 2
#3.9.2 Could not find easy vector-based way to do this, so i loop of participant IDs drawn at random
# and get all the values of y that match that id
y1.boot=y2.boot=name1.boot=name2.boot=c()
for (k in k1.subject)
{
y1.boot=c(y1.boot, y12[id12==k])
name1.boot=c(name1.boot,names12[id12==k])
}
for (k in k2.subject)
{
y2.boot=c(y2.boot, y12[id12==k])
name2.boot=c(name2.boot,names12[id12==k])
}
df1.bootk1=run.once(y1.boot,name1.boot) #Compute the cohen d's
df2.bootk2=run.once(y2.boot,name2.boot) #Compute the cohen d's
rk.subject[simk]=cor(df1.bootk1$d, df2.bootk2$d) #Store the correlation between d1 and d2 in this sample
#Counter every 500
if (simk%%10==0) cat("...",simk)
}
#3.9.3 Results of the bootstrap
# Mean correlation, and p-value, bootstrapping from merging from waves
r.expected=mean(rk.subject) #Average correlation under the null, is very similar, r=.795 instead of .79
r.expected
hist(rk.subject,col='gray77',border=F,breaks=35,ylim=c(0,100),
xlab='Correlation: r(Wave 1, wave 2))',
main="Distribution of Correlations between Wave 1 and Wave 2 under null")
abline(v=r.obs,col='red4',lwd=4)
abline(v=mean(rk12),col='gray44',lwd=3,lty=2)
text(.76,550,adj=1,"Observed correlation r=.77",col='red4')
text(r.expected,550,adj=0,"Expected correlation r=.795",col='gray44')
#p-values
mean(abs(rk.subject) <= abs(r.obs))