---
title: "Post-fitting inference, etc. with glmmTMB models"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{model_evaluation}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

The purpose of this vignette is to describe (and test) the
functions in various downstream packages that are available for summarizing
and otherwise interpreting `glmmTMB` fits.
Some of the packages/functions discussed below may
not be suitable for inference on parameters of
the zero-inflation or dispersion
models, but will be restricted to the conditional-mean model.

```{r pkgs,message=FALSE}
library(glmmTMB)
library(car)
library(emmeans)
library(effects)
library(multcomp)
library(MuMIn)
library(DHARMa)
## library(broom)
## library(broom.mixed)
library(dotwhisker)
library(ggplot2); theme_set(theme_bw())
## retrieve slow stuff
L <- load(system.file("vignette_data","model_evaluation.rda",package="glmmTMB"))
```

A couple of example models:

```{r fit_model,eval=FALSE}
owls_nb1 <- glmmTMB(SiblingNegotiation ~ FoodTreatment*SexParent +
                        (1|Nest)+offset(log(BroodSize)),
                    contrasts=list(FoodTreatment="contr.sum",
                                   SexParent="contr.sum"),
                    family = nbinom1,
                    zi = ~1, data=Owls)
```

```{r fit_model2,cache=TRUE}
data("cbpp",package="lme4")
cbpp_b1 <- glmmTMB(incidence/size~period+(1|herd),
                   weights=size,family=binomial,
                   data=cbpp)
## simulated three-term Beta example
set.seed(1001)
dd <- data.frame(z=rbeta(1000,shape1=2,shape2=3),
                 a=rnorm(1000),b=rnorm(1000),c=rnorm(1000))
simex_b1 <- glmmTMB(z~a*b*c,family=beta_family,data=dd)
```

## model checking and diagnostics

### DHARMa

The `DHARMa` package provides diagnostics for hierarchical models.

After running
```{r dharma_fake,eval=FALSE}
owls_nb1_simres <- simulateResiduals(owls_nb1)
```
```{r dharma,cache=TRUE}
system.time(sr <- simulateResiduals(owls_nb1))
```

```{r dharmaplot,fig.width=8,fig.height=4}
plot(sr)
```

#### issues

See warning. `DHARMa` Will only work for models using families for which a `simulate` method has been implemented (in `TMB`, and appropriately reflected in `glmmTMB`).

## Inference

### car::Anova

We can use `car::Anova()` to get traditional ANOVA-style tables from `glmmTMB` fits. A few limitations/reminders:

