![cc](pix/cc-attrib-nc.png)
Licensed under the 
[Creative Commons attribution-noncommercial license](http://creativecommons.org/licenses/by-nc/3.0/).
Please share \& remix noncommercially, mentioning its origin.

## Preliminaries

```{r setup,echo=FALSE,message=FALSE}
## ignore this stuff ...
library("knitr")
opts_chunk$set(fig.align="center",
               fig.width=5,fig.height=5,tidy=FALSE,message=FALSE)
knit_hooks$set(basefig=function(before, options, envir) {
                   if (before) {
                       par(bty="l",las=1)
                   } else { }
               })
```

![cc](pix/cc-attrib-nc.png)
<!---
(http://i.creativecommons.org/l/by-nc-sa/3.0/88x31.png)
--->
Licensed under the 
[Creative Commons attribution-noncommercial license](http://creativecommons.org/licenses/by-nc-sa/2.5/ca/).
Please share \& remix noncommercially, mentioning its origin.

Load packages:
```{r pkgs}
## modeling packages
library("lme4")     ## basic (G)LMMs
library("nlme")     ## more LMMs
library("afex")     ## helper functions
library("emmeans")

## graphics
library("ggplot2"); theme_set(theme_bw())
## squash panels together
zero_margin <- theme(panel.spacing=grid::unit(0,"lines")) 
library("lattice")
library("dotwhisker")  ## coefficient plots

## data manipulation
library("broom")
library("broom.mixed")
```

<!-- (trying to fix heading display) -->

## Linear mixed model: starling example

Data from Toby Marthews, [r-sig-mixed-models mailing list](http://tolstoy.newcastle.edu.au/R/e13/help/11/01/1364.html): you can find it [here](https://github.com/bbolker/upr_2016/tree/master/lab/data).

- `subject`: individual bird
- `roostsitu` (tree, nest-box, inside, other): location of nest box
- `mnth` (Nov, Jan): month
- `stmass`: mass of bird

### Graphics: spaghetti plots

```{r starling1,fig.height=4,fig.width=7}
## cosmetic tweaks
load("../data/starling.RData")   ## loads "dataf"
ggplot(dataf,aes(x=mnth,y=stmass))+
  geom_point()+
  geom_line(aes(group=subject))+  ## connect subjects with a line
  facet_grid(.~roostsitu)+        ## 1 row of panels by roost
  zero_margin                         ## squash together
```

You should also try
```{r altplot,fig.keep="none"}
ggplot(dataf,aes(mnth,stmass,colour=roostsitu))+
    geom_point()+
    geom_line(aes(group=subject))
```

for comparison: what are the tradeoffs?

Overkill for this data set, but sometimes it can be useful to put every individual in its own facet:

```{r altplot2,fig.keep="none"}
## reorder individuals by magnitude of difference
dataf <- transform(dataf,subject=reorder(subject,stmass,FUN=diff))
ggplot(dataf,aes(mnth,stmass,colour=roostsitu,group=subject))+
    geom_point()+geom_line()+facet_wrap(~subject)+
    zero_margin
```

It's pretty obvious that the starting (November) mass varies among roost situations (tree/nest-box/etc.), and that mass increases from November to January, but we might like to know if the slopes differ among situations. That means our fixed effects would be `~roostsitu*mnth`, with our attention focused on the `roostsitu:mnth` (interaction) term.  For random effects, we can allow both intercept (obviously) and slope (maybe) to vary among individuals, via `(1+mnth|subject)` or equivalent ... in this case, because measurements are only taken in two months, we can also write the random term as `(1|subject/mnth)`.

**However**, it turns out that we can't actually estimate random slopes for this model, because every individual is only measured twice. That means that the variance in the slopes would be completely confounded with the residual variance.

If we fit this with lme4
```{r lmerfit1,error=TRUE}
lme1 <- lmer(stmass~mnth*roostsitu+(1|subject/mnth),data=dataf)
```

We get an error. We *can* actually do this (`lmerControl(check.nobs.vs.nlev="ignore")`; see `?lmerControl` for details), and it will actually work, but it will pick an arbitrary division between the residual variance and the `month:subject` variance.

But now let's forget about including the (unidentifiable) random slopes.
```{r lmer2}
lmer2 <- lmer(stmass~mnth*roostsitu+(1|subject),data=dataf)
```
We can now get estimates, although the subject-level random effects are *very* uncertain: see `confint(lmer2)`

Walking through the output:
```{r echo=FALSE}
lmer2out <- capture.output(summary(lmer2))
cat(lmer2out[1:5],sep="\n")
```

- This gives us the formula, reminds us that the model was fitted by restricted maximum likelihood, and gives us the "REML criterion" (equivalent to -2 log likelihood for a model fitted by maximum likelihood).

```{r lmer2sum2,echo=FALSE}
cat(lmer2out[6:9],sep="\n")
```

This is a quick reminder of the scaled residuals; these should be
approximately unskewed (median $\approx 0$) and the extremes should
be somewhere around $\pm 2$ (bigger in a big data set \ldots)

```{r lmer2sum3,echo=FALSE}
cat(lmer2out[11:14],sep="\n")
```

- This tells us about the random effect - 
it tells us the variances and standard deviations
of the random effect 
(the standard deviations
are *not* an estimate of the uncertainty of the estimate -- the estimate
itself *is* a variance, or standard deviation), and the
variance and standard deviation
of the residual variance.  We can see that the standard deviation
of the intercept among subjects is much less than the 
standard deviation of the residual error ...

```{r lmer2sum4,echo=FALSE}
cat(lmer2out[15],sep="\n")
```

- The listed number of observations and groups is *very* useful for double-checking that the random effects grouping specification is OK. Are the numbers of groups what we expected from the experimental design?

```{r lmer2sum5,echo=FALSE}
cat(lmer2out[17:25],sep="\n")
```

- the standard fixed-effect output. `lmer` does *not* tell us the denominator degrees of freedom for the test (although we can get a rough idea of importance/significance fro the $t$ statistics; e.g. $t>2.5$ will be significant at $p<0.05$ for 6 or more degrees of freedom). We'll come back to this in the inference section.

## Diagnostic plots

Basic diagnostic plots: fitted vs. residuals, coloured according to the `roostsitu` variable (`col=dataf$roostsitu`), point type according to month (`pch=`).

```{r lmediag}
plot(lmer2,col=dataf$roostsitu,pch=dataf$mnth,id=0.05)
```

The `id` argument specifies to mark points that are beyond the specified Normal
confidence limits (i.e. $\alpha$-values) (these are not outliers in any formal
statistical sense but might be points you want to look at if they are otherwise
interesting). By default they are labeled by random-effect group; use `idLabels`
to change the default labels, and `idLabels=~.obs` to label by observation number.

You can get a scale-location plot by specifying

```{r lmediag2,fig.keep="none"}
plot(lmer2,sqrt(abs(residuals(.))) ~ fitted(.),type=c("p","smooth"))
```

Adding the smoothed line is helpful because uneven sampling can
influence your perception of the pattern.

Q-Q plot (a straight line indicates normality)
```{r qqnorm}
qqmath(lmer2,col=dataf$roostsitu,pch=dataf$mnth)
```

We can also do this with `broom.mixed::augment()` and `stat_qq()` (although there are also some some limitations - `colour` argument doesn't actually work...)
```{r ggqq,fig.keep="none",cache=TRUE}
ggplot(augment(lmer2),
       aes(sample=.resid/sd(.resid)))+  ## scale to variance=1
    stat_qq(aes(group=1,colour=roostsitu))+
    geom_abline(intercept=0,slope=1)
```

There are some deviations here, but not enough that I would worry very much.
In particular, the distribution is slightly *thin-tailed* (the smallest residuals are largest than expected, and the largest residuals are smaller than expected),
which would make the results slightly conservative (a fat-tailed distribution would make them anticonservative).

Boxplots of residuals subdivided by `roostsitu` (you have to put the grouping variable on the *left* side of the formula here):

```{r diagbox}
plot(lmer2,roostsitu~resid(.))
```

Might be easier with `augment`:

```{r diagbox_gg,fig.keep="none"}
aa <- augment(lmer2)
ggplot(aa,aes(roostsitu,.resid))+
    geom_boxplot()+coord_flip()
```

Plot random effects to look for outliers:

```{r ranef}
dotplot(ranef(lmer2,condVar=TRUE))
```

Or with `tidy` + `ggplot`:
```{r ranef2}
tt <- tidy(lmer2,effects="ran_vals")
ggplot(tt,aes(level,estimate))+
    geom_pointrange(aes(ymin=estimate-1.96*std.error,
                        ymax=estimate+1.96*std.error))+
    coord_flip()
```

In general, try plotting residuals and $\sqrt{\textrm{abs}(r)}$ both as a
function of the fitted values and as a function of particularly important predictors.

- if you have your own example data set handy that you can write a simple model for, run these diagnostics and interpret them
- the `hatvalues()` function can also be useful (although all hat values are identical here)
- also see the `influence.ME` package.  Try `plot(influence(lmer2,group="subject"))` and interpret the results ...

### Try it yourself: tundra carbon example

- Download the file `tundra_agg.rda` (from the data directory given above). The `GS.NEE` records the net ecosystem exchange (in grams of carbon/m^2/year) in a given year and site; `cYear` is a centered version of the year; `n` records the number of observations taken in a given site/year (the NEE value given is an average of these observations).
- Draw some plots (with `ggplot2` if you can, otherwise with base R graphics) to inspect the data
- We are interested in the change in NEE over time (i.e., a random-slopes model). Run an `lmer` fit (note this data set does not have the issue of the starling data set), draw some preliminary conclusions, and inspect the plot diagnostics.

## Inference

One good way to assess the results of a model fit is to look at a coefficient plot:

```{r coefplot}
dwplot(lmer2)+geom_vline(xintercept=0,lty=2)
```

(Note that the intercept is automatically dropped; this is usually the right choice. Use `show_intercept=TRUE` to have the intercept included.) For more complicated coefficient plots, you can use the combination of the `broom.mixed` package to get the parameters in a useful form and `dwplot` to draw the plot.

Stop and explain to yourself what these parameters mean.  If you're not sure, try resetting the base level of the `roostsitu` factor: `dataf2 <- transform(dataf,roostsitu=relevel(roostsitu,ref="other"))`, predict what will happen to the results, and re-run the analysis.)

**Exercise**: This is actually a slightly trivial example, because there are only two measurements for each individual. Thus we can actually analyze the slopes and get the same answers by using a paired analysis, i.e. by computing the mass *difference* in subjects and then analyzing with a single-level linear model. (In the simplest case this would reduce to a paired $t$-test, but in this case it is a 1-way ANOVA on `roostsitu`.

Rearrange the data to get differences by subject:
```{r diffmass}
dataf2 <- aggregate(stmass~roostsitu+subject,
               data=dataf,
               FUN=function(x) x[2]-x[1])
```
This says to aggregate the data frame `dataf` grouped by `roostsitu` and `subject`; for each group, compute the difference of the masses.

Draw a picture (boxplot+beeswarm plot):
```{r plotdiffs,message=FALSE}
ggplot(dataf2,aes(x=roostsitu,y=stmass))+geom_boxplot()+
   geom_dotplot(binaxis="y",stackdir="center",fill="red",alpha=0.5,
                binwidth=0.5)
```
As you can see, `geom_dotplot()` adds a horizontal dot-plot for each group
(see the documentation [`?geom_dotplot`] for more details).

- Analyze the data with `lm` and convince yourself that the estimates (fixed-effect coefficients, $t$ statistics, etc.) are equivalent to those found from the previous analysis.
- It is also possible to rearrange and summarize the data to test the difference in intercepts, or to estimate the among-individual variance in intercepts (how?)

```{r lmcheat,echo=FALSE,results="hide"}
summary(lm(stmass~roostsitu,dataf2))
```

## Inference part 2

What about $p$-values or more sophisticated confidence intervals?
The default confidence intervals shown by `dwplot()` are just Normal
confidence intervals (not taking into account any of the effects that
make these approximate for mixed models).

We can get profile confidence intervals via `confint(lmer2)`,
but the confidence intervals are almost identical to those
from `confint(lmer2,method="Wald")` (the default shown by
`dwplot()`.

(Insert long discussion of degrees of freedom and $p$-values here.)

We can use `lmerTest` or `afex::mixed` (which both wrap functions
from the `pbkrtest` package) to get "denominator degrees of freedom"/$p$-values.

- `lmerTest` wraps the `lmer` function to compute some additional information,
and modifies the `summary()` function:
```{r lmertest,warning=FALSE}
library("lmerTest")
lmer2R <- lmer(stmass~mnth*roostsitu+(1|subject),data=dataf)
detach("package:lmerTest")
```
Capturing just the coefficient table from `summary(lmer2R)`:
```{r lmer2Rsum_out,echo=FALSE}
s1 <- summary(lmer2R,digits=3,ddf="Satterthwaite")
lmer2Rout <- capture.output(s1)
cat(lmer2Rout[18:27],sep="\n")
```
By default `lmerTest` uses the Satterthwaite approximation to the
degrees of freedom.  You can also use `summary(lmer2R,ddf="Kenward-Roger")`
to get a slightly more accurate (and slower) estimate, but in this case
the answers are basically identical.

We conclude that the interactions are not doing much,
but there's definitely an effect of the roosts and months.
This agrees with the picture from `dwplot` (above).

However, we probably want to test the overall effect of the interactions, not the individual levels.
Here are the type II (sequential) ANOVA results:
```{r anovalme2}
anova(lmer2)
```
`lmer` doesn't give denominator degrees of freedom ...

If we want to evaluate the *marginal* sums of squares, i.e. dropping one term at a time from the model, when there are interactions,
we usually want to change the model to use sum-to-zero contrasts: 

```{r contrlme}
lmer2B <- update(lmer2,            
          contrasts=list(mnth="contr.sum",roostsitu="contr.sum"))
```

The alternative approach is to use `options(contrasts=c("contr.sum","contr.poly"))`, then refit the model, but I prefer to use the `contrasts` argument because it is more explicit.

- Use `afex::mixed` to do "type 3" tests of all effects in
a model (using Kenward-Roger to get df by default).
```{r mixed}
afex::mixed(stmass~mnth*roostsitu+(1|subject),data=dataf)
```
- Use `anova()` to do a likelihood ratio test on
individual pairs of models.
- Use `drop1()` to drop terms from the model one at a time

**be careful when doing "marginal"/"type 3" tests in the presence of interactions!**

In this case the results ($F$ values) are identical **because the original design is balanced (hence, orthogonal)**.  Not true if the data are (1) unbalanced (which is often true of ANOVA [categorical-predictor] designs, and almost always true of regression designs) or (2) GLMM or nonlinear.

The explicit model-comparison approach uses a likelihood ratio test rather than an $F$ test (i.e., it does not correct for the fact that the denominator sums of squares is estimated with error). In this case it hardly matters.

```{r testmodels}
lmer2C <- update(lmer2B,REML=FALSE)
lmer2D <- update(lmer2C,. ~ . - mnth:roostsitu)
anova(lmer2C,lmer2D)
```

If we now want to use the model-comparison approach on the reduced (no-interaction) model to test the significance of `roostsitu`, we can use `update` to drop the `roostsitu` effect, but we also have to make sure to update the `contrasts` argument so that it only refers to predictors that remain in the reduced model (otherwise, we get an error).

```{r test2,warning=FALSE}
drop1(lmer2D,test="Chisq") ## ignore warnings
```

If we want to test the random effect, we would in principle remove the random effect and test with `anova`, but this is a bit problematic here: `lmer` can't fit a model without any random effects.

Let's try this with `lm`:
```{r anovacmp}
lm1 <- lm(stmass~mnth*roostsitu,data=dataf)
(a1 <- anova(lmer2C,lm1))
```

We can also deduce the number of degrees of freedom from
standard rules:
```{r datatab}
with(dataf,table(mnth,roostsitu))
```
If we think about this in terms of the paired $t$-test,
there are 40 comparisons and 4 parameters (one for each
roost site)= 36 df.

If you wanted to compute the $p$-values by hand, you could:
```{r pvalues}
## here we use lme4:::anova.merMod to force R to use the anova
## method from lme4 rather than the one from lmerTest
a2 <- lme4:::anova.merMod(lmer2B)
fvals <- a2[,"F value"]
numdf <- a2[,"Df"]
dendf <- 36
pf(fvals,numdf,dendf,lower.tail=FALSE)
```

Alternatively, we can try a parametric bootstrap: the `pbkrtest` package can do this, or we can just set it up by hand:

```{r pboot,echo=FALSE,cache=TRUE,warning=FALSE}
lmer3 <- lmer(stmass~mnth*roostsitu+(1|subject),data=dataf)
lmer4 <- lmer(stmass~mnth+roostsitu+(1|subject),data=dataf)

pboot <- function(m0,m1) {
  s <- simulate(m0)
  L1 <- logLik(refit(m1,s[[1]]))
  L0 <- logLik(refit(m0,s[[1]]))
  2*(L1-L0)
}
pboot(lmer4,lmer3)
boothist <- replicate(1000,pboot(lmer4,lmer3))
library(plyr)
boothist <- rlply(1000,pboot(lmer4,lmer3))
## can use .progress="text" to get a progress bar ...
boothist <- unlist(boothist)
hist(boothist,breaks=50,col="gray")

obsval <- 2*(logLik(lmer3)-logLik(lmer4))
abline(v=obsval,col=2,lwd=2)
mean(boothist>obsval)
```

```{r pbkrtest,warning=FALSE}
library(pbkrtest)
t1 <- system.time(PBmodcomp(lmer3,lmer4,nsim=500))
```

(takes `r round(t1["elapsed"])` seconds)

In this case, the Kenward-Roger correction appropriately does nothing different -- we have a classical balanced design and no correction is actually necessary.  But it does give a denominator df and a $p$-value for this lmer model, which is handy ...

**Exercise**: repeat for the tundra data

## Predictions and plotting

In `lmer`, `predict` has a `re.form` argument that specifies which random effects should be included (`NA` or `~0`=none, population level; `NULL` (=all) or `~subject`=prediction at the subject level; more complex models, might have additional nested levels).

```{r predictplot}
dataf$pred <- predict(lmer2,re.form=NA)  ## population level
dataf$pred1 <- predict(lmer2) ## individual level
g0 <- ggplot(dataf,aes(mnth,stmass))+
    geom_point()+
    geom_line(aes(group=subject))+
    facet_grid(.~roostsitu)+
    zero_margin
g0 +   geom_line(colour="gray",aes(y=pred1,group=subject)) +
    geom_line(colour="red",aes(y=pred,group=subject))
```

There is so much shrinkage (the among-individual variance is very small) that we can barely see the individual-level predictions (gray lines) behind the population-level predictions (red lines).

Unfortunately computing confidence intervals for the predictions is a little tricky: again, there is some code on the [GLMM FAQ](https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html) for this.

```{r predplot2}
ggplot(dataf,aes(mnth,pred1))+
    geom_line(aes(group=subject,x=as.numeric(mnth)),colour="gray")+
    facet_wrap(~roostsitu,scale="free_y",nrow=1)+
    geom_line(aes(y=pred,x=as.numeric(mnth)),colour="red")
```

For most cases you will want to set up a new data frame to do prediction rather than just using the covariates from the original data (e.g. if the data are sampled irregularly, or sparsely), and use the `newdata` argument of `predict`.  The `expand.grid` function is handy in this context too.

### History

* originally developed for NCEAS summer institute, July 2013
* updated, Alaska ASA workshop, August 2014
* updated, Ottawa workshop, Jan 2016 (removed JAGS material)
* updated, Morelia workshop (moved MCMCglmm to standalone)
