June 21, 2016 Filed in: Interactive | Psych | R | Science | Statistics | Teaching

“If you torture the data long enough, it will confess.”

This aphorism, attributed to Ronald Coase, sometimes has been used in a disrespective manner, as if it was wrong to do creative data analysis.

In fact, the art of creative data analysis has experienced despicable attacks over the last years. A small but annoyingly persistent group of second-stringers tries to denigrate our scientific achievements. They drag psychological science through the mire.

These people propagate stupid method repetitions; and what was once one of the supreme disciplines of scientific investigation – a creative data analysis of a data set – has been crippled to conducting an empty-headed step-by-step pre-registered analysis plan. (Come on: If I lay out the full analysis plan in a pre-registration, even an *undergrad* student can do the final analysis, right? Is that really the high-level scientific work we were trained for so hard?).

They broadcast in an annoying frequency that p-hacking leads to more significant results, and that researcher who use *p*-hacking have higher chances of getting things published.

What are the consequence of these findings? The answer is clear. Everybody should be equipped with these powerful tools of research enhancement!

Some researchers describe a performance-oriented data analysis as “data-dependent analysis”. We go one step further, and call this technique *data-optimal analysis (DOA)*, as our goal is to produce the optimal, most significant outcome from a data set.

I developed an **online app that allows to practice creative data analysis and how to polish your ***p*-values. It’s primarily aimed at young researchers who do not have our level of expertise yet, but I guess even old hands might learn one or two new tricks! It’s called “The p-hacker” (please note that ‘hacker’ is meant in a very positive way here. You should think of the cool hackers who fight for world peace). You can use the app in teaching, or to practice *p*-hacking yourself.

Please test the app, and give me feedback! You can also send it to colleagues: http://shinyapps.org/apps/p-hacker

The full R code for this Shiny app is on Github.

Here’s a quick walkthrough of the app. Please see also the quick manual at the top of the app for more details.

First, you have to run an initial study in the “New study” tab:

When you ran your first study, inspect the results in the middle pane. Let’s take a look at our results, which are quite promising:

Sometimes outlier exclusion is not enough to improve your result.

Now comes the magic. Click on the “Now: p-hack!” tab – **this gives you all the great tools to improve your current study**. Here you can fully utilize your data analytic skills and creativity.

In the following example, we could not get a significant result by outlier exclusion alone. But after adding 10 participants (in two batches of 5), controlling for age and gender, and focusing on the variable that worked best – voilà!

Do you see how easy it is to craft a significant study?

Now it is important to **show even more productivity**: Go for the next conceptual replication (i.e., go back to Step 1 and collect a new sample, with a new manipulation and a new DV). Whenever your study reached significance, click on the *Save* button next to each DV and the study is saved to your stack, awaiting some additional conceptual replications that show the robustness of the effect.

Honor to whom honor is due: Find the best outlet for your achievements!

My friends, let’s stand together and Make Psychological Science Great Again! I really hope that the *p*-hacker app can play its part in bringing psychological science back to its old days of glory.

Best regards,

PS: A similar app can be found on FiveThirtyEight: Hack Your Way To Scientific Glory

Comments (5) | Trackback

April 15, 2016 Filed in: R | Statistics

Several recent discussions on the Psychological Methods Facebook group surrounded the question whether an optional stopping procedure leads to biased effect size estimates (see also this recent blog post by Jeff Rouder). Optional stopping is a rather new technique, and potential users wonder about the potential down-sides, as these (out-of-context) statements demonstrate:

- “… sequential testing appears to inflate the observed effect size”
- “discussion suggests to me that estimation is not straight forward?“
- “researchers who are interested in estimating population effect sizes should not use […] optional stopping“
- “we found that truncated RCTs provide biased estimates of effects on the outcome that precipitated early stopping” (Bassler et al., 2010)

Hence, the concern is that the usefulness of optional stopping is severely limited, because of this (alleged) bias in parameter estimation.

Here’s a (slightly shortened) scenario from a Facebook discussion:

Given the recent discussion on optional stopping and Bayes, I wanted to solicit opinions on the following thought experiment.Researcher A collects tap water samples in a city, tests them for lead, and stops collecting data once a t-test comparing the mean lead level to a “safe” level is significant atp<.05. After this optional stopping, researcher A computes a Bayesian posterior (with weakly informative prior), and reports the median of the posterior as the best estimate of the lead level in the city.Researcher B collects the same amount of water samples but with a pre-specified N, and then also computes a Bayesian estimate.Researcher C collects water samples from every single household in the city (effectively collecting the whole population).Hopefully we can all agree that the best estimate of the mean lead level in the city is obtained by researcher C. But do you think that the estimate of researcher B is closer to the one from researcher C and should be preferred over the estimate of researcher A? What – if anything – does this tell us about optional stopping and its influence on Bayesian estimates?

