Internal meta-analysis on 4 studies, 50% of which are true effects

4study_meta_50%_true_effects.R

if(!require(meta)){install.packages('meta')}
library(meta)

nSims <- 1000000 #number of simulated experiments
numberstudies<-4 # nSim/numberofstudies should be whole number
p <-numeric(nSims) #set up empty container for all simulated p-values
metapran <-numeric(nSims/numberstudies) #set up empty container for all simulated p-values for random effects MA
metapfix <-numeric(nSims/numberstudies) #set up empty container for all simulated p-values for fixed effects MA
heterog.p<-numeric(nSims/numberstudies) #set up empty container for test for heterogeneity
d <-numeric(nSims) #set up empty container for all simulated d's
n=50 #set sample size in each condition
D<-0.43 #set effect size for studies with true effects

for(i in 1:(nSims/2)){ #first half of simulations, true effect
  x<-rnorm(n = n, mean = D, sd = 1) #produce 100 simulated participants
  #with mean=100 and SD=20
  y<-rnorm(n = n, mean = 0, sd = 1) #produce 100 simulated participants
  #with mean=100 and SD=20
  z<-t.test(x,y, var.equal=TRUE) #perform the t-test
  p[i]<-z$p.value #get the p-value and store it
  d[i] <-  (((mean(x)-mean(y))) / sqrt((((n - 1)*((sd(x)^2))) +  ((n - 1)*((sd(y)^2))))/((n * 2) -2)))  
}
for(i in (nSims/2):nSims){ #for each simulated experiment
  x<-rnorm(n = n, mean = 0, sd = 1) #produce simulated participants
  y<-rnorm(n = n, mean = 0, sd = 1) #produce simulated participants
  z<-t.test(x,y, var.equal=TRUE) #perform the t-test
  p[i]<-z$p.value #get the p-value and store it
  d[i] <-  (((mean(x)-mean(y))) / sqrt((((n - 1)*((sd(x)^2))) +  ((n - 1)*((sd(y)^2))))/((n * 2) -2)))  #calculate d
}

J<-1-3/(4*(n+n-2)-1) #correction for bias

matd<-matrix(d, ncol=numberstudies, byrow=FALSE) #create columns for sets of numberofstudies
matd.se<-sqrt((((n+n)/(n*n))+(matd^2/(2*(n+n))))*J^2) #calc SE

for(i in 1:(nSims/numberstudies)){ #for each simulated experiment
  es.d<-matd[i,]*J #get effect sizes in each row, multiply by correction for bias
  d.se<-matd.se[i,]
  meta<-metagen(es.d, d.se) #perform meta-analysis
  metapran[i]<-meta$pval.random
  metapfix[i]<-meta$pval.fixed
  heterog.p[i]<-(1-pchisq(meta$Q, meta$df.Q)) #calculate p-value for heterogeneity test
}

(sum(metapran < 0.05)/(nSims/numberstudies)) #power random effects
(sum(metapfix < 0.05)/(nSims/numberstudies)) #power fixed effects

(sum(p < 0.05)/nSims) #power for all studies

(sum(heterog.p < 0.05)/nSims) #power test for heterogeneity


matp<-matrix(p, ncol=numberstudies, byrow=FALSE)
sum(apply(matp, 1, function(x) all(x < 0.05^(1/numberstudies))))/(nSims/numberstudies)#number of sets of studies both all studies significant