# Set the working directory. CHANGE THIS TO THE FOLDER WITH THE DRYAD STUFF IN IT setwd("~/Desktop/Dryad stuff") # Load the pre-made datasets (see R script "R script to filter and classify the p value data") experimental.equals<-read.csv("Exact p value dataset.csv")[,-1] prop.sig<-read.csv("Proportion of significant p values per paper.csv")[,-1] experimental.equals\$blind <- as.factor(experimental.equals\$blind) prop.sig\$blind <- as.factor(prop.sig\$blind) ########### DESCRIPTIVE STATISTICS # How many papers are there in the two datasets? length(unique(experimental.equals\$paperID)) # 7644 length(unique(prop.sig\$ID)) # 12710 # And how many p values are there in each? nrow(experimental.equals) #55688 sum(prop.sig[,1]) + sum(prop.sig[,2]) # 126060 # What proportion of papers are blind in the two datasets? uniqueIDs <- unique(experimental.equals\$paperID); tally <- rep(0, length(uniqueIDs)) for(i in 1:length(uniqueIDs)) tally[i] <- as.numeric(experimental.equals\$blind[experimental.equals\$paperID==uniqueIDs[i]][1]) sum(tally - 1) / length(tally) # 14.82% are blind in the p equals dataset xx<-tapply(prop.sig\$blind, prop.sig\$blind, length); xx[2] / (xx[1]+xx[2]) # 13.38% are blind in the prop.sig dataset # List the sample size (in number of papers) and proportion of papers that are blind for each FoR category (z score dataset) # This information is in Figure 2 mean.z.per.paper <- tapply(experimental.equals\$z, experimental.equals\$paperID, mean) mean.z.per.paper <- data.frame(z = as.numeric(mean.z.per.paper), paperID = names(mean.z.per.paper), num.authors = as.numeric(tapply(experimental.equals\$num.authors, experimental.equals\$paperID, function(x) x[1])), year = as.numeric(tapply(experimental.equals\$year, experimental.equals\$paperID, function(x) x[1])), category = as.character(tapply(experimental.equals\$category, experimental.equals\$paperID, function(x) as.character(x[1]))), blind = as.numeric(tapply(experimental.equals\$blind, experimental.equals\$paperID, function(x) x[1]))) mean.z.per.paper\$blind <- mean.z.per.paper\$blind - 1 cat.z<-data.frame(tapply(mean.z.per.paper\$z, list(mean.z.per.paper\$category, mean.z.per.paper\$blind), length)) cat.z[is.na(cat.z)] <- 0; cat.z\$n <- cat.z[,1] + cat.z[,2]; cat.z\$prop <- cat.z[,2] / (cat.z[,1] + cat.z[,2]); cat.z <- cat.z[,3:4]; cat.z\$category <- rownames(cat.z); rownames(cat.z)<-NULL cat.z<-cat.z[with(cat.z, order(-n, -prop)), ] cat.z # List the sample size (in number of papers) and proportion of papers that are blind for each FoR category (prop.sig dataset) # This information is in Figure 3 cat.prop<-data.frame(tapply(prop.sig\$Sig, list(prop.sig\$category, prop.sig\$blind), length)) cat.prop[is.na(cat.prop)] <- 0; cat.prop\$n <- cat.prop[,1] + cat.prop[,2]; cat.prop\$prop <- cat.prop[,2] / (cat.prop[,1] + cat.prop[,2]); cat.prop <- cat.prop[,3:4]; cat.prop\$category <- rownames(cat.prop); rownames(cat.prop)<-NULL cat.prop<-cat.prop[with(cat.prop, order(-n, -prop)), ] cat.prop ########### STATISTICAL MODELS #### 1. Literature review of evolutionary biology evolution.data <- read.csv("Evolution literature review data.csv") evolution.data\$Blind <- as.factor(evolution.data\$Blind) # Is the effect sisze of the non-blind study higher than the blind one more often than expected? diff <- evolution.data\$g[evolution.data\$Blind==0] - evolution.data\$g[evolution.data\$Blind==1] binom.test(length(diff[diff > 0]), length(diff)) # Yes, p = 0.0152 # Descriptive statistics on the difference in effect size between non-blind and blind studies in each paie mean(diff) median(diff) sd(diff) / sqrt(length(diff)) # SE tapply(evolution.data\$Sample.size, evolution.data\$Blind, mean, na.rm=T) # Mean sample size of blind and non-blind studies tapply(evolution.data\$Sample.size, evolution.data\$Blind, function(x) sd(x[!is.na(x)])/sqrt(length(x[!is.na(x)]))) # Associated SE tapply(evolution.data\$Sample.size, evolution.data\$Blind, function(x) length(x[!is.na(x)])) # Number of sample sizes # Is there a difference in sample size between the blind and non-blind papers in the literature review? Answer: not significantly so (p = 0.0535) wilcox.test(evolution.data\$Sample.size[evolution.data\$Blind==0], evolution.data\$Sample.size[evolution.data\$Blind==1]) # Now a meta-analytical approach to the literature review (this takes account of the variance associated with each effect size - "var.g") library(metafor) meta<-rma.mv(yi = g, V = var.g., mods = ~Blind, random = ~Blind|PairID, method="REML", data = evolution.data) summary(meta) # the results # Reproduce Figure 1 (later cleaned up in Inkscape software) evol.means <- data.frame(Blind = names(tapply(evolution.data\$g, evolution.data\$Blind, mean)), g = as.numeric(tapply(evolution.data\$g, evolution.data\$Blind, mean)),PairID = "mean") ggplot(data=evolution.data) + coord_cartesian(ylim = c(-1.5, 6)) + geom_line(alpha=0.5, aes(x = Blind, y = g, colour = factor(PairID), group=PairID, size=1,alpha=0.5)) + geom_line(data=evol.means, aes(x = Blind, y = g, group=1), colour='black', size=2.5) +geom_point(aes(x = Blind, y = g, size=1/var.g.)) + theme_bw() ####### 2. Text mined p values # Load required packages library(arm) # needed for "standardize" function library(MuMIn) # needed for Model averaging # Set global option for NAs (needed for the call to 'dredge' below) options(na.action = "na.fail") ###### Analysis of the exact p-values (expressed as z scores) # Declare the full linear mixed model, for the z score dataset z.score.lmer <- standardize(lmer(z ~ blind*category + year + num.authors + I(num.authors^2) + (1|paperID), data = experimental.equals, REML=F)) # Compare all models, from the null to the full model, by their AICc scores ('dredge' builds the model selection table - look at it by typing "z.dredge" into the console) z.dredge<-dredge(z.score.lmer) # Perform model averaging on models with Akaike weight > 0.05 # Note - we use model averaging because several models have decent delta-AIC scores in the dredge output # These are the data in Table 1 model.avg(subset(z.dredge,weight>0.05))\$avg.model # Get the importance values for each parameters importance(subset(z.dredge,weight>0.05)) ###### Analysis of the proportion of significant p-values (expressed as z scores) # First, we need to "hack" the quasibinomial family function to make it produce the info needed for QAIC model ranking (this code was borrowed from the author of the MuMIn package, Kamil Barton) dfun <- function(object){with(object, sum((weights *residuals^2)[weights > 0]) / df.residual)} # Get deviance parameter, termed c-hat, from a model x.quasibinomial <- function(...) { # Modified quasibinomial family that has AIC attached res <- quasibinomial(...) res\$aic <- binomial(...)\$aic res} # Express the number of significant and non-significant results as a two-column matrix resp <- cbind(prop.sig\$Sig, prop.sig\$Not) # Specify the global model prop.sig.GLM <- standardize(glm(resp ~ blind*category + year + num.authors + I(num.authors^2), data = prop.sig, family = x.quasibinomial)) # Compare all models, from the null to the full model, by their QAIC scores (using the dispersion parameter of the global model) dredge.prop <- dredge(prop.sig.GLM, rank="QAIC", chat = dfun(prop.sig.GLM)) # Best model (dAIC is 4.58, weight = 0.91) # Because delta-AIC is so high for the best model, let's just look at the results of that one. It has all parameters except the blind x category interaction # These are the data in Table 2 summary(update(prop.sig.GLM, ~. -c.blind:category)) # Estimating the 95% confidence limits in Table 2 coeffs <- as.data.frame(summary(update(prop.sig.GLM, ~. -c.blind:category))\$coefficients) coeffs\$lowerCI <- coeffs[,1] - 1.96*coeffs[,2] coeffs\$upperCI <- coeffs[,1] + 1.96*coeffs[,2] ###### Finally, let's see what measured variable affect the probability that a paper is blind # Run a binomial GLM, with FoR category, year of publication, and author number as predictors, and blindness as the response variable summary(glm(blind~category+year+num.authors+ I(num.authors^2), family=binomial, data = prop.sig)) # Do a likelihood ratio test to get a p value for the effect of category mod1 <- glm(blind~category+year+num.authors+ I(num.authors^2), family=binomial, data = prop.sig) mod2 <- glm(blind~year+num.authors+ I(num.authors^2), family=binomial, data = prop.sig) anova(mod1, mod2, test="Chisq") # Highly significant