Let’s simulate the scenario (R code provided below) with the following settings:

- The true lead level in the city has a mean of 3 with a SD of 2
- The “safe” lead level is defined at 2.7 (or below)

We run 10,000 studies with strategy A and save the sample mean of the lead level along with the final sample size. We run 10,000 studies with strategy B (sample sizes matched to those of the 10,000 A-runs) and save the sample mean of the lead level along with the sample size.

In contrast to the quoted scenario above, I will *not* compute a Bayesian posterior, because the usage of a prior *will* bias the estimate. (A side note: When using a prior, this bias is deliberately accepted, because a small bias is traded in for a reduction in variance of the estimates, as extreme and implausible sample estimates are shrunken towards more realistic numbers). Here, we simply take the plain sample mean, because this is an unbiased estimator – at least in the typical textbook-case of fixed sample sizes. But what happens with optional stopping?

For strategy A, we compute the mean across all significant Monte Carlo simulations (which were ~11%), which is 1.42.

This is much less than the true value of 3! When we look at the 10,000 fixed-*n* studies with the same sample sizes, we get a mean lead level of 3.00, which is exactly the true value.

The impact of optional stopping seems devastating – it screws up my effect size estimates, and leads to an underestimation of the true lead level!

Does it really?

The naive analysis, however, ignores two crucial points:

**If effect sizes from samples with different sample sizes are combined, they must be meta-analytically weighted according to their sample size (or precision).**Optional stopping (e.g., based on*p*-values, but also based on Bayes factors) leads to a conditional bias: If the study stops very early, the effect size must be overestimated (otherwise it wouldn’t have stopped with a significant*p*-value). But early stops have a small sample size, and in a meta-analysis, these extreme early stops will get a small weight.**The determination of sample size (fixed vs. optional stopping) and the presence of publication bias are separate issues.**By comparing strategy A and B, two issues are (at least implicitly) conflated: A does optional stopping*and*has publication bias, as she only reports the result if the study hits the threshold. Non-significant results go into the file drawer. B, in contrast, has a fixed sample size, and reports all results,*without*publication bias. You can do optional stopping without publication bias (stop if significant, but also report result if you didn’t hit the threshold before reaching n_max). Likewise, if B samples a fixed sample size, but only reports trials in which the effect size is close to a foregone conclusion, it will be very biased as well.

As it turns out, the overestimations from early terminations are perfectly balanced by *under*estimations from late terminations (Schönbrodt, Wagenmakers, Zehetleitner, & Perugini, 2015). Hence, optional stopping leads to a *conditional* bias (i.e., conditional on early or late termination), but it is *unconditionally* unbiased.

Hence, let’s keep these factors separate in our analysis and look at the 2 (publication bias: yes vs. no) x 2 (optional stopping vs. fixed sample size) x 2 (naive average vs. meta-analytically weighted average) combinations.

For this purpose, we have to update the simulation:

- The strategies A and B
*without publication bias*report all outcomes - Strategy A
*with*publication bias reports only the studies, which are significantly lower than the safe lead level - Strategy B
*with*publication bias reports only studies which show a sample mean which is smaller than the safe lead level (regardless of significance)

This is the sampling distribution of the sample means, across all 10,000 replications:

The distribution from strategy B (fixed-*n*) is well-behaved and symmetric. The distribution from strategy A (optional stopping) shows a bump at small effect sizes (these are the early stops with a small lead level).

Another way to look at this is to plot the single study estimates by sample size:

Early terminations in the sequential design underestimate the true level – but the late terminations at n=50 overestimate on average in the sequential design. This is the conditional bias – underestimation in early stops (because the optional stopping favored small lead levels), but overestimation in late stops. In the fixed-*n* design there is no conditional bias.

Here are the compute mean levels in our 8 combinations (true value = 3):

Sampling plan | PubBias | Naive mean | Weighted mean |
---|---|---|---|

sequential | FALSE | 2.85 | 3.00 |

fixed | FALSE | 3.00 | 3.00 |

sequential | TRUE | 1.42 | 1.71 |

fixed | TRUE | 2.46 | 2.55 |

If you selectively only report studies with a desired outcome (rows 3 & 4), the estimates cannot be trusted – all of them are way below the true value. Or, as Joachim Vandekerckhove put it: “I think it’s obvious that you can’t actively bias your data and expect magic to happen”.

If you report all studies (no publication bias), they must be properly weighted if they are combined. And then it does not matter whether sample sizes are fixed or optionally stopped! Both sampling plans lead to unbiased estimates.

(To be precise: it does not matter with respect to the unbiasedness of effect size estimates. It does matter concerning other properties, like the variance of estimates or the average sample size).

A more detailed analysis of the impact of sequential testing on parameter estimates can be found in our paper “Sequential Hypothesis Testing with Bayes Factors“. Finally, I want to quote a paragraph from our recent paper on Bayes Factor Design Analysis (Schönbrodt & Wagenmakers, 2016), which also summarizes the discussion and provides some more references:

Concerning the sequential procedures described here, some authors have raised concerns that these procedures result in biased effect size estimates (e.g., Bassler et al., 2010, J. Kruschke, 2014). We believe these concerns are overstated, for at least two reasons.First, it is true that studies that terminate early at the H1 boundary will, on average, overestimate the true effect. This conditional bias, however, is balanced by late terminations, which will, on average, underestimate the true effect. Early terminations have a smaller sample size than late terminations, and consequently receive less weight in a meta-analysis. When all studies (i.e., early and late terminations) are considered together, the bias is negligible (Berry, Bradley, & Connor, 2010; Fan, DeMets, & Lan, 2004; Goodman, 2007; Schönbrodt et al., 2015). Hence, the sequential procedure is approximately unbiased overall.Second, the conditional bias of early terminations is conceptually equivalent to the bias that results when only significant studies are reported and non-significant studies disappear into the file drawer (Goodman, 2007). In all experimental designs –whether sequential, non-sequential, frequentist, or Bayesian– the average effect size inevitably increases when one selectively averages studies that show a larger-than-average effect size. Selective publishing is a concern across the board, and an unbiased research synthesis requires that one considers significant and non-significant results, as well as early and late terminations.Although sequential designs have negligible unconditional bias, it may nevertheless be desirable to provide a principled “correction” for the conditional bias at early terminations, in particular when the effect size of a single study is evaluated. For this purpose, Goodman (2007) outlines a Bayesian approach that uses prior expectations about plausible effect sizes. This approach shrinks extreme estimates from early terminations towards more plausible regions. Smaller sample sizes are naturally more sensitive to prior-induced shrinkage, and hence the proposed correction fits the fact that most extreme deviations from the true value are found in very early terminations that have a small sample size (Schönbrodt et al., 2015).

library(ggplot2)

library(dplyr)

library(htmlTable)

# set seed for reproducibility

set.seed(0xBEEF)

trueLevel <- 3

trueSD <- 2

safeLevel <- 2.7

maxN <- 50

minN <- 3

B <- 10000 # number of Monte Carlo simulations

res <- data.frame()

for (i in 1:B) {

print(paste0(i, "/", B))

maxSample <- rnorm(maxN, trueLevel, trueSD)

# optional stopping

for (n in minN:maxN) {

t0 <- t.test(maxSample[1:n], mu=safeLevel, alternative="less")

#print(paste0("n=", n, "; ", t0$estimate, ": ", t0$p.value))

if (t0$p.value <= .05) break;

}

finalSample.seq <- maxSample[1:n]

# now construct a matched fixed-n

finalSample.fixed <- rnorm(n, trueLevel, trueSD)

# ---------------------------------------------------------------------

# save results in long format

# sequential design

res <- rbind(res, data.frame(

id = i,

type = "sequential",

n = n,

p.value = t0$p.value,

selected = t0$p.value <= .05,

empMean = mean(finalSample.seq)

))

# fixed design

res <- rbind(res, data.frame(

id = i,

type = "fixed",

n = n,

p.value = NA,

selected = mean(finalSample.fixed) <= safeLevel, # some arbitrary publication bias selection

empMean = mean(finalSample.fixed)

))

}

save(res, file="res.RData")

# load("res.RData")

# Figure 1: Sampling distribution

ggplot(res, aes(x=n, y=empMean)) + geom_jitter(height=0, alpha=0.15) + xlab("Sample size") + ylab("Sample mean") + geom_hline(yintercept=trueLevel, color="red") + facet_wrap(~type) + theme_bw()

# Figure 2: Individual study estimates

ggplot(res, aes(x=empMean)) + geom_density() + xlab("Sample mean") + geom_vline(xintercept=trueLevel, color="red") + facet_wrap(~type) + theme_bw()

# the mean estimate of all late terminations

res %>% group_by(type) %>% filter(n==50) %>% summarise(lateEst = mean(empMean))

# how many strategy A studies were significant?

res %>% filter(type=="sequential") %>% .[["selected"]] %>% table()

# Compute estimated lead levels

est.noBias <- res %>% group_by(type) %>% dplyr::summarise(

bias = FALSE,

naive.mean = mean(empMean),

weighted.mean = weighted.mean(empMean, w=n)

)

est.Bias <- res %>% filter(selected==TRUE) %>% group_by(type) %>% dplyr::summarise(

bias = TRUE,

naive.mean = mean(empMean),

weighted.mean = weighted.mean(empMean, w=n)

)

est <- rbind(est.noBias, est.Bias)

est

# output a html table

est.display <- txtRound(data.frame(est), 2, excl.cols=1:2)

t1 <- htmlTable(est.display,

header = c("Sampling plan", "PubBias", "Naive mean", "Weighted mean"),

rnames = FALSE)

t1

cat(t1)

library(dplyr)

library(htmlTable)

# set seed for reproducibility

set.seed(0xBEEF)

trueLevel <- 3

trueSD <- 2

safeLevel <- 2.7

maxN <- 50

minN <- 3

B <- 10000 # number of Monte Carlo simulations

res <- data.frame()

for (i in 1:B) {

print(paste0(i, "/", B))

maxSample <- rnorm(maxN, trueLevel, trueSD)

# optional stopping

for (n in minN:maxN) {

t0 <- t.test(maxSample[1:n], mu=safeLevel, alternative="less")

#print(paste0("n=", n, "; ", t0$estimate, ": ", t0$p.value))

if (t0$p.value <= .05) break;

}

finalSample.seq <- maxSample[1:n]

# now construct a matched fixed-n

finalSample.fixed <- rnorm(n, trueLevel, trueSD)

# ---------------------------------------------------------------------

# save results in long format

# sequential design

res <- rbind(res, data.frame(

id = i,

type = "sequential",

n = n,

p.value = t0$p.value,

selected = t0$p.value <= .05,

empMean = mean(finalSample.seq)

))

# fixed design

res <- rbind(res, data.frame(

id = i,

type = "fixed",

n = n,

p.value = NA,

selected = mean(finalSample.fixed) <= safeLevel, # some arbitrary publication bias selection

empMean = mean(finalSample.fixed)

))

}

save(res, file="res.RData")

# load("res.RData")

# Figure 1: Sampling distribution

ggplot(res, aes(x=n, y=empMean)) + geom_jitter(height=0, alpha=0.15) + xlab("Sample size") + ylab("Sample mean") + geom_hline(yintercept=trueLevel, color="red") + facet_wrap(~type) + theme_bw()

# Figure 2: Individual study estimates

ggplot(res, aes(x=empMean)) + geom_density() + xlab("Sample mean") + geom_vline(xintercept=trueLevel, color="red") + facet_wrap(~type) + theme_bw()

# the mean estimate of all late terminations

res %>% group_by(type) %>% filter(n==50) %>% summarise(lateEst = mean(empMean))

# how many strategy A studies were significant?

res %>% filter(type=="sequential") %>% .[["selected"]] %>% table()

# Compute estimated lead levels

est.noBias <- res %>% group_by(type) %>% dplyr::summarise(

bias = FALSE,

naive.mean = mean(empMean),

weighted.mean = weighted.mean(empMean, w=n)

)

est.Bias <- res %>% filter(selected==TRUE) %>% group_by(type) %>% dplyr::summarise(

bias = TRUE,

naive.mean = mean(empMean),

weighted.mean = weighted.mean(empMean, w=n)

)

est <- rbind(est.noBias, est.Bias)

est

# output a html table

est.display <- txtRound(data.frame(est), 2, excl.cols=1:2)

t1 <- htmlTable(est.display,

header = c("Sampling plan", "PubBias", "Naive mean", "Weighted mean"),

rnames = FALSE)

t1

cat(t1)

February 16, 2016 Filed in: Open Science | Science

The Psychology Department at LMU Munich continues to** change the incentive structure towards reproducible and open science**. The internal distribution of funding now partly is based on transparency criteria: Publications with open data, open material and pre-registrations get bonus points which directly translate into larger money allocations for that research unit.

Related posts: