stats_07_linear-mixed-effects-models.Rmd

---
title: "Linear Mixed Effects Models"
author: "Marius 't Hart"
output:
  html_notebook: default
  html_document:
    df_print: paged
  pdf_document: default
---

```{r setup, cache=FALSE, include=FALSE}
library(knitr)
opts_chunk$set(comment='')
```

Sometimes you have lots of missing data and no way to impute this easily, but still would like to do an ANOVA. ANOVA's don't handle missing data, but there is an alternative: linear-mixed effects models (LME). First, we'll get some data from our own work, that we analyzed with an LME for publication.

We need some packages to do LME. The `lme4` package actually fits the models, and the `lmerTest` package gets some p-values.

```{r}
#install.packages(c('lme4','lmerTest'))
library('lme4')
library('lmerTest')
```

Do we need this? Yes. Yes, we do.

```{r}
default.contrasts <- options('contrasts')
options(contrasts=c('contr.sum','contr.poly'))
```

We also need some example data:

```{r}
exp <- read.csv('data/exposure_localization.csv', stringsAsFactors = FALSE)

# some columns need to be factors:
exp$participant   <- factor(exp$participant)
exp$rotated_b     <- factor(exp$rotated_b)
exp$passive_b     <- factor(exp$passive_b)
exp$handangle_deg <- factor(exp$handangle_deg)

```

This is from our paper [Motor Learning Without Moving](https://doi.org/10.1371/journal.pone.0221861) and it's the data used to create Fig 3A. So what's in the data set?

```{r}
str(exp)
```

So it's a data frame with 6 columns and 588 rows. The `group` column is not so meaningful, as I'm only using this one group for this notebook. What else?

A brief explanation of the data: this is a measurement of where people (mis)localize their unseen hand before and after a visuomotor adaptation task (the data in the `taperror_deg` column). In this case, the trained hand was moved by a robot in such a way that the cursor that representing the hand position would always go straight to the target: there are no motor errors to learn from. Importantly, there is also only 1 target: at 45 degrees. In the rotated session, using the same single target, the cursor did move in a different direction than the hand. Or rather, the hand didn't really move to the target, but the cursor did. The hand moved 30 degrees counter clockwise from the target. This creates a visuo-proprioceptive discrepancy that should recalibrate people's felt hand position, or what we sometime call "proprioception". Here we measure that with what is called "hand localization". So after training, the right hand is moved out, but there is no cursor. Instead, an arc appears, and people use their visible, untrained hand on a touch screen to indicate where they think their right hand is: they localize their right hand. If this localization measure reflects proprioception, it should shift with the training (we code whether the measure was done after rotated training or not in the `rotated_b` column). Interestingly, the movement of the right hand in the localization task (not the training) can be done voluntarily (actively) so that people have an efferent-based prediction, as well as their proprioception. Or it can be done passively, so that only proprioception should be available (the `passive_b` column). In this task we don't expect the availability of efferent-based predictions to matter. Since the active movements were in voluntarily chosen directions, and the passive movements went to the exact same spots, we had to estimate the hand localization at a set of actual hand locations (the `handangle_deg` column) by some form of interpolation. In this case we used kernel smoothing interpolation. However, not all participants (the `participant` column) spread out their reaches as much as the others, so there is some missing data. The effects should be larger at the trained target and then fall of on either side. This means that the missing data is hard to impute or extrapolate.

For most analyses in the paper, we did not use the 15 degree point, as it had perhaps too much missing data so that a few participants might have biased the analyses. Maybe that was a bad choice, but it's what we did.

There should probably be a figure here, showing the data.

```{r fig.width=7.5, fig.height=3}
source('R/exposureLocalization.R')
plotExposureLocalization(target='inline', remove15=FALSE)
```

So it looks like the manipulation had an effect: after rotated exposure training, the localization is shifted relative to after aligned training. I'm not really sure there is a difference between the localization shifts in active and passive. But... the first step is to verify statistically that the training had an effect. That's _"all"_ that will be done in this tutorial.

# ANOVA

You could skip this first section, as it merely shows why using an ANOVA would not be the best use of the data.

Let's see what happens if we _do_ do an ANOVA, using the ez package:

```{r error=TRUE}
#install.packages('ez')
# library(ez)
# ezANOVA(dv=taperror_deg, wid=participant, within=c(rotated_b, passive_b, handangle_deg), data=exp)
library(afex)
afex::aov_ez(dv='taperror_deg', id='participant', within=c('rotated_b', 'passive_b', 'handangle_deg'), data=exp)
```

That throws a big error because of the missing values (NAs) in taperror_deg:

```{r}
length(which(is.na(exp$taperror_deg)))
```

There are 29 NAs... which is less than 5% of the data. Should be acceptable. We could remove participants if one of their data points is a missing value:

```{r}
expPA <- exp[which(!exp$participant %in% unique(exp$participant[which(is.na(exp$taperror_deg))])),]
ezANOVA(dv=taperror_deg, wid=participant, within=c(rotated_b, passive_b, handangle_deg), data=expPA)
```

Now we get output! That is great.

```{r}
cat(sprintf('Removed %d rows, or %0.2f%% of the data.\n',
            dim(exp)[1] - dim(expPA)[1] - length(which(is.na(exp$taperror_deg))),
            100 - (dim(expPA)[1]/(dim(exp)[1] - length(which(is.na(exp$taperror_deg))))) * 100
            )
    )
```

However, even though there were only 29 missing values, we had to remove 167 additional rows which is about 30% of the data. So what is that ANOVA really telling us? Maybe we can do better?

Let's remove all the hand angles that have some missing data instead:

```{r}
expHA <- exp[which(!exp$handangle_deg %in% unique(exp$handangle_deg[which(is.na(exp$taperror_deg))])),]

cat(sprintf('Removed %d rows, or %0.2f%% of the data.\n',
            dim(exp)[1] - dim(expHA)[1] - length(which(is.na(exp$taperror_deg))),
            100 - (dim(expHA)[1] / (dim(exp)[1] - length(which(is.na(exp$taperror_deg))))) * 100
            )
   )
```

That's even more data that gets remove: about 40%. Hmmm. Those NA's are not missing at random, it's because we can't really extrapolate from the present data. Also, there are only data points missing at 15, 25 and 35 degrees. Maybe we should remove the most egregious hand angle (15 degrees) as we did for some analyses in the paper, and then check which participants still have missing data. Perhaps then we won't have to remove so much data?

```{r}
expH <- exp[-which(exp$handangle_deg == 15),]
expHP <- expH[which(!expH$participant %in% unique(expH$participant[which(is.na(expH$taperror_deg))])),]

cat(sprintf('Removed %d rows, or %0.2f%% of the data.\n',
            dim(exp)[1] - dim(expHP)[1] - length(which(is.na(exp$taperror_deg))),
            100 - (dim(expHP)[1] / (dim(exp)[1] - length(which(is.na(exp$taperror_deg))))) * 100
            )
   )
```

That's the least amount of data removed. Here is the ANOVA:

```{r}
ezANOVA(dv=taperror_deg, wid=participant, within=c(rotated_b, passive_b, handangle_deg), data=expHP)
```

Some might think that is acceptable, but it does add up to over 1/5th of the good data that had to be removed. And this might skew our results: it's 1 of 7 hand angles and 3 of 21 participants that no longer contribute any data at all. Not super-kosher...

So we should probably do some other analysis.

# Linear Mixed Effects models

Linear mixed effects models to the rescue!

First, if you want to exactly reproduce the analysis that was published, you should set `remove15` to `TRUE` in the following chunk, but it works fine with all the data as well:

```{r}
remove15 <- FALSE
if (remove15) {
  exp <- exp[-which(exp$handangle_deg == 15),]
}
```

Let's fit the linear mixed effects model:

```{r}
exp_model_lmer <- lmerTest::lmer(taperror_deg ~ rotated_b * passive_b * handangle_deg - (1|participant),
                       na.action = na.exclude,
                       data = exp,
                       REML = TRUE,
                       control = lmerControl(optimizer ="Nelder_Mead")
                       )
```

Let's break this command down.

In LMEs there are two types of `effects`: fixed effects (which are usually the factors you'd be interested in in an ANOVA) and random effects (which you'd usually not care about much - but this may depend on your question and data). In this case the "formula" is:

`taperror_deg ~ rotated_b * passive_b * handangle_deg - (1|participant)`

- `taperror_deg` is the dependent variable, that the model should explain
- `rotated_b * passive_b * handangle_deg` are the fixed effects and the `*` tells the `lmer()` function that we want all combinations of their levels: these factors may interact 
- `(1|participant)` says that participant is a random effect: we get (intercept) coefficients for each, but they don't interact with the fixed effects, this is akin to (but not the same as) doing a repeated-measures ANOVA

We also tell the function what data frame to use, to exclude NA's, and to use REML (not ML: maximum-likelihood).

Finally, I've used a simple Nelder-Mead optimizer, because using the standard optimizer gives a warning that the model did not converge (try it), although the result is more or less the same.

First, we should probably look at this model:

```{r}
summary(exp_model_lmer)
```

And it's residuals:

```{r}
plot(exp_model_lmer)
```

In this case there might be some systematic errors in the model, as the residuals seem close to zero around 0 and around -8, but in between they might be a tiny bit higher. I won't check this here, but it would mean that the fixed effects (predictors / factors) don't fully explain the data, maybe there is some additional effect going on, or we should put the data on some other scale. However, by and large it seems to work reasonably well. Much better than no model.

This is pretty cool, but there are terms for the combinations of every level in the factors. Maybe this provides a lot more information, but it's not so comparable to the trusty, old ANOVA.

LME doesn't give you ANOVA-like output by design. The `lme4` package used to give ANOVA-like output, but it turns out that there is no good 1:1 mapping from LME to ANOVA-like output. There are several methods to approximate it, but opinions differ on which one is the best, the most appropriate to certain situations, or even which ones are good. So the `lme4` package dropped it. You can however still get that kind of output, using the `lmerTest` package, but as of yet this might not be widely accepted, as in: some reviewers may disagree with using this at all.

The `lmerTest` package allows using the "Satterthwaite" method and overwrites (or extends?) the `anova()` function:

```{r}
print(anova(exp_model_lmer,ddf='Satterthwaite',type=3))
```

Another popular method is "Kenward-Roger" which depends on the `pbkrtest` package:

```{r}
#install.packages('pbkrtest')
print(anova(exp_model_lmer,ddf='Kenward-Roger',type=3))
```

The exact numbers are slightly different, but the general pattern is the same.

As a third method to convert an LME model to ANOVA-like output, we'll try out a Chi-squared method, which doesn't use the `lme4` and `lmerTest` package. It uses the `lme()` function from the `nlme` package, and then get's the Chi-squared approximation by using the `Anova()` function from the `car` package.

```{r}
#install.packages(c('nlme','car'))
library(nlme)
library(car)
```

You will see a warning that something from `lme4` is overwritten. This is why you should use namespaces.

```{r}
print(car::Anova(nlme::lme(taperror_deg ~ rotated_b * passive_b * handangle_deg, random = ~1|participant, na.action=na.exclude, data=exp), type=3))

```

Again, this provides the same overall pattern of p-values, even when the exact values are somewhat different.

Back to `lme4`. If you don't want to use the Nelder-Mead optimizer, but want to do something a little more advanced, you can grab some different optimizers from the `optimx` package:

```{r}
#install.packages('optimx')
library(optimx)
```

First, let's see that warning without the Nelder-Mead optimizer:

```{r}
exp_model_lmer <- lmerTest::lmer(taperror_deg ~ rotated_b * passive_b * handangle_deg - (1|participant),
                       na.action = na.exclude,
                       data = exp,
                       REML = TRUE
                       )
```

That sounds scary, but **in this case** it still gives more or less the same output:

```{r}
print(anova(exp_model_lmer,ddf='Satterthwaite',type=3))
```

No guarantee for other cases though! If you get a warning like that, you should see if using different optimizers or changing the control settings gives you very different models. If so, you should probably not trust the non-converging model.

Now, let's use an optimizer from `optimx` to see if that makes a difference:

```{r}
exp_model_lmer <- lmerTest::lmer(taperror_deg ~ rotated_b * passive_b * handangle_deg - (1|participant),
                                 na.action = na.exclude,
                                 data = exp,
                                 REML = TRUE,
                                 control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B'))
)
```

This gives no warnings.

```{r}
print(anova(exp_model_lmer,ddf='Satterthwaite',type=3))
```

And still the same pattern of results.

Now we still need to interpret this... If I'm correct, then using the right contrast options (that I sneakily set in the third chunk) and using type 3 sums of squares, you _can_ interpret main effects in the presence of interactions. You just have to be careful.

Source for that bold statement: [Matt's Stats n stuff](https://mcfromnz.wordpress.com/2011/03/02/anova-type-iiiiii-ss-explained/).

That would mean that providing rotated feedback during exposure training has an effect on hand localization, so does the angle where the hand actually is, and these two factors also interact. However, for this test, I just wanted to say that first thing: does the training with rotated feedback change the hand localizations? For this data, I don't actually care if that happens at all hand positions or in both tasks, as long as it changes at some hand angle (preferably close to the trained target, and more so for the active localization instead of the passive). If it happens at all, we can then calculate difference scores (localization after rotated exposure - localization before rotated exposure, for all hand angles and active and passive localization) and then look at the _magnitude_ of that change across the conditions. In the model I wanted to include all sorts of other possible effects, but I can do the model without interactions:

```{r}
exp_model_lmer <- lmerTest::lmer(taperror_deg ~ rotated_b + passive_b + handangle_deg - (1|participant),
                                 na.action = na.exclude,
                                 data = exp,
                                 REML = TRUE,
                                 control = lmerControl(optimizer ='optimx', optCtrl=list(method='L-BFGS-B'))
)

print(anova(exp_model_lmer,ddf='Satterthwaite',type=3))
```

Can I still draw the same conclusions as before, now that there are main effects without an interaction? Maybe not, as these main effects might be an artefact of the simpler model.

However, given all the tests done here, I'm pretty sure that exposure training did have an effect on localization responses. That means I can continue looking at the size of that effect in other analyses, but this is where the tutorial stops.