### 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