First you need to download the datafile from
https://github.com/chartgerink/2016statcheck_data/blob/master/statcheck_dataset.csv

Right click on <Raw> and save file.

When you are done, provide Path where R can find the file.

# Provide Path
GetPath =  <path>

# give file name
fn = “statcheck_dataset.csv”

# read datafile
d = read.csv(paste0(GetPath,fn))

# get t-values
t = d$Value
t[d$Statistic != “t”] = 0
summary(t)

#convert t-values into absolute z-scores
z.val.t = qnorm(pt(abs(t),d$df2,log.p=TRUE),log.p=TRUE)
z.val.t[z.val.t > 20] = 20
z.val.t[is.na(z.val.t)] = 0
summary(z.val.t)
hist(z.val.t[z.val.t < 6 & z.val.t > 0],breaks=30)
abline(v=1.96,col=”red”,lwd=2)
abline(v=1.65,col=”red”,lty=3)

#get F-values
F = d$Value
F[d$Statistic != “F”] = 0
F[F > 20] = 400
summary(F)

#convert F-values into absolute z-scores
z.val.F = qnorm(pf(abs(F),d$df1,d$df2,log.p=TRUE),log.p=TRUE)
z.val.F[z.val.F > 20] = 20
z.val.F[z.val.F < 0] = 0
z.val.F[is.na(z.val.F)] = 0
summary(z.val.F)
hist(z.val.F[z.val.F < 6 & z.val.F > 0],breaks=30)
abline(v=1.96,col=”red”,lwd=2)
abline(v=1.65,col=”red”,lty=3)

#get z-scores and convert into absolute z-scores
z.val.z = abs(d$Value)
z.val.z[d$Statistic != “Z”] = 0
z.val.z[z.val.z > 20] = 20
summary(z.val.z)
hist(z.val.z[z.val.z < 6 & z.val.z > 0],breaks=30)
abline(v=1.96,col=”red”,lwd=2)
abline(v=1.65,col=”red”,lty=3)

#check results
summary(cbind(z.val.t,z.val.F,z.val.z))

#get z-values for t,F, and z-tests
z.val = z.val.t + z.val.F + z.val.z

#check median absolute z-score by test statistic
tapply(z.val,d$Statistic,median)

##### save as r data file and reuse

### run analysis for specific author

# provide an author name as it appears in the authors column of the data file
author = “Stapel”
author.found = grepl(pattern = author,d$authors,ignore.case=TRUE)
table(author.found)

# give title for graphic
Name = paste0(“StatCheck “,author)

# select z-scores of author
z.val.sel = z.val[author.found == “TRUE”]

#set limit of y-axis of graph
ylim = 1

#create histogram
hist(z.val.sel[z.val.sel < 6 & z.val.sel > 0],xlim=c(0,6),ylab=”Density”,xlab=”|z| scores”,breaks=30,ylim=c(0,ylim),freq=FALSE,main=Name)
#add line for significance
abline(v=1.96,col=”red”,lwd=2)
#add line for marginal significance
abline(v=1.65,col=”red”,lty=3)

#compute median observed power
mop = median(z.val[z.val > 2 & z.val < 4])
# add bias to move normal distribution of fitted function to match observed distribution
bias = 0
# change variance of normal distribution to model heterogeneity (1 = equal power for all studies)
hetero = 1

### add fitted model curve to the plot
par(new=TRUE)
curve(dnorm(x,mop-bias,hetero),0,6,col=”red”,ylim=c(0,ylim),xlab=””,ylab=””)

### when satisfied with fit compute power
power = length(z.val.sel[z.val.sel > 1.96 & z.val.sel < 4]) * pnorm(mop – bias,1.96) + length(z.val.sel[z.val.sel > 4])
power = power / length(z.val.sel[z.val.sel > 1.96])

### add power estimate to the figure
text(3,.8,pos=4,paste0(“Power: “,round(power,2)))