---
title             : "Guidelines for developing reproducible workflows"
shorttitle        : "Reproducible workflows"

author: 
  - name          : "Shravan Vasishth"
    affiliation   : "1"
    corresponding : yes    # Define only one corresponding author
    address       : "University of Potsdam"
    email         : "vasishth@uni-potsdam.de"

affiliation:
  - id            : "1"
    institution   : "University of Potsdam"

authornote: |
    Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project number 317633480 – SFB 1287, project Q.
 
abstract: |
  These are some suggested guidelines on developing a reproducible workflow for carrying out research in any experiimental domain. The guidelines are not meant to be hard and fast rules, but rather are intended to be suggestions. The ultimate aim of these guidelines is to allow the researcher planning a study to have some informed basis for deciding what sample size to choose, to allow anyone in the field to examine and reanalyze your data and code, and to allow anyone to use your research results for carrying out further analyses, such as meta-analyses. 
  
  <!-- https://tinyurl.com/ybremelq -->
  
keywords          : "data analysis;  open science; transparency; power analysis"
wordcount         : "5352 words"

bibliography      : ["bibliography.bib"]

floatsintext      : no
figurelist        : no
tablelist         : no
footnotelist      : no
linenumbers       : no
mask              : no
draft             : no
toc               : true
documentclass     : "apa6"
classoption       : "doc"
output            : papaja::apa6_pdf
header-includes: |
  \usepackage{gb4e}\noautomath
  \usepackage{todonotes}
  \usepackage[utf8]{inputenc}
  \usepackage{fancyvrb}
---

```{r setup, include = FALSE}
library("papaja")
```


```{r analysis-preferences}
# Seed for random number generation
set.seed(42)
knitr::opts_chunk$set(cache.extra = knitr::rand_seed,
                      echo=TRUE)
library(brms)
library(lme4)
library(tidyverse)
library(ggplot2)
library(lattice)
```

# Revision history and online repository

- First version: 8 April 2017 following comments from Nieuwland and others.
- Revision: 10 Feb 2020. Added power analyses.

All code and data are available from: 

https://vasishth.github.io/ReproducibleWorkflows/

# The motivation for these guidelines

As a PhD advisor, I often encounter the following situation: the student carries out a study, and presents their experimental results to me. The code is often not well organized, and almost never documented. Students writing their first paper will often hard-code numerical values such as model estimates in their papers; inevitably, what happens is that something in the data or modeling changes as a result of the reviewing process or due to some other reason (coding errors are common), which then necessitates recopying all the code and data into the paper again. This is very error-prone, not to mention painfully inefficient. It is possible to automate the entire process, as I show below.

As a reviewer and action editor for journals, one problem I face regularly is that I review a paper, but I can't understand what statistical models the authors actually fit. Often I can't even figure out what the estimates of the effect were, because people usually don't even reveal basic information about the effect (such as the estimated coefficient and the standard error of the estimate; or what the parameter estimates of the different variance components were). 

As a researcher, I carry out a lot of meta-analyses. Whenever I ask a researcher for the data and code behind the paper, I am usually met with silence, or some reason why the data or code cannot be released. In my experience, only about 25% of researchers are able to release their data and code. **Not providing the data and code that lie behind a paper is a waste of research funds**. It leads to loss of information, and prevents us from incrementally building up our knowledge about a research problem. To see how valuable meta-analyses are, see some of our recent work: @NicenboimPreactivation2019, @JaegerEngelmannVasishth2017, @NicenboimRoettgeretal, @BuerkiEtAl2020. You can also watch a talk I recently gave that spells out the importance of meta-analysis in computational modeling:  https://youtu.be/UfB6JYaIY9I.

This tutorial aims to provide some suggestions on how to develop a reasonably reproducible workflow, which will allow future you and other researchers to revisit your work, and to build on it. As an example, I will use a replication attempt by Nieuwland et al of a famous EEG study from Marta Kutas' lab.

# Introduction: DeLong et al 2005 replication

Recently, @nieuwland2018large ([see here](https://elifesciences.org/articles/33468)) carried out an interesting series of replication attempts of the @DeLong2005 study ([see here](http://search.proquest.com/openview/f549a27edaad0cedad040e0df278f9ec/1?pq-origsite=gscholar&cbl=44706)) Nature Neuroscience paper. 

DeLong and colleagues examined the effect of predicting an upcoming noun at the determiner region that precedes the noun. Consider the sentences below:

\begin{exe}
\ex \label{ex:delong}
\begin{xlist}
\ex The day was breezy so the boy went outside to fly \textbf{\underline{a} kite}. \label{ex:delongP}
\ex The day was breezy so the boy went outside to fly \textbf{\underline{an} airplane}. \label{ex:delongU}
\end{xlist}
\end{exe}

Participants were shown sentences ending with a predictable noun phrase, such as *'a kite'* in (\@ref(ex:delongP)), or an unpredictable one, such as *'an airplane'* in (\@ref(ex:delongU)). 

In both examples, the determiners preceding the critical noun have the same meaning, so  there should  be no difference between the fit of *'a'* and *'an'* to the semantic context (i.e., both determiners should  incur the same integration costs). DeLong and colleagues showed that the amplitude of the negativity was smaller (becomes more positive) with increasing cloze probability at the determiner (and also at the noun; but that effect is expected given prior work, so not so surprising). 

Nieuwland et al report a failure to replicate the original effect. Here, we focus only on the article data from the Nieuwland et al study. 

Incidentally, DeLong et al never released their data, even though they have it available in a form that could be made available publicly. This sort of situation should never happen; data and code should be made available routinely with a published paper. I believe that leading journals like the *Journal of Memory and Language* now require data and code release, unless there is a compelling reason (e.g., privacy issues)  not to make the data public.

Nieuwland et al's work is exemplary in this respect; they  publicly released their data and code even before the paper was accepted. First, let's look at how Nieuwland et al made their data and code available. 

## Data and code from Nieuwland et al

Here are Nieuwland et al's data and code. 

- All materials were made available on the Open Science Foundation (OSF) repository: https://osf.io/eyzaq/
- The tree structure of the directory was as follows when I downloaded it in 2017:

```
.
|- 500ms_baseline
      |- art_b5_P.txt
      |- art_b5_R.txt
      |- noun_b5_P.txt
      |- noun_b5_R.txt
|- BFs
      |- Art_Bayes_Factor_Replication_Delong.txt
      |- Nouns_Bayes_Factor_Replication_Delong.txt
      |- orig_art.txt
      |- orig_noun.txt
      |- rep_art.txt
      |- rep_noun.txt
|- Nieuwland_etal_elife_accepted.pdf
|- ReplicationFunctionsCorrelation.R
|- artfinalP.txt
|- artfinalR.txt
|- cor_data_art.txt
|- cor_data_noun.txt
|- nounfinalP.txt
|- nounfinalR.txt
|- public_article_data.txt
|- trial_level_data
    |- art_replication_analysis_500msbaseline.txt
    |- art_replication_analysis_original.txt
    |- nouns_replication_analysis_500msbaseline.txt
    |- nouns_replication_analysis_original.txt
```

**What is missing here is a README file. Many of the filenames are essentially self-explanatory, but not all.**

## The R code file

The R code is presented as a file called public_script.txt.

Some suggestions for improvement:

- An R markdown file would have allowed the reader to see the precomputed output. One can always extract the R code from an R Markdown file called, say, yourfile.Rmd, by doing:

```
library(knitr)
purl("yourfile.Rmd")
```

- Running SessionInfo() provides some minimal information on the packages used and their versions; this can help in diagnosing problems if the code fails to run due to lack of backward compatibility.
- The call to the file has a hard-coded absolute path. This will fail to run on your computer:

```
#read article data
articles <-read.delim("~/Desktop/DELONGR/public_article_data.txt", quote="")
```

**Suggestion**:  give relative path names (relative to current working directory).

```{r loaddata}
articles <-read.delim("NieuwlandEtAl2018/public_article_data.txt", 
                      quote="")
head(articles)
```

- Avoid using T and F for TRUE and FALSE. This can play havoc with your code if T and F are bound variables. RStudio allows auto-completion, so one can just type  T and hit tab, and one should get TRUE (same for F and FALSE).

```
clozemean <- mean( articles$cloze, na.rm=T )
clozesd <- sd( articles$cloze, na.rm=T )
```

-  I could not find a column called n400, so this line doesn't run:

```
articles$base100 <- as.numeric(as.character(articles$n400))
```


Here is the code that I had in 2017 (it may have been updated since by Nieuwland et al):


```{r originalcode}
# R script for the DeLong replication data.
# base100 is the N400 dependent measure with the pre-registered 100 ms baseline, which will be used throughout this script to keep it short
# base200 has the 200 ms baseline, base500 the 500 ms baseline
# filt01 has the 0.01 hz filtered data, and pre-stim the -500 to -100 ms data 
#read article data: loaded above by SV.
#articles <-read.delim("~/Desktop/DELONGR/public_article_data.txt", quote="")

# make sure data is right format
articles$item <-as.factor(articles$item)
articles$cloze <- as.numeric(as.character(articles$cloze))
## no column called n400:
##articles$base100 <- ##as.numeric(as.character(articles$n400))

# Z-transform cloze and store the mean and SD
articles$zcloze <- scale(articles$cloze, center = TRUE, scale = TRUE)
clozemean <- mean( articles$cloze, na.rm=T )
clozesd <- sd( articles$cloze, na.rm=T )
```

Four linear mixed models are fit for the determiner. Most of the models fail to converge for me.

```{r articleanalysis,cache=TRUE}
# article models
model1a <- lmer(base100 ~ lab*zcloze + 
                  ( zcloze | subject) + (zcloze  | item),
                contrasts=list( lab=contr.sum(9) ), data = articles , 
                control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )

model2a <- lmer(base100 ~ lab/zcloze + 
                  ( zcloze | subject) + (zcloze  | item), 
                contrasts=list( lab=contr.sum(9) ), data = articles , 
                control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )

model3a <- lmer(base100 ~  zcloze + ( zcloze | subject) +
                  (zcloze  | item),  
                data = articles,
                control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )

model4a <- lmer(base100 ~         ( zcloze | subject) +
                  (zcloze  | item),  data = articles ,
                control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )

# compare models
anova(model1a,model2a)
```


Linear mixed models for the noun:

```{r nounanalysis,cache=TRUE,eval=FALSE}
########################
#### NOUNS

#nouns <-read.delim("~/Desktop/DELONGR/public_noun_data.txt", quote="")
nouns <-read.delim("NieuwlandEtAl2018/public_noun_data.txt", quote="")

nouns$item <-as.factor(nouns$item)
nouns$cloze <- as.numeric(as.character(nouns$cloze))
nouns$n400 <- as.numeric(as.character(nouns$n400))

nouns$zcloze <- scale(nouns$cloze, center = TRUE, scale = TRUE)
clozemean <- mean( nouns$cloze, na.rm=T )
clozesd <- sd( nouns$cloze, na.rm=T )

model1n <- lmer(n400 ~  lab*zcloze + 
                  ( zcloze | subject) + ( zcloze  | item), contrasts=list( lab=contr.sum(9) ),
                data = nouns , 
                control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )
model2n <- lmer(n400 ~  lab + zcloze + 
                  (zcloze | subject) + (zcloze  | item), contrasts=list( lab=contr.sum(9) ), 
                data = nouns , 
                control=lmerControl(optCtrl=list(maxfun=1e9)), REML = FALSE )
model3n <- lmer(n400 ~  zcloze + 
                  (zcloze | subject) + (zcloze  | item), data = nouns , 
                control=lmerControl(optCtrl=list(maxfun=1e9)), REML = FALSE )
model4n <- lmer(n400 ~   
                  (zcloze | subject) + (zcloze  | item), data = nouns , 
                control=lmerControl(optCtrl=list(maxfun=1e9)), REML = FALSE )

# compare models
anova(model1n,model2n)

```

## Results for the determiner

Let's focus on just one of the models:

```{r model3a}
model3a <- lmer(base100 ~  zcloze + 
                  ( zcloze | subject) + (zcloze  | item),  
                data = articles, 
                control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )
model3aNULL <- lmer(base100 ~  1 + 
                      ( zcloze | subject) + (zcloze  | item),  
                    data = articles, 
                    control=lmerControl(optCtrl=list(maxfun=1e5)), REML = FALSE )
anova(model3a,model3aNULL)
```

There doesn't seem to be evidence for zcloze affecting the response.

One useful thing to look at is the u0 and u1 values for each participant and item:

```{r dotplotmodel3a,fig.height=10,fig.width=8}
print(dotplot(ranef(model3a,condVar=TRUE)))
```

There are more sophisticated ways to plot individual level effects. See our Bayesian modeling lecture notes [here](https://vasishth.github.io/Bayes_CogSci/);  see section [5.1.3](https://vasishth.github.io/Bayes_CogSci/a-hierarchical-normal-model-the-n400-effect.html#sec:uncorrelated).

# Planning a future study based on existing data

This assumes a frequentist analysis.

We will:

- begin with the existing data from Nieuwland et al 2018
- first visualize the data
- then use all available data to compute an estimate of the effect (meta-analysis)
- compute sample size needed to achieve 80% power using simulation
- carry out new study
- an alternative Bayesian approach would be to use the estimate of the effect from the previous studies as an informative prior in a Bayesian analysis.

## Simulating data is a very important skill

One skill we will learn here is how to simulate data that reflects the underlying generative process we are assuming from our experiment. Simulating data is a very important skill for the data analyst because:

- one can determine whether the model we have specified for data analysis is sensible (generates reasonable data)
- one can determine whether the model we have specified can in principle recover the parameters of interest under repeated sampling
- one can easily do power calculations (although this can take some time to compute for complex models)

## Case study: The Niewland et al study data

### Subset the existing data

First, choose the relevant columns:

```{r subsetdata}
dat<-dplyr::select(articles,subject,item,lab,zcloze,base100)
```

If we are using  all the data, it might be nicer if we re-name all the subjects by numerical id's instead of lab id plus the  subject  id.

```{r renamesubjects}
## get subject names
subjnames<-unique(dat$subject)
## convert to numerical values:
subjects_numerical<-data.frame(subjname=subjnames,subj=factor(1:length(subjnames)))

## create new column with numerical id's:
dat2<-merge(dat,subjects_numerical,by.x="subject",by.y="subjname")
dat<-dat2
head(dat)
```

### Visualize the data by subject and by item

There is data from 9 labs. We can look at the by-subjects data for each lab. For example, here is the data from Birmingham.

```{r bysubject,fig.width=7,fig.height=7}
## save lab names:
labnames<-as.character(unique(dat$lab))

## Choose first lab's data:
current_lab<-filter(dat,lab==labnames[1])

xyplot(base100~zcloze|subj,current_lab,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.abline(lm(y ~ x))
       },main="By subjects data")
```

```{r byitem,fig.width=7,fig.height=7}
xyplot(base100~zcloze|item,
       current_lab,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.abline(lm(y ~ x))
       },,main="By items data")
```

One can do this in ggplot as well. First, we define a function that mimics the lattice plot. I got this code from someone on twitter but I lost the source.

```{r definexyplotfunction}
gg_xyplot <- function(x, y, formula, shape, size, 
                      xlabel="zcloze", 
                      ylabel="base100 (microvolts)",
                      data=current_lab){
    ggplot(data = data, aes(x = data[,x],
                            y = data[,y]))  +
    facet_wrap(formula) +
    geom_smooth(method="lm")+
    geom_point(color = "blue", 
               shape = shape, size = size) +
    theme(panel.grid.minor = element_blank()) +
    theme_bw() +
    ylab(ylabel) +
    xlab(xlabel)
}
```

```{r ggxyplotsubj,fig.width=7,fig.height=7}
gg_xyplot(x = "zcloze", y = "base100",  ~ subj,  
          shape = 1, size = 3, 
          data = current_lab)
          
```

```{r ggxyplotitem,fig.width=7,fig.height=7}
gg_xyplot(x = "zcloze", y = "base100",  ~ item,  
          shape = 1, size = 3, 
          data = current_lab)
          
```

One problem with the above visualizations is that we would have to plot nine labs' plots separately and not get an overall impression of the by-subject variability in the zcloze effect. 

```{r lmList}
## fit separate linear models for each subject: 
m_lmlist<-lmList(base100~zcloze|subj,dat)
lmlist_coef<-summary(m_lmlist)$coefficients
## extract by subject slopes:
slopes<-lmlist_coef[,,2]
means<-slopes[,1]
lower<-means-2*slopes[,2]
upper<-means+2*slopes[,2]

subjects<-unique(dat$subj)

slopes_summary<-data.frame(subj=subjects,
                       means=means,
                       lower=lower,upper=upper)

## reorder means by magnitude:
slopes_summary<-slopes_summary[order(slopes_summary$means),]

p<-ggplot(slopes_summary,aes(x=subjects, y=means)) + 
  geom_line() +
  geom_point()+
  geom_errorbar(aes(ymin=lower, ymax=upper), width=.2,
                 position=position_dodge(0.05))+theme_bw()
p
```

Notice that there is a lot of variability across the subjects (each subject's data is being analyzed separately). Some show negative effects, some show the expected positive effects.  This kind of oscillation is characteristic of low power experiments; the overall experiment by Nieuwland et al might not be low power, but analyzing each subject separately definitely is going to lead to an underpowered design! See @gelmancarlin for discussion on how this kind of oscillation of the effects around 0 happens.

### Power analyses

In order to understand the simulation approach I take here, it is important to understand how linear mixed models work. Here, we are going to switch to the slides [accompanying this file, on linear mixed models introduction](https://github.com/vasishth/ReproducibleWorkflows/blob/master/LinearMixedModels/LinearMixedModels.pdf).

Incidentally, in order to avoid messages about convergence failures, I am going to add the following to the lmer function:

```
control=lmerControl(calc.derivs=FALSE)
```

This is just for convenience; normally, we will have to worry about convergence problems. See @barr2013, @hannesBEAP, @BatesEtAlParsimonious.

#### Computing power distribution for the current sample size

##### Step 1: Extract estimates

First, we extract the estimated parameter values from the model  that was fit to the whole data-set:

```{r extractestimates}
## extract estimates of fixed-effects parameters:
beta<-summary(model3a)$coefficients[,1]
## extract standard deviation estimate:
sigma_e<-attr(VarCorr(model3a),"sc")
## assemble variance covariance matrix for subjects:
subj_ranefsd<-attr(VarCorr(model3a)$subj,"stddev")
subj_ranefcorr<-attr(VarCorr(model3a)$subj,"corr")
Sigma_u<-round(SIN::sdcor2cov(stddev=subj_ranefsd,
                        corr=subj_ranefcorr),6)

## assemble variance covariance matrix for items:
item_ranefsd<-attr(VarCorr(model3a)$item,"stddev")
item_ranefcorr<-attr(VarCorr(model3a)$item,"corr")
Sigma_w<-SIN::sdcor2cov(stddev=item_ranefsd,
                        corr=item_ranefcorr)
```

##### Step 2: Generate simulated data using estimates, compute power (or Type I error):

Load a function to generate simulated data. I explain the function in the appendix.

```{r generatefakedata}
source("R/gen_sim_norm.R")
```

```{r producesimulatedvalues,eval=FALSE}
## Example simulated data, generated several times:
nsim<-10
sim_values<-matrix(rep(NA,nsim*dim(dat)[1]),ncol = nsim)
for(i in 1:nsim){
simdat<-gen_sim_norm(dat,
               alpha=beta[1],beta=beta[2],
               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
               sigma_e=sigma_e)
sim_values[,i]<-simdat$simbase100
}
head(sim_values)
```

It is useful to compare real and simulated data, to check if the generative model produces  realistic data.

```{r comparerealsim,fig.height=5,fig.weight=5}
## Compare fake and simulated data:
## load precomputed simulated values:
load("data/sim_values.Rda")
## the observed data:
hist(dat$base100,freq=FALSE,
     ylim=c(0,0.06),
     main="Real vs simulated data",
     xlab="base100")
## simulated data:
for(i in 1:ncol(sim_values)){
lines(density(sim_values[,i]),lty=2)
}
```

Here, it looks like the model is producing realistic data.

A critical question here is: what should the estimate of the effect be? Normally, for the problems I study, the estimate would be derived from a computational model or from a meta-analysis of existing data [@NicenboimPreactivation2019, @JaegerEngelmannVasishth2017, @NicenboimRoettgeretal, @BuerkiEtAl2020]. For now, we will take the estimate (mean and SE) from the data; this is just a convenient starting point. We will improve on this later in an exercise.

```{r}
b1<-round(summary(model3a)$coefficients[2,1],3)
b1se<-round(summary(model3a)$coefficients[2,2],3)
```

Given these estimates, we can compute the lower and upper bounds of the estimate:

```{r}
b1-2*b1se; b1+2*b1se
```

We could repeatedly (100 times) simulate estimates from the distribution Normal(`r round(b1,2)`,`r round(b1se,2)`). Each time that we take a sample, we generate simulated data 100 times, and then compute the proportion of times that the null hypothesis is rejected. This gives us 100 estimates of power; we can plot these estimates as a distribution.

But the above approach will take a lot of time, so for this demonstration, I would compute the power estimate based only on the estimated mean effect. One could also compute power based on the lower bound and upper bound of the estimated effect (i.e., mean $\pm$ 2$\times$SE).  

```{r poweranalysismean,cache=TRUE}
## store for power:
nsim<-100
tvals<-c()
for(i in 1:nsim){
  #print(paste("i=",i,sep=""))
simdat<-gen_sim_norm(dat,
               alpha=beta[1],
               beta=b1, ## using the estimated slope
               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
               sigma_e=sigma_e)
m<-lmer(simbase100~zcloze+(1+zcloze|subj)+
          (1+zcloze|item),simdat,
        control=lmerControl(calc.derivs=FALSE))
tvals[i]<-summary(m)$coefficients[2,3]
}
power_mean<-mean(abs(tvals)>2)
power_mean
save(power_mean,file="data/power_mean.Rda")
```

Power for the estimated mean from the existing studies is `r round(power_mean,2)`.

I provide the code for computing power using the lower and upper bounds of the mean estimate, but I don't run the code here (by setting eval to FALSE), to save time. 

```{r poweranalysislower,cache=TRUE,eval=FALSE}
## store for power:
nsim<-100
tvals<-c()
for(i in 1:nsim){
  #print(paste("i=",i,sep=""))
simdat<-gen_sim_norm(dat,
               alpha=beta[1],
               beta=b1-2*b1se, ## using the lower bound
               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
               sigma_e=sigma_e)
m<-lmer(simbase100~zcloze+(1+zcloze|subj)+(1+zcloze|item),
        simdat,
        control=lmerControl(calc.derivs=FALSE))
tvals[i]<-summary(m)$coefficients[2,3]
}
power_lower<-mean(abs(tvals)>2)
power_lower
```

```{r poweranalysisupper,cache=TRUE,eval=FALSE}
## store for power:
nsim<-100
tvals<-c()
for(i in 1:nsim){
  #print(paste("i=",i,sep=""))
simdat<-gen_sim_norm(dat,
               alpha=beta[1],
               beta=b1+2*b1se, ## using the upper bound
               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
               sigma_e=sigma_e)
m<-lmer(simbase100~zcloze+(1+zcloze|subj)+(1+zcloze|item),
        simdat,
        control=lmerControl(calc.derivs=FALSE))
tvals[i]<-summary(m)$coefficients[2,3]
}
power_upper<-mean(abs(tvals)>2)
power_upper
```

In summary, the power estimate for the mean is `r round(power_mean,2)`. If you run the estimates for the lower and upper bound, you will find that there is a lot  of  uncertainty on the power estimate; this uncertainty is coming from the fact that  we are unsure about the magnitude of the true effect.

<style>
div.blue { background-color:#e6f0ff; border-radius: 5px; padding: 20px;}
</style>
<div class = "blue">
#### Exercise

Write a function computepower that can produce a power estimate given a particular effect size. You should be able to write something like

computepower(b=0.10)

and it  should return the power estimate for the effect size  estimate 0.10 microvolts.

</div>


#### Computing sample size needed for 80% power

Next, we want to know how many subjects we will need to obtain 80% power. To do this, we will  replicate our existing data frame n times, where n can 2, 3, 4, etc.

For illustration, suppose our data frame is only the Birmingham data. 

We can use the data.table package to quickly replicate the data frame. Here is an example of how 

```{r datareplicationexample}
library(data.table)
df <- data.frame(a = c(1,2,3), b = c(1,2,3))
dt <- as.data.table(df)
n <- 3
dt[rep(dt[, .I], n)]
```

Select the relevant columns:

```{r selectbirmdata}
birm_dat<-dplyr::select(current_lab,
                        subj,item,zcloze)
head(birm_dat)
```

Generate 3 replicates of the data frame for generating a larger set of simulated data. Note that one has to update the subject ids.

```{r replicatebirmdata}
birm_dat<-as.data.table(birm_dat)
subjid<-birm_dat$subj
n <- 3
repl_dat<-birm_dat[rep(birm_dat[, .I], n)]
dim(birm_dat)[1]*3
dim(repl_dat)[1]

## create new subject vector extended n times
add_id<-rep(seq(100,100*n,by=100),each=dim(birm_dat)[1])
repl_dat$subj<-rep(as.numeric(as.character(birm_dat$subj)),n)+add_id
head(repl_dat)
length(unique(birm_dat$subj))*3
length(unique(repl_dat$subj))
```

Then we simulate data to compute power for a given sample size:

```{r computesamplesize,cache=TRUE}
nsim<-100
tvals<-c()
for(i in 1:nsim){
simdat<-gen_sim_norm(repl_dat,
               alpha=beta[1],
               beta=b1, ## using the estimated slope
               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
               sigma_e=sigma_e)
m<-lmer(simbase100~zcloze+(1+zcloze|subj)+(1+zcloze|item),simdat,
        control=lmerControl(calc.derivs=FALSE))
tvals[i]<-summary(m)$coefficients[2,3]
}
## number of subjects:
length(unique(simdat$subj))
power_mean<-mean(abs(tvals)>2)
power_mean
```

Obviously, 126 subjects is too  little here! We will need a *lot* more to achieve 80% power!  One can now re-run the above code with increasing numbers of subjects to get a power curve as a function of  the subject sample size.

We can write a function to simplify the process of computing power as a function of subject sample size.

```{r computepowerfunction}
compute_power<-function(dat=birm_dat,
                        replicates=3,nsims=100){
dat<-as.data.table(dat)
subjid<-dat$subj
nrep <- replicates
repl_dat<-dat[rep(dat[, .I], nrep)]

## create new subject vector extended n times
add_id<-rep(seq(100,100*nrep,by=100),each=dim(dat)[1])
repl_dat$subj<-rep(as.numeric(as.character(birm_dat$subj)),
                   nrep)+add_id
nsim<-nsims
tvals<-c()
for(i in 1:nsim){
simdat<-gen_sim_norm(repl_dat,
               alpha=beta[1],
               beta=b1, ## using the estimated slope
               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
               sigma_e=sigma_e)
m<-lmer(simbase100~zcloze+(1+zcloze|subj)+(1+zcloze|item),
        simdat,
        control=lmerControl(calc.derivs=FALSE))
tvals[i]<-summary(m)$coefficients[2,3]
}
power_mean<-mean(abs(tvals)>2)
## return result with sample size:
res<-c(length(unique(simdat$subj)),
                power_mean)
res
}
```

```{r computesamplesize2,cache=TRUE}
repl4<-compute_power(dat=birm_dat,
              replicates = 4,
              nsims=10)
repl4
repl6<-compute_power(dat=birm_dat,
              replicates = 6,
              nsims=10)
repl6
```

Now, with the above code, in principle we can display the effect of increasing subject sample size on power. Here is a simplified version of such a table:

```{r powersamplesizeplot}
results<-rbind(repl4,repl6)
results<-data.frame(results)
colnames(results)<-c("nsubj","power estimate")
results
```

Don't be surprised if the power doesn't go up much with sample size, or if it even goes down with increasing sample size. Because power is being computed using simulated data, we will see some variation in the power calculation from one simulation to another.

<style>
div.blue { background-color:#e6f0ff; border-radius: 5px; padding: 20px;}
</style>
<div class = "blue">
**Exercise (to be  done after the workshop)**

This exercise should be done overnight! Compute power for 2,4,6,8,10  replicates of the Birmingham data, (use 100 simulations, not 10) and plot power against sample size.

</div>

## Some comments on power

Note here that power computed from the data you already have is called **post-hoc power**, and is pointless to compute; this is because once the p-value is known, power can be computed analytically [@hoenigheisey].  In other words, post-hoc power is just a transform of the observed p-value.

What we are doing here is using existing data to estimate **variance  components**. Then, we use simulation to compute a power *distribution* for a range of plausible values of the effect of interest. 

The key problem here is figuring out what a plausible range of values of the effect is. Here, we have simply used all the available data to get an estimate of the effect. This makes our power analysis look like a post-hoc power calculation. But I used the observed estimates of the effect here just as an example.

In realistic data, a better way to obtain estimates of the effect is by doing a meta-analysis. See @NicenboimPreactivation2019 for an estimate based on a meta-analysis. The estimate based on all publicly available data is: 0.11 microvolts, 95% credible intervals [0.05,0.16]. One can use that to compute power. 

<style>
div.blue { background-color:#e6f0ff; border-radius: 5px; padding: 20px;}
</style>
<div class = "blue">
### Exercise

Use the existing data and the above meta-analysis estimates to work out estimates of power for a sample size of 334 subjects. Use the mean, lower, and upper bounds of the meta-analysis estimate to work out the power estimates.  

</div>

# Suggestions for making data and code public

- Develop and release a data management plan (required by the DFG). Example: https://osf.io/f9dqk/
- Use R Markdown to produce commented and compiled code output so that the reader doesn't need to re-run code.
- Once the analysis is complete, put all the code and data on OSF or github (or link github to osf). Example: https://github.com/vasishth/StatSigFilter
- Make sure you don't depend on data and functions loaded via a hidden .Rprofile file that is local to your computer.
- Avoid Excel files for storing data that will be loaded into R; save data in .csv or .txt files.
- The directory structure for the data release should have a structure that the outsider can easily work out where everything is. The structure can be based on R package structures (discussed below), but doesn't have to be.
- Create a README file that users can look at for guidance. Sometimes names of columns in data-frames are not easy to interpret. E.g., in eyetracking data, the user cnanot know  whether RP is regression probability or regression path duration.  
- At the end of the R Markdown file, provide session information using SessionInfo().
- Make sure that your R Markdown file runs! Ideally, download all the code and data onto another machine and check that it runs.
-  If you have very large files that can't be saved on github, osf, etc., then store them on a (university) server and link to the files. Example: https://osf.io/reavs/
- Fitting complicated models can lead to problems when compiling an Rmd file. One way around that problem is to save the results (or a summary) of the resulting model offline and to load it when needed. 
- One irritating problem with R Markdown is that the error messages on compilation failure are often not very informative. One solution to this is to incrementally build the file, using separate ``child'' files. The R Markdown online documentation explains how to do this. Despite the frustration that R Markdown causes, the long-term advantages of using it are great.

# Appendix 1  

## Explanation of the function in gen_sim_norm.R

The way I generate data is by traversing the observed data's data frame row by row, and then generating simulated data for each row, using the following pieces of information:

- the subject id
- the item id
- the parameter estimates from the model that was fit on the whole data-set

This section presuppose knowledge of how a linear mixed model is ``assembled''. The underlying model is as follows.

i indexes subjects, j items.

\begin{equation}
y_{ij} = \alpha + u_{0i} + w_{0j} + (\beta + u_{1i} + w_{1j}) * zcloze_{ij} + \varepsilon_{ij}
\end{equation}

where $\varepsilon_{ij} \sim Normal(0,\sigma)$ and 

\begin{equation}\label{eq:covmat2}
\Sigma _u
=
\begin{pmatrix}
\sigma _{u0}^2  & \rho _{u}\sigma _{u0}\sigma _{u1}\\
\rho _{u}\sigma _{u0}\sigma _{u1}    & \sigma _{u1}^2\\
\end{pmatrix}
\quad 
\Sigma _w
=
\begin{pmatrix}
\sigma _{w0}^2  & \rho _{w}\sigma _{w0}\sigma _{w1}\\
\rho _{w}\sigma _{w0}\sigma _{w1}    & \sigma _{w1}^2\\
\end{pmatrix}
\end{equation}

\begin{equation}\label{eq:jointpriordist2}
\begin{pmatrix}
  u_0 \\ 
  u_1 \\
\end{pmatrix}
\sim 
\mathcal{N} \left(
\begin{pmatrix}
  0 \\
  0 \\
\end{pmatrix},
\Sigma_{u}
\right),
\quad
\begin{pmatrix}
  w_0 \\ 
  w_1 \\
\end{pmatrix}
\sim 
\mathcal{N}\left(
\begin{pmatrix}
  0 \\
  0 \\
\end{pmatrix},
\Sigma_{w}
\right)
\end{equation}

Given such a mathematical model, we can generate simulated data row by row, given subject and item ids, and given point value estimates for 

- $\alpha$
- $\beta$
- $\sigma$
- $\sigma_{u0},\sigma_{u1},\sigma_{w0},\sigma_{w1}$
- $\rho_u, \rho_w$

The function first creates subject and item random effects:

```{r generateranefs}
nsubj<-length(unique(dat$subj))
nitem<-length(unique(dat$item))
u<-mvrnorm(n=nsubj, # number of subjects
             mu=c(0,0),Sigma=Sigma_u)
w<-mvrnorm(n=nitem, # number of items
             mu=c(0,0),Sigma=Sigma_w)

## add subject and item id's:
u<-data.frame(subjid=unique(dat$subj),u)
w<-data.frame(itemid=unique(dat$item),w)

```

Each row comtains randomly generated intercept and slope adjustments for each of the subjects (items):

```{r}
head(u)
head(w)
```

If we want to know what subject 1's intercept adjustment is, we write:

```{r}
u[1,2]
```

If we want to know what subject 1's slope adjustment is, we write:

```{r}
u[1,3]
```

Now, for the first row of the  data frame ```dat```: 

```{r}
dat[1,]
```
we can produce randomly generated data as follows. In this row, we have subject id 1, and item id 102.  This can be extracted from the data frames u and w.

```{r}
## intercept adjustment:
u[u$subjid==1,2]
## slope adjustment:
u[u$subjid==1,3]
```

These simply need to be added to  the fixed effects intercept and slope values:

```{r}
beta[1] + u[u$subjid==1,2] + 
  w[w$itemid==102,2]+ 
  (beta[2] + 
     u[u$subjid==1,3] +
     w[w$itemid==102,3])*dat[1,]$zcloze +
  rnorm(1,0,sigma_e)
```

Now, if we repeat this for every row of the data-frame, we can produce simulated data. Obviously, we cannot write the subject and item id by hand for each row. But we can automatically extract it from the data frame ```dat```:

```{r}
dat[1,]$subj
dat[1,]$item
```

Thus, for each row, we can extract the subject and item id by indexing the row id with i ranging from 1 to the total number of rows N, and then plugging in that value into the equation below. I show what happens with i=1:

```{r}
i<-1
current_subjid<-dat[i,]$subj
current_itemid<-dat[i,]$item

beta[1] + u[u$subjid==current_subjid,2] +
  w[w$itemid==current_itemid,2]+ 
  (beta[2] + u[u$subjid==current_subjid,3] +
      w[w$itemid==current_itemid,3])*dat[i,]$zcloze + rnorm(1,0,sigma_e)
```

Now, if we change the index i to 2 (this is an exercise for you), we get the simulated value for row 2. 

One can now generalize this, and write a for-loop that traverses every row of the data frame:

```{r}
## number of rows
N<-dim(dat)[1]
# save vector of simulated data:
simbase100 <- rep(NA,N)
for(i in 1:N){
current_subjid<-dat[i,]$subj
current_itemid<-dat[i,]$item
simbase100[i] <-  beta[1] + u[u$subjid==current_subjid,2] +
  w[w$itemid==current_itemid,2]+ 
  (beta[2] + u[u$subjid==current_subjid,3] +
      w[w$itemid==current_itemid,3])*dat[i,]$zcloze + rnorm(1,0,sigma_e)
}
```

These data look pretty similar to the observed data:

```{r}
hist(dat$base100)
lines(density(simbase100))
```

## Why I don't use ready-made packages for computing power

There are ready-made packages in R that generate simulated data from a fitted lmer model. However, I never teach how to use those packages, because the procedure I show you above helps the researcher appreciate the underlying generative process that produced the data. It also allows for more flexible variants on the above approach. For example, if we think like a Bayesian and have a prior distribution on each parameter, we can sample from a prior distribution for the effect size or indeed for any of the variance components. This only takes a small modification of the above code and data. With a ready-made package, you only get what the package developer decided to allow you to do.

# Session information

```{r sessioninfo}
sessionInfo()
```

# References
```{r create_r-references,echo=FALSE}
r_refs(file = "bibliography.bib")
```

\begingroup
\setlength{\parindent}{-0.5in}
\setlength{\leftskip}{0.5in}

<div id = "bibliography"></div>
\endgroup