- these tables are using Wald $\chi^2$ statistics for comparisons (neither likelihood ratio tests nor $F$ tests)
- they apply to the fixed effects of the conditional component of the model only (other components *might* work, but haven't been tested at all)
- as always, if you want to do type 3 tests, you should probably set sum-to-zero contrasts on factors and center numerical covariates (see `contrasts` argument above)

```{r Anova}
Anova(owls_nb1)  ## default type II
Anova(owls_nb1,type="III")
```

### effects

```{r effects,fig.width=8,fig.height=4}
(ae <- allEffects(owls_nb1))
plot(ae)
plot(allEffects(simex_b1))
```

### emmeans

```{r emmeans}
emmeans(owls_nb1, poly ~ FoodTreatment | SexParent)
```

### drop1

`stats::drop1` is a built-in R function that refits the model with various terms dropped. In its default mode it respects marginality (i.e., it will only drop the top-level interactions, not the main effects):

```{r drop1,cache=TRUE}
system.time(owls_nb1_d1 <- drop1(owls_nb1,test="Chisq"))
```

```{r print_drop1}
print(owls_nb1_d1)
```

In principle, using `scope = . ~ . - (1|Nest)` should work to execute a "type-3-like" series of tests, dropping the main effects one at a time while leaving the interaction in (we have to use `- (1|Nest)` to exclude the random effects because `drop1` can't handle them). However, due to the way that R handles formulas, dropping main effects from an interaction of *factors* has no effect on the overall model. (It would work if we were testing the interaction of continuous variables.)

#### issues

The `mixed` package implements a "true" type-3-like parameter-dropping mechanism for `[g]lmer` models. Something like that could in principle be applied here.

### MuMIn

Model selection and averaging.

We can run `MuMIn::dredge(owls_nb1)` on the model to fit all possible submodels.
Since this takes a little while (45 seconds or so), we'll instead load some previously computed results:

```{r dredge1}
owls_nb1_dredge
```

```{r plot_dredge1,fig.width=8,fig.height=8}
op <- par(mar=c(2,5,14,3))
plot(owls_nb1_dredge)
par(op) ## restore graphics parameters
```

Model averaging:

```{r model_avg}
model.avg(owls_nb1_dredge)
```

#### issues

- may not work for Beta models because the `$family` component ("beta") is not identical to the name of the family function (`beta_family()`)?
(Kamil Bartoń, pers. comm.)

### multcomp

Multiple comparisons and *post hoc* tests.

```{r glht_def}
glht_glmmTMB <- function (model, ..., component="cond") {
    glht(model, ...,
         coef. = function(x) fixef(x)[[component]],
         vcov. = function(x) vcov(x)[[component]],
         df = NULL)
}
modelparm.glmmTMB <- function (model, coef. = function(x) fixef(x)[[component]],
                               vcov. = function(x) vcov(x)[[component]],
                               df = NULL, component="cond", ...) {
    multcomp:::modelparm.default(model, coef. = coef., vcov. = vcov.,
                        df = df, ...)
}
```

```{r glht_ex}
g1 <- glht(cbpp_b1, linfct = mcp(period = "Tukey"))
summary(g1)
```

#### issues

It is possible to make `multcomp` work in a way that (1) actually
uses the S3 method structure and (2) doesn't need access to
private `multcomp` methods (i.e. accessed by `multcomp:::`) ? Not sure,
but both of the following hacks should work. (The `glht_glmmTMB` solution below is clunky because it isn't a real S3 method; the `model.parm.glmmTMB` solution can't be included in the package source code as-is because `:::` is not allowed in CRAN package code.)

### broom etc.

The `broom`, `broom.mixed` packages are designed to extract information from a broad range of models in a convenient (tidy) format; the `dotwhisker` package builds on this platform to draw elegant coefficient plots.

```{r dwplot,fig.width=8, eval=FALSE}
(t1 <- broom.mixed::tidy(owls_nb1, conf.int = TRUE))
if (packageVersion("dotwhisker")>"0.4.1") {
    ## to get this version (which fixes various dotwhisker problems)
    ## use devtools::install_github("bbolker/broom.mixed") or
    ## wait for pull request acceptance/submission to CRAN/etc.
    dwplot(owls_nb1)+geom_vline(xintercept=0,lty=2)
} else {
    owls_nb1$coefficients <- TRUE  ## hack!
    dwplot(owls_nb1,by_2sd=FALSE)+geom_vline(xintercept=0,lty=2)
}
```

#### issues

(these are more general `dwplot` issues)

- use black rather than color(1) when there's only a single model, i.e. only add `aes(colour=model)` conditionally?
- draw points even if std err / confint are `NA` (draw `geom_point()` as well as `geom_pointrange()`? need to apply all aesthetics, dodging, etc. to both ...)
- for `glmmTMB` models, allow labeling by component? or should this be done by manipulating the tidied frame first? (i.e. `tidy(.) %>% tidyr::unite(term,c(component,term))`)

```{r save_output,echo=FALSE,eval=FALSE}
## store time-consuming stuff
save("owls_nb1",
     "owls_nb1_simres",
     "owls_nb1_dredge",
     file="../inst/vignette_data/model_evaluation.rda")
```

## to do

- add: `texreg` (regression tables)? `sjPlot` ? `memisc`, `huxtable`, ... ?
