Created
April 30, 2016 16:29
-
-
Save Lakens/88f7046b6b22a6eab827c79dfa82f0c2 to your computer and use it in GitHub Desktop.
Internal meta-analysis on 4 studies, 50% of which are true effects
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment