### clear the workspace
rm(list = ls())

### install packages if not already installed
install.packages("zcurve")
install.packages("readxl")

### load the zcurve package to create a zcurve
library(zcurve)

### load the excel package to read excel files
library(readxl)

### you can find many datasets here 
"https://www.dropbox.com/s/wj82c0wy7f08gpb/"

### in theory you should be able to load data directly from the web
### but I am having problems with this file 
### for now, just download the files 
#link = "https://www.dropbox.com/s/wj82c0wy7f08gpb/Ulrich%20Schimmack-z.wos.csv?dl=0"
#data = read.csv(url(link))
#dim(data)

### set the directory to where the data are stored
setwd("c:\\Users\\ulric\\Dropbox\\PHPCurves\\Meta-Meta\\Author")

### get a list of file names in the directory 
file.names = list.files(full.names = FALSE);  
file.names = file.names[pos = unlist(gregexpr(pattern = "z.wos.csv",file.names,ignore.case=TRUE)) > 0]
file.names

author = "Ulrich Schimmack"
file.name = file.names[pos = unlist(gregexpr(pattern = author,file.names,ignore.case=TRUE)) > 0]
file.name

### read the data from csv file
data = data.frame(read.csv(file.name))

### read the data from xlsx file
data = data.frame(readxl(file.name))

### see how many rows and columns there are 
dim(data)

### read the column/variable names
colnames(data)

### see how many rows have pvalues vs. no data (NA)
table(!is.na(data$z))

### remove rows with missing data 
data = data[!is.na(data$z),]

### if the data are p-values, transform p-values to z-values
### using log.p allows to convert even very small p-values 
data$z = qnorm(log(1 - data$p/2), log.p=TRUE)

### replace very high z-scores and infinity (Inf) with a maximum value of your choice (any value greater than 6 is ok)
data$z[data$z > 10] = 10

### get summary statistics for the transformed z-scores
summary(data$z)


### conduct z-curve analysis and store results in a data frame "results"

### if numer of z-scores (k) is large, run first without bootstrap
results = zcurve(data$z, bootstrap = 0); plot.zcurve(results); summary(results)

### show full z-curve including predictd distribution for non-significant results 
plot.zcurve(results,extrapolate = TRUE)


### if numer of z-scores (k) is small or for final analysis, run with boostrap default minimum 500
results = zcurve(data$z, bootstrap = 500); plot.zcurve(results); summary(results)

### show full z-curve including predictd distribution for non-significant results and confidence interval  
plot.zcurve(results,extrapolate = TRUE, CI = TRUE)

### add annotation with key results to the plot
plot.zcurve(results,extrapolate = TRUE, CI = TRUE, annotation = TRUE)


### compute Observed Discovery Rate (ODR); percentage of significant results
alpha = .05
odr = table(data$z > qnorm(1-alpha/2)); odr = odr/sum(odr); odr

### compute the false discovery risk 
fdr = round((1/coef(summary(results))[2,]-1)*(alpha/(1-alpha))*100)[c(1,3,2)]; fdr

### fit z-curve with different alpha level, lower alpha to decrease false positive risk 
alpha = .01
### lower.limit.of.fitted.interval
ll.int = 1.96
results = zcurve(data$z, bootstrap = 500, control=list(a = ll.int, sig_level = alpha)); plot.zcurve(results); summary(results)
fdr = round((1/coef(summary(results))[2,]-1)*(alpha/(1-alpha))*100)[c(1,3,2)]; fdr
plot.zcurve(results,extrapolate = TRUE, CI = TRUE, annotation = TRUE)
odr = table(data$z > qnorm(1-alpha/2)); odr = odr/sum(odr); odr