09-multiplecomparisons.Rmd

# Understanding Type I error inflation using simulation {#ch:MultComp} 


```{r echo=FALSE}
library(knitr)
library(brms)
library(bcogsci)
library(papaja) 
library(MASS)
library(hypr)
library(lme4)
library(lingpsych)
library(ggplot2)
```


The simplest situation we encountered in the earliest chapters was to carry out a single hypothesis test, having set Type I error to 0.05. In the linear mixed model setting, even this simple scenario is fraught with danger; Type I error can become startlingly large depending on the kind of linear mixed model you fit. We begin this chapter by illustrating this problem through an example of a simple two-condition experiment. However, there are at least two other sources of Type I error inflation: multiple comparison, and model misspecification.  We illustrate these issues below with an example from a $2\times 2\times 2$ repeated measures design.

## Overly simple random effects structure in LMMs inflate Type I error

We return here to the data-set we saw earlier, the @grodner self-paced reading study (their Experiment 1) on English relative clauses.
First, load the data, fit a model with a full variance-covariance matrix for subjects and items, and then extract the parameters from the fitted model, just like we did in the simulation chapter.

```{r}
library(lingpsych)
## Grodner and Gibson 2005 Expt 1 data:
data("df_gg05e1")
gg05e1<-df_gg05e1
## look at the data frame:
head(gg05e1)
## convert raw RTs to log RTs:
gg05e1$logrt<-log(gg05e1$rawRT)
## sum-coding of predictor:
gg05e1$so<-ifelse(gg05e1$condition=="objgap",1/2,-1/2)
## fit model to log RTs:
m<-lmer(logrt ~ so + 
          (1+so|subject) + 
          (1+so|item), 
        data=gg05e1,
        ## "switch off" warnings:
        control=lmerControl(calc.derivs=FALSE))

## extract parameter estimates:
beta<-round(summary(m)$coefficients[,1],4)
sigma_e<-round(attr(VarCorr(m),"sc"),2)
subj_ranefsd<-round(attr(VarCorr(m)$subject,"stddev"),4)
subj_ranefcorr<-round(attr(VarCorr(m)$subject,"corr"),1)

## assemble variance-covariance matrix for subjects:
Sigma_u<-SIN::sdcor2cov(stddev=subj_ranefsd,
                        corr=subj_ranefcorr)

item_ranefsd<-round(attr(VarCorr(m)$item,"stddev"),4)
item_ranefcorr<-round(attr(VarCorr(m)$item,"corr"),1)

## assemble variance matrix for items:

## choose some intermediate values for correlations:
corr_matrix<-(diag(2) + matrix(rep(1,4),ncol=2))/2

Sigma_w<-SIN::sdcor2cov(stddev=item_ranefsd,
                        corr=corr_matrix)
```

We are now ready to compute Type I error using simulation.

### Type I error with a varying intercepts-only model

Compute Type I error for a varying intercepts model using log reading times:

```{r message=FALSE,warning=FALSE}
## null is true:
beta[2]<-0
nsim<-100
significant<-rep(NA,nsim)

for(i in 1:nsim){
#generate sim data:
simdat<-gen_sim_lnorm2(nitem=16,
                         nsubj=42,
                       beta=beta,
                       Sigma_u=Sigma_u,
                       Sigma_w=Sigma_w,
                      sigma_e=sigma_e)

## fit model to sim data:
m<-lmer(log(rt)~so+(1|subj)+(1|item),simdat,
        control=lmerControl(calc.derivs=FALSE))
## record whether t-value is significant:
significant[i]<-ifelse(abs(summary(m)$coefficients[2,3])>2,1,0)
}

## Type I error is quite high:
mean(significant)
```

We find that Type I error is quite high. This is the problem with fitting varying intercepts-only models that @barr2013 point out: if the data-generative process is such that varying slopes should have been included in the model, fitting a varying intercepts-only model generally will lead to Type I error inflation. Of course, the extent of the inflation will vary from data-set to data-set; no simulation can give us a general answer to the question of how much inflation will happen.

### Type I error with a varying intercepts and varying slopes model

Next, let's compute Type I error for a model with varying intercepts and varying slopes, but without correlations for the random effects (because these would be difficult to estimate with the number of subjects (42) and items (16) used in @grodner). The dependent measure is still log rt.

```{r warning=FALSE,message=FALSE,results="asis"}
significant<-rep(NA,nsim)

for(i in 1:nsim){
#generate sim data:
simdat<-gen_sim_lnorm2(nitem=16,
                         nsubj=42,
                       beta=beta,
                       Sigma_u=Sigma_u,
                       Sigma_w=Sigma_w,
                      sigma_e=sigma_e)

## fit model to sim data:
m<-lmer(log(rt)~so+(1+so|subj)+(1+so|item),simdat,
        control=lmerControl(calc.derivs=FALSE))
## is the t-value significant?
significant[i]<-ifelse(abs(summary(m)$coefficients[2,3])
                       >2,1,0)
}

## Type I error is the expected value:
mean(significant)
```

Now we find that Type I error is near the expected 5%. As @barr2013 pointed out in their paper, when varying slopes are included, Type I error will have the expected value. 

### Type I error inflation due to model mis-specification

Although this point is not discussed in @barr2013, an inappropriate likelihood function, like the Normal likelihood in the case where the true generative process involves a log-normal distribution, can also lead to inflated Type I error. This can happen even if the model has a full variance-covariance matrix specified for it. This is because of extreme values generated by the log-normal distribution being excessively influential in generating Type I errors.

```{r echo=FALSE,warning=FALSE,message=FALSE,results="asis"}
nsim<-100
intercept_log<-effect_estlog<-
  effect_estraw<-significantlog<-
  significantraw<-
  se_estlog<-se_estraw<-rep(NA,nsim)

for(i in 1:nsim){
#generate sim data:
simdat<-gen_sim_lnorm2(nitem=16,
                         nsubj=42,
                       beta=beta,
                       Sigma_u=Sigma_u,
                       Sigma_w=Sigma_w,
                      sigma_e=sigma_e)

## fit model to log rts:
mlog<-lmer(log(rt)~so+(1+so|subj)+(1+so|item),simdat,
        control=lmerControl(calc.derivs=FALSE))

## fit model to sim data but on raw RTs:
mraw<-lmer(rt~so+(1+so|subj)+(1+so|item),simdat,
        control=lmerControl(calc.derivs=FALSE))

## is the effect significant?
significantlog[i]<-ifelse(abs(summary(mlog)$coefficients[2,3])
                       >2,1,0)
significantraw[i]<-ifelse(abs(summary(mraw)$coefficients[2,3])
                       >2,1,0)

#store estimates:
intercept_log[i]<-summary(mlog)$coefficients[1,1]
effect_estlog[i]<-summary(mlog)$coefficients[2,1]
effect_estraw[i]<-summary(mraw)$coefficients[2,1]
se_estlog[i]<-summary(mlog)$coefficients[2,2]
se_estraw[i]<-summary(mraw)$coefficients[2,2]
}
```

Type I error is inflated even with a maximal model if the likelihood is mis-specified:

```{r}
mean(significantlog)
mean(significantraw)
```

One can also compare the estimates (along with their 95% confidence intervals) that were significant in the raw and log ms scales:

```{r echo=FALSE}
logsig<-which(significantlog==1)

## backtransform log estimates to raw ms:
b_log<-exp(intercept_log[logsig]+
                 effect_estlog[logsig]*0.5)-
  exp(intercept_log[logsig]-effect_estlog[logsig]*0.5)
b_log_upper<-exp(intercept_log[logsig]+
                 (effect_estlog[logsig]+2*se_estlog[logsig])*0.5)-
  exp(intercept_log[logsig]-(effect_estlog[logsig]+2*se_estlog[logsig])*0.5)
b_log_lower<-exp(intercept_log[logsig]+
                 (effect_estlog[logsig]-2*se_estlog[logsig])*0.5)-
  exp(intercept_log[logsig]-(effect_estlog[logsig]-2*se_estlog[logsig])*0.5)

rawsig<-which(significantraw==1)
b_raw<-effect_estraw[rawsig]

b_raw_upper<-b_raw+2*se_estraw[rawsig]
b_raw_lower<-b_raw-2*se_estraw[rawsig]

rawests<-data.frame(transform="raw",mn=b_raw,lower=b_raw_lower,upper=b_raw_upper)
logests<-data.frame(transform="log",mn=b_log,lower=b_log_lower,upper=b_log_upper)
ests<-rbind(rawests,logests)

n<-dim(ests)[1]
ests$n<-1:n

ests$width<-ests$upper-ests$lower
```

```{r}
pd <- position_dodge(.3)
ggplot(ests, aes(x=n, y=mn, colour=transform, group=transform)) + geom_errorbar(aes(ymin=lower, ymax=upper),
width=.2, size=0.25, colour="black", position=pd) + geom_line(position=pd) +
geom_point(position=pd, size=2.5)+theme_bw()+theme(axis.title.x = element_blank())

with(ests,tapply(width,transform,mean))
```

We see that the effect estimates that were significant were based on  over-/mis-estimates on the raw scale (the true effect is 0 ms).

## Type I error inflation due to multiple comparisons

Next, we consider a case where the design is more complex than a two-condition experiment. The data are from an eyetracking reading study conducted by @Dillon-EtAl-2013; this is a $2\times 2 \times 2$ repeated measures design. The dependent measure is total reading time at a particular word in the sentence. Because there are eight conditions, we can carry out at most seven comparisons. 

The question we want to answer here is: if we carry out multiple (here, seven) hypothesis tests, what will be the probability of rejecting at least one of them incorrectly? Analytically, the answer is easy to compute. Assume that the Type I error probability for each hypothesis test is 0.05. To keep things simple, suppose we are testing two null hypotheses.  The probability of getting both the two null hypothesis tests correctly accepted is $0.95^2$; this assumes that the tests are independent. It follows that the probability of getting at least one null hypothesis test incorrectly rejected is $1-0.95^2=0.0.0975$. That's higher than our desired Type I error of 0.05. 

\includegraphics{figures/probspace.png}

Next, let's see what simulation gives us when we fit a varying intercepts and slopes model for the case where we are testing seven hypotheses.

First, load the data, and obtain the estimates from a fitted linear mixed model with full variance-covariance matrices for subjects and items. Because fitting such a model with `lmer` will lead to convergence errors, we fit the model using a Bayesian linear mixed model, and use the estimates from this model to simulate data. 

```{r warning=FALSE,message=FALSE}
library(lingpsych)
data("dillonE1ttnested")
library(brms)
postD_mE1ttnested<-posterior_samples(dillonE1ttnested)
beta1<-postD_mE1ttnested$b_Intercept
beta2<-postD_mE1ttnested$"b_dep"
beta3<-postD_mE1ttnested$"b_gram"
beta4<-postD_mE1ttnested$"b_dep:gram"
beta5<-postD_mE1ttnested$"b_int.rg"
beta6<-postD_mE1ttnested$"b_int.ag"
## the two estimates of interest:
beta7<-postD_mE1ttnested$"b_int.ru"
beta8<-postD_mE1ttnested$"b_int.au"
betameans<-c(beta1,beta2,beta3,beta4,beta5,beta6,beta7,beta8)

cor<-posterior_samples(postD_mE1ttnested,"^cor")
sd<-posterior_samples(postD_mE1ttnested,"^sd")
sigma<-posterior_samples(postD_mE1ttnested,"sigma")

item_re<-posterior_samples(postD_mE1ttnested,"^r_item")
subj_re<-posterior_samples(postD_mE1ttnested,"^r_subj")

library(SIN)
sds<-colMeans(sd)
cors<-colMeans(cor)
sig<-mean(sigma$sigma)
## assume intermediate correlation:
cormat<-matrix(rep(0.5,8*8),ncol=8)
diagmat<-diag(rep(0.5,8))
cormat<-cormat+diagmat
Sigma_u<-SIN::sdcor2cov(stddev=sds[8+(c(1,2,3,8,7,5,6,4))],
                          corr=cormat)
Sigma_w<-SIN::sdcor2cov(stddev=sds[c(1,2,3,8,7,5,6,4)],
                          corr=cormat)
```

Finally, we set all the slopes to zero, and compute Type I error under different conditions.

```{r}
## set all effects to 0:
betameans[2:8]<-0
```

We will fit a model for a full variance-covariance matrix for both subjects and items. We avoid fitting the correlation parameters, because these will be difficult to estimate with the sample size (40 subjects and 48 items) used in the @@Dillon-EtAl-2013 study. To illustrate the effect of mis-specification of the likelihood function, we will fit the simulated data to both log and raw reading times.

```{r fullmodelDillonEtAl2013,eval=FALSE}
## full model:
nsim<-100
significantlog<-significantraw<-matrix(0,nsim,ncol=7)
## record failed convergences, remove these later:
failedlog<-failedraw<-rep(0,nsim)

for(i in 1:nsim){
##  print(paste("********simulation ",i,"********",sep=""))
  simdat<-gen_fake_lnorm2x2x2(beta = betameans,
                               Sigma_u=Sigma_u,Sigma_w=Sigma_w,
                               sigma_e=sig,nitem=48,
                               nsubj=40)

  m_fakelog<-lmer(log(rt)~Dep+Gram+DepxGram+Int_gram_refl+Int_gram_agr+
                 Int_ungram_refl+Int_ungram_agr+
                 (1+Dep+Gram+DepxGram+Int_gram_refl+Int_gram_agr+
                    Int_ungram_refl+Int_ungram_agr||subj)+
                 (1+Dep+Gram+DepxGram+Int_gram_refl+Int_gram_agr+
                    Int_ungram_refl+Int_ungram_agr||item),simdat,
               control=lmerControl(optimizer="nloptwrap", calc.derivs = FALSE))

    m_fakeraw<-lmer(rt~Dep+Gram+DepxGram+Int_gram_refl+Int_gram_agr+
                 Int_ungram_refl+Int_ungram_agr+
                 (1+Dep+Gram+DepxGram+Int_gram_refl+Int_gram_agr+
                    Int_ungram_refl+Int_ungram_agr||subj)+
                 (1+Dep+Gram+DepxGram+Int_gram_refl+Int_gram_agr+
                    Int_ungram_refl+Int_ungram_agr||item),simdat,
               control=lmerControl(optimizer="nloptwrap", 
                                   calc.derivs = FALSE))

  ## ignore failed trials (log rts model)
    if(any( grepl("failed to converge", m_fakelog@optinfo$conv$lme4$messages) ) ||
       any( grepl("singular", m_fakelog@optinfo$conv$lme4$messages) )){
      failedlog[i]<-1
    } else{
  ## first slope:
  if(abs(summary(m_fakelog)$coefficients[2,3])>2){
    significantlog[i,1]<-1
  } else {significantlog[i,1]<-0}
  ## second slope:
  if(abs(summary(m_fakelog)$coefficients[3,4])>2){
    significant[i,2]<-1
  } else {significantlog[i,2]<-0}
  ## third slope:
  if(abs(summary(m_fakelog)$coefficients[4,4])>2){
    significantlog[i,3]<-1
  } else {significantlog[i,3]<-0}
  ## fourth slope:
  if(abs(summary(m_fakelog)$coefficients[5,4])>2){
    significant[i,4]<-1
  } else {significantlog[i,4]<-0}
  ## fifth slope:
  if(abs(summary(m_fakelog)$coefficients[6,4])>2){
    significantlog[i,5]<-1
  } else {significantlog[i,5]<-0}
  ## sixth slope:
  if(abs(summary(m_fakelog)$coefficients[7,4])>2){
    significantlog[i,6]<-1
  } else {significantlog[i,6]<-0}
  ## seventh slope:
  if(abs(summary(m_fakelog)$coefficients[8,4])>2){
    significantlog[i,7]<-1
  } else {significantlog[i,7]<-0}
    }
    
    ## ignore failed trials (raw rts model)
    if(any( grepl("failed to converge", m_fakeraw@optinfo$conv$lme4$messages) ) ||
       any( grepl("singular", m_fakeraw@optinfo$conv$lme4$messages) )){
      failedraw[i]<-1
    } else{
  ## first slope:
  if(abs(summary(m_fakeraw)$coefficients[2,3])>2){
    significantraw[i,1]<-1
  } else {significantraw[i,1]<-0}
  ## second slope:
  if(abs(summary(m_fakeraw)$coefficients[3,4])>2){
    significant[i,2]<-1
  } else {significantraw[i,2]<-0}
  ## third slope:
  if(abs(summary(m_fakeraw)$coefficients[4,4])>2){
    significantraw[i,3]<-1
  } else {significantraw[i,3]<-0}
  ## fourth slope:
  if(abs(summary(m_fakeraw)$coefficients[5,4])>2){
    significant[i,4]<-1
  } else {significantraw[i,4]<-0}
  ## fifth slope:
  if(abs(summary(m_fakeraw)$coefficients[6,4])>2){
    significantraw[i,5]<-1
  } else {significantraw[i,5]<-0}
  ## sixth slope:
  if(abs(summary(m_fakeraw)$coefficients[7,4])>2){
    significantraw[i,6]<-1
  } else {significantraw[i,6]<-0}
  ## seventh slope:
  if(abs(summary(m_fakeraw)$coefficients[8,4])>2){
    significantraw[i,7]<-1
  } else {significantraw[i,7]<-0}
  }  
}

## In 100 simulations, what is the proportion of times that *at least one* of the effects is significant?
## Log rt model:
save(significantlog,file="maxmodelsiglog.rda")
countslog<-table(rowSums(significantlog)>0)
## Type I error under multiple comparisons:
countslog[2]/sum(countslog)

## Raw rt model:
save(significantraw,file="maxmodelsigraw.rda")
countsraw<-table(rowSums(significantraw)>0)
## Type I error under multiple comparisons:
countsraw[2]/sum(countsraw)
```

```{r echo=FALSE}
load("data/maxmodelsiglog.rda")
load("data/maxmodelsigraw.rda")
```

What we see below is that (a) with the log reading time model, we get an inflation of Type I error when we compare multiple hypotheses, (b) model mis-specification (using a Normal likelihood when a log-normal should have been used) will further inflate Type I error.

```{r}
countslog<-table(rowSums(significantlog)>0)
## Type I error under multiple comparisons:
countslog[2]/sum(countslog)

countsraw<-table(rowSums(significantraw)>0)
## Type I error under multiple comparisons + model mis-specification:
countsraw[2]/sum(countsraw)
```

## The practical implications   

What all this means is that when we have multiple hypothesis within a linear mixed model, we cannot reasonably assume that Type I error is 0.05 for each of the hypothesis tests.

## Summary

This chapter demonstrates, through simulations, how Type I error can get inflated due to (a) an overly simplified random effects structure in linear mixed models, (b) mis-specification of the likelihood, and (c) multiple comparisons.

## Further reading

@barr2013 is an important reference on the topic of Type I error inflation through overly simplified specification of the random effects structure of the linear mixed model, although also see @hannesBEAP and @BatesEtAlParsimonious. @MalsburgAngele2016 investigate the problem of Type I error through multiple comparisons in linear mixed models, focusing on eyetracking reading data.