Skip to content

Commit c9f945c

Browse files
committed
Code clean up
- add 7-methodDominance.R - move all secondary data files to /dataFiles - deleted old/unnecessary files
1 parent b142c59 commit c9f945c

23 files changed

+182
-210
lines changed

.gitignore

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ vignettes/*.pdf
1414

1515
analysisParts/*.RData
1616

17-
# for the moment: do not upload simParts (1.2GB)
17+
# for the moment: do not upload simParts (1.2GB overall)
1818
simParts
1919

2020
# shell script for automatically adding files > 50 MB:
@@ -23,3 +23,8 @@ simParts
2323
res.final.RData
2424
res.wide.RData
2525
res.wide.red.RData
26+
.git/objects/pack/pack-04e34f809f18eeb4d75e1cc330cb24b24fd26ce0.pack
27+
dataFiles/res.final.RData
28+
dataFiles/res.wide.RData
29+
dataFiles/res.wide.red.RData
30+

1-simFramework.R

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,8 @@
1+
## ======================================================================
2+
## This file generates the simulated data sets, which are stored in separate
3+
## files in folder /simParts
4+
## ======================================================================
5+
16
# run this file:
27
# source("1-simFramework.R", echo=TRUE)
38

@@ -19,7 +24,7 @@ k_set <- c(10, 30, 60, 100) # number of studies in each MA
1924
delta_set <- c(0, .2, .5, .8) # true mean of effect sizes
2025
qrpEnv_Set <- c("none", "med", "high") # QRP environment
2126
selProp_set <- c(0, .6, .9) # publication bias
22-
tau_set <- c(0, .2, .4) # heterogeneity
27+
tau_set <- c(0, .2, .4) # heterogeneity; assumed to follow a normal distribution
2328

2429
# params stores all possible combinations of experimental factors
2530
params <- expand.grid(k=k_set, delta=delta_set, qrpEnv=qrpEnv_Set, selProp=selProp_set, tau=tau_set)
@@ -81,17 +86,14 @@ for (j in 1:nrow(param)) {
8186
res0 <- cbind(
8287
batch = batch,
8388
replication = i,
84-
condition = j,
89+
condition = j,
8590

8691
# settings of the condition
8792
p,
8893

8994
# results of the computation
9095
as.matrix(MA1))
9196

92-
# collect results in the matrix
93-
#res[counter:(counter+nrow(MA1)-1), ] <- res0
94-
#counter <- counter+nrow(MA1)
9597
res <- rbind(res, res0)
9698
} # of b-loop
9799

@@ -102,9 +104,6 @@ for (j in 1:nrow(param)) {
102104
sim <- sim %>% mutate(id=1000*(batch*10^(floor(log10(max(replication))+1)) + replication) + condition)
103105
save(sim, file=paste0("simPartsDemo/simData_condition_", j, ".RData"), compression="gzip")
104106

105-
# send a push notification after each finished condition:
106-
# userkey <- "uY7zyarxM2HoNaTLeX8HXjWvpFA4Cp" #Define user key
107-
# send_push(userkey, paste0("Condition ", j, "/", nrow(params), " finished"))
108107
} # of j (loop through parameter combinations)
109108

110109

2-analysisFramework.R

Lines changed: 3 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,6 @@ simDatFiles <- list.files("simParts", pattern=".*\\.RData", full.names=TRUE)
2121

2222
library(gtools)
2323
simDatFiles <- mixedsort(simDatFiles)
24-
# f <- simDatFiles[[33]]
2524

2625

2726
# loop through all simParts files
@@ -66,11 +65,11 @@ for (f in simDatFiles) {
6665
pc_skew(t=MAdat$t, df=MAdat$N-2, long=TRUE),
6766
pcurveEst(t=MAdat$t, df=MAdat$N-2, progress=FALSE, long=TRUE, CI=FALSE),
6867
puniformEst(t.value=MAdat$t, n1=MAdat$n1, n2=MAdat$n2),
68+
TPSM.est(t=MAdat$t, n1=MAdat$n1, n2=MAdat$n2, long=TRUE)#,
6969
#topN(MAdat$d, MAdat$v, MAdat$n1, MAdat$n2, est="fixed", fixed.effect=0.3),
7070
#topN(MAdat$d, MAdat$v, MAdat$n1, MAdat$n2, est="rma"),
71-
#topN(MAdat$d, MAdat$v, MAdat$n1, MAdat$n2, est="PEESE"),
72-
TPSM.est(t=MAdat$t, n1=MAdat$n1, n2=MAdat$n2, long=TRUE),
73-
betaSM.est(d=MAdat$d, v=MAdat$v, long=TRUE)
71+
#topN(MAdat$d, MAdat$v, MAdat$n1, MAdat$n2, est="PEESE"),
72+
#betaSM.est(d=MAdat$d, v=MAdat$v, long=TRUE)
7473
)
7574

7675

@@ -92,10 +91,4 @@ for (f in simDatFiles) {
9291
} # of dopar
9392

9493
save(res, file=paste0("analysisParts/analysis_", basename(f)), compress="gzip")
95-
96-
# send a push notification after each 50 finished conditions:
97-
if (which(simDatFiles == f) %% 50 == 0) {
98-
userkey <- "uY7zyarxM2HoNaTLeX8HXjWvpFA4Cp" #Define user key
99-
send_push(userkey, paste0("Condition ", which(simDatFiles == f), " finished"))
100-
}
10194
} # of "f in simDatFiles"

3-resultsFramework.R

Lines changed: 8 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -26,8 +26,8 @@ res.final <- bind_rows(res_list)
2626
str(res.final)
2727

2828
# final data set in long format:
29-
save(res.final, file="res.final.RData")
30-
#load(file="res.final.RData")
29+
save(res.final, file="dataFiles/res.final.RData")
30+
#load(file="dataFiles/res.final.RData")
3131

3232

3333
# Show conditions
@@ -58,8 +58,8 @@ print(tab2, n=54)
5858

5959
res.wide <- inner_join(res.wide, tab2)
6060

61-
save(res.wide, file="res.wide.RData", compress="gzip")
62-
#load(file="res.wide.RData")
61+
save(res.wide, file="dataFiles/res.wide.RData", compress="gzip")
62+
#load(file="dataFiles/res.wide.RData")
6363

6464
# ---------------------------------------------------------------------
6565
# save a reduced version that applies some selections
@@ -75,9 +75,6 @@ res.wide.red[res.wide.red$method == "3PSM" & is.na(res.wide.red$b0_p.value), c("
7575
## RULE 3: Ignore p-uniform when it doesn't provide a lower CI (very rare cases)
7676
res.wide.red[res.wide.red$method == "puniform" & is.na(res.wide.red$b0_conf.low), c("b0_estimate", "b0_conf.low", "b0_conf.high", "b0_p.value")] <- NA
7777

78-
# TODO: Skip this??
79-
# RULE X: set pcurve and puniform estimates to NA for all conditions which have less than 500/1000 successful meta-analyses
80-
#res.wide.red[res.wide.red$method %in% c("pcurve.evidence", "pcurve.hack", "pcurve.lack", "pcurve", "puniform") & res.wide.red$nMA.with.kSig.larger.3 < 500, c("b0_estimate", "b0_conf.low", "b0_conf.high", "b0_p.value", "skewtest_p.value")] <- NA
8178

8279
# ---------------------------------------------------------------------
8380
# For hypothesis test: Add H0.rejection rule
@@ -97,8 +94,8 @@ res.wide.red$H0.reject[is.na(res.wide.red$p.value)] <- NA
9794
res.wide.red$H0.reject.wrongSign <- (res.wide.red$p.value < .05) & (is.na(res.wide.red$b0_estimate) | res.wide.red$b0_estimate < 0)
9895
res.wide.red$H0.reject.wrongSign[is.na(res.wide.red$p.value)] <- NA
9996

100-
save(res.wide.red, file="res.wide.red.RData", compress="gzip")
101-
#load(file="res.wide.red.RData")
97+
save(res.wide.red, file="dataFiles/res.wide.red.RData", compress="gzip")
98+
#load(file="dataFiles/res.wide.red.RData")
10299

103100

104101
# ---------------------------------------------------------------------
@@ -135,8 +132,8 @@ print(summ, n=50)
135132

136133
# summ contains the full summary of the simulations. This object can then be used to build tables, plots, etc.
137134
library(rio)
138-
export(summ, file="summ.csv")
139-
save(summ, file="summ.RData")
135+
export(summ, file="dataFiles/summ.csv")
136+
save(summ, file="dataFiles/summ.RData")
140137

141138
# also export into Shiny app
142139
save(summ, file="Shiny/MAexplorer/summ.RData")

4-EstimationPlot.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ library(dplyr)
88
library(grid)
99
library(gridExtra)
1010

11-
load("summ.RData")
11+
load("dataFiles/summ.RData")
1212

1313
# ---------------------------------------------------------------------
1414
# Plot settings

5-hypTestPlot.R

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,7 @@ library(grid)
88
library(gridExtra)
99

1010
#setwd("C:/Users/evan.c.carter/Documents/Meta-analysis showdown")
11-
load(file="summ.RData")
11+
load(file="dataFiles/summ.RData")
1212

1313
# ---------------------------------------------------------------------
1414
# SETTINGS
@@ -34,7 +34,7 @@ hyp.wide <- inner_join(hyp.H0, hyp.H1)
3434
hyp.wide$rejectionRatio <- hyp.wide$Power/hyp.wide$TypeI
3535
hyp.wide$errorSum <- (1-hyp.wide$Power) + hyp.wide$TypeI
3636

37-
save(hyp.wide, file="hyp.wide.RData")
37+
save(hyp.wide, file="dataFiles/hyp.wide.RData")
3838
save(hyp.wide, file="Shiny/MAexplorer/hyp.wide.RData")
3939

4040

6-coveragePlot.R

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ library(dplyr)
77
library(grid)
88
library(gridExtra)
99

10-
load(file="summ.RData")
10+
load(file="dataFiles/summ.RData")
1111

1212
# ---------------------------------------------------------------------
1313
# SETTINGS

7-methodDominance.R

Lines changed: 139 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,139 @@
1+
library(ggplot2)
2+
library(dplyr)
3+
library(grid)
4+
library(gridExtra)
5+
6+
load(file="dataFiles/summ.RData")
7+
head(summ)
8+
9+
s2 <- summ %>%
10+
ungroup() %>%
11+
select(condition, method, ME.pos, RMSE.pos) %>%
12+
filter(!method %in% c("pcurve.evidence", "pcurve.hack", "pcurve.lack", "PET.rma", "PEESE.rma", "PETPEESE.rma"))
13+
14+
# get experimental factors
15+
conditions <- s2 %>% filter(method=="reMA") %>% select(1:6)
16+
17+
18+
# ---------------------------------------------------------------------
19+
# pairwise dominance
20+
21+
# reshape results into a data set with binary comparisons
22+
# This probably could be made more elegant, but I couldn't figure out how.
23+
24+
n.methods <- length(unique(s2$method))
25+
methods <- unique(s2$method)
26+
27+
res <- data.frame()
28+
for (C in unique(s2$condition)) {
29+
print(C)
30+
for (i in 1:n.methods) {
31+
for (j in 1:n.methods) {
32+
if (i > j)
33+
res <- rbind(res, data.frame(
34+
condition = C,
35+
method1 = methods[i],
36+
method2 = methods[j],
37+
ME.pos1 = as.numeric(s2[s2$condition == C & s2$method==methods[i], "ME.pos"]),
38+
ME.pos2 = as.numeric(s2[s2$condition == C & s2$method==methods[j], "ME.pos"]),
39+
RMSE.pos1 = as.numeric(s2[s2$condition == C & s2$method==methods[i], "RMSE.pos"]),
40+
RMSE.pos2 = as.numeric(s2[s2$condition == C & s2$method==methods[j], "RMSE.pos"])
41+
))
42+
}
43+
}
44+
}
45+
46+
save(res, file="dataFiles/dominanceScore.RData")
47+
#load(file="dataFiles/dominanceScore.RData")
48+
49+
res2 <- res %>% na.omit()
50+
51+
# define points, winners and draws
52+
53+
res2$winner <- ""
54+
res2$draw1 <- ""
55+
res2$draw2 <- ""
56+
57+
winner1 <- (abs(res2$ME.pos1) < abs(res2$ME.pos2)) & (res2$RMSE.pos1 < res2$RMSE.pos2)
58+
res2$winner[winner1] <- as.character(res2$method1[winner1])
59+
60+
winner2 <- (abs(res2$ME.pos1) > abs(res2$ME.pos2)) & (res2$RMSE.pos1 > res2$RMSE.pos2)
61+
res2$winner[winner2] <- as.character(res2$method2[winner2])
62+
63+
draw <- ((abs(res2$ME.pos1) > abs(res2$ME.pos2)) & (res2$RMSE.pos1 < res2$RMSE.pos2)) | ((abs(res2$ME.pos1) < abs(res2$ME.pos2)) & (res2$RMSE.pos1 > res2$RMSE.pos2))
64+
res2$draw1[draw] <- as.character(res2$method1[draw])
65+
res2$draw2[draw] <- as.character(res2$method2[draw])
66+
67+
68+
## give 1 point for "draw"
69+
# scores <- res2 %>% group_by(condition) %>% summarise(
70+
# RE = sum(winner=="reMA")*2 + sum(draw1=="reMA") + sum(draw2=="reMA"),
71+
# TF = sum(winner=="TF")*2 + sum(draw1=="TF") + sum(draw2=="TF"),
72+
# PET = sum(winner=="PET.lm")*2 + sum(draw1=="PET.lm") + sum(draw2=="PET.lm"),
73+
# PEESE = sum(winner=="PEESE.lm")*2 + sum(draw1=="PEESE.lm") + sum(draw2=="PEESE.lm"),
74+
# PETPEESE = sum(winner=="PETPEESE.lm")*2 + sum(draw1=="PETPEESE.lm") + sum(draw2=="PETPEESE.lm"),
75+
# pcurve = sum(winner=="pcurve")*2 + sum(draw1=="pcurve") + sum(draw2=="pcurve"),
76+
# puniform = sum(winner=="puniform")*2 + sum(draw1=="puniform") + sum(draw2=="puniform"),
77+
# TPSM = sum(winner=="3PSM")*2 + sum(draw1=="3PSM") + sum(draw2=="3PSM")
78+
# )
79+
80+
scores <- res2 %>% group_by(condition) %>% summarise(
81+
RE = sum(winner=="reMA")*2 ,
82+
TF = sum(winner=="TF")*2 ,
83+
PET = sum(winner=="PET.lm")*2,
84+
PEESE = sum(winner=="PEESE.lm")*2,
85+
PETPEESE = sum(winner=="PETPEESE.lm")*2,
86+
pcurve = sum(winner=="pcurve")*2 ,
87+
puniform = sum(winner=="puniform")*2 ,
88+
TPSM = sum(winner=="3PSM")*2
89+
)
90+
91+
scores.long <- melt(scores, id.vars="condition")
92+
93+
scores$winner <- ""
94+
for (i in 1:nrow(scores)) {
95+
scores$winner[i] <- paste0(colnames(scores)[which(scores[i, 2:9] == max(scores[i, 2:9])) + 1], collapse = ", ")
96+
}
97+
98+
scores <- merge(scores, conditions, by="condition")
99+
scores.long <- merge(scores.long, conditions, by="condition")
100+
101+
scores$winner <- factor(scores$winner, levels=names(sort(table(scores$winner), decreasing=TRUE)))
102+
103+
sort(table(scores$winner), decreasing=TRUE)
104+
105+
ggplot(scores, aes(y=qrpEnv, x=factor(delta), shape=factor(winner))) + geom_point(size=6) + facet_grid(selProp~k+tau)
106+
107+
# That looks ugly ...
108+
scores.long$loop <- paste0(scores.long$delta, ":", scores.long$qrpEnv, ":", scores.long$tau)
109+
ggplot(scores.long, aes(y=value, x=loop, color=factor(variable), group=factor(variable))) + geom_line() + facet_grid(k~selProp)
110+
111+
# ---------------------------------------------------------------------
112+
# strong (absolute) dominance
113+
114+
s2 <- summ %>%
115+
ungroup() %>%
116+
select(condition, k, delta, qrpEnv, selProp, tau, method, ME.pos, RMSE.pos) %>%
117+
filter(!method %in% c("pcurve.evidence", "pcurve.hack", "pcurve.lack", "PET.rma", "PEESE.rma", "PETPEESE.rma"))
118+
119+
# which method is best in ME?
120+
ME.matrix <- dcast(s2, condition ~ method, value.var="ME.pos")
121+
ME.matrix$ME.winner <- colnames(ME.matrix)[apply(ME.matrix, 1, which.min)]
122+
123+
# which method is best in RMSE?
124+
RMSE.matrix <- dcast(s2, condition ~ method, value.var="RMSE.pos")
125+
RMSE.matrix$RMSE.winner <- colnames(RMSE.matrix)[apply(RMSE.matrix, 1, which.min)]
126+
127+
fullDominance <- cbind(ME.matrix[, c("condition", "ME.winner")], RMSE.winner = RMSE.matrix$RMSE.winner)
128+
fullDominance$winner <- ""
129+
fullDominance$winner[fullDominance$ME.winner == fullDominance$RMSE.winner] <- fullDominance$ME.winner[fullDominance$ME.winner == fullDominance$RMSE.winner]
130+
131+
sort(table(fullDominance$winner[fullDominance$winner != ""]), decreasing = TRUE)
132+
133+
s3 <- merge(fullDominance, conditions, by="condition")
134+
135+
136+
ggplot(s3 %>% filter(winner != ""), aes(y=qrpEnv, x=factor(delta), shape=factor(winner))) + geom_point(size=6) + facet_grid(selProp~k+tau)
137+
138+
139+
ggplot(s3 %>% filter(winner != ""), aes(y=qrpEnv, x=factor(delta), shape=factor(winner))) + geom_point(size=6) + facet_grid(selProp~k+tau)

Empirical n and ES distributions/fit_Fraley_2014.R

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -7,12 +7,9 @@ r <- rtrunc(1, spec="beta", a=.08, b=.14, 1.34, 5.03 )
77
dtrue <- 2*r/sqrt(1-r^2) # this will select effect sizes that hit the social psych "sweet spot." it runs from 20th to 40th percentiles. This is roughly from the mode up a bit.
88

99

10-
11-
1210
#### Get sample sizes. ####
1311
# I pulled the raw data from Fraley & Vazire's n-pact factor paper, then selected the social psych journals (JESP, the first 2 sections of JPSP, and SPPS)
1412
# these were total N per study. To be generous, I assumed that they were all 2 group designs, then used this function to randomly sample a per-group n from that distribution
1513

16-
1714
rtrunc(n=1, spec="nbinom", a=10, b=Inf, size=2.3, mu=48)
1815

Empirical n and ES distributions/tau.R

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,8 +16,12 @@ prop.table(table(sel$tau < .2))
1616
prop.table(table(sel$tau < .4))
1717

1818

19-
# Investigate I^2 across all meta-analyses
19+
# ---------------------------------------------------------------------
20+
# Investigate I^2 across all meta-analyses
21+
22+
# Compute I^2 from Q
2023
dat$I2 <- (dat$Q - (dat[, "# of effect sizes"] - 1))/dat$Q
2124
dat$I2[dat$Q <= (dat[, "# of effect sizes"] - 1)] <- 0
2225

23-
prop.table(table(dat$I2 > .25))
26+
prop.table(table(dat$I2 > .25))
27+
prop.table(table(dat$I2 > .75))

0 commit comments

Comments
 (0)