\newcommand{\bm}[1]{#1}

<style>
.refs {
   font-size: 16px;
}
h2 { 
 color: #3399ff;		
}
h3 { 
 color: #3399ff;		
}
.title-slide {
   background-color: #55bbff;
}
</style>
<!--    content: url(https://i.creativecommons.org/l/by-sa/4.0/88x31.png)
>
<!-- Limit image width and height -->
<style type="text/css">
img {     
  max-height: 560px;     
  max-width: 800px; 
}
</style>

```{r knitr,include=FALSE}
require("knitr")
knit_hooks$set(crop=hook_pdfcrop)
##opts_chunk$set(fig.width=5,fig.height=4,
##               out.width="0.6\\textwidth",
##               fig.align="center",echo=FALSE)
opts_chunk$set(echo=FALSE)
```

```{r pkgs,message=FALSE}
library(lattice)
library(plotrix) ## for axis.break
library(reshape2)
library(ggplot2)
library(RColorBrewer)
library(mvbutils) ## for foodweb()
library(lme4)
library(MASS)
library(nlme)
library(grid)
zmargin <- theme(panel.spacing=unit(0,"lines"))
library(scales) ## for 'scientific', 'trans_format'
theme_set(theme_bw())
```

## license

![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.


## (Generalized) linear mixed models

(G)LMMs: a statistical modeling framework incorporating:

-   combinations of categorical and continuous predictors,\
    and interactions
-   (some) non-Normal responses\
    (e.g. binomial, Poisson, and extensions)
-   (some) nonlinearity\
    (e.g. logistic, exponential, hyperbolic)
-   non-independent (grouped) data

---

![image](pix/models-glmm.png)

## Coral protection from seastars (*Culcita*) by symbionts [@mckeon_multiple_2012]

```{r culcita1,cache=TRUE,message=FALSE,warning=FALSE}
## Culcita fits
stierdat <- read.csv("../../data/culcitalogreg.csv")
stierdat$ttt <- as.factor(stierdat$ttt) ## treatment should be a factor
## contrast matrix
cmat <- matrix(c(0,1,1,1,0,1,-1,0,0,1,1,-2),ncol=3,
               dimnames=list(NULL,c("symb","crab.vs.shr","addsymb")))
contrasts(stierdat$ttt) <- cmat
stierdat$block <- as.factor(stierdat$block) ## not really necessary but logical
stierdat <- stierdat[,1:4]   ## don't really need the rest
mod0 <- glm(predation~ttt+block,binomial,data=stierdat)
mod1 <- glm(predation~ttt,binomial,data=stierdat)
library(lme4)
mod2 <- glmer(predation~ttt+(1|block),family=binomial,data=stierdat)
## AGQ 
mod2B <- glmer(predation~ttt+(1|block),family=binomial,
               nAGQ=20,data=stierdat)
mod5=MASS::glmmPQL(predation~ttt,random=~1|block,family=binomial,data=stierdat)
levels(stierdat$ttt) <- c("none","shrimp","crabs","both")
m <- melt(stierdat[,1:3],id.vars=1:2)
m2 <- dcast(m,ttt~variable,fun.aggregate=mean)
m3 <- dcast(m,ttt+block~variable,fun.aggregate=sum)
p <- with(m3,table(predation,ttt))
``` 


```{r culcita2,warning=FALSE,message=FALSE}
load("../../data/culcita.RData")
library(lme4)
cmod1 <- glmer(predation~ttt+(1|block),data=culcita_dat,
               family=binomial,nAGQ=10)
cmod3 <- glm(predation~ttt,data=culcita_dat,family=binomial)
## jump through hoops to get GLM to estimate without complete-separation blocks 7,8,9 messing everything up!
cmod4 <- suppressWarnings(glm(predation~ttt+block,data=transform(culcita_dat,block=relevel(block,2)),
             family=binomial,
             start=c(2,-4,-4,-6,-2,0,0,4,5,30,30,30,5)))
             ## contrasts=list(block=contr.sum),
             ## start=rep(0,)
             ## start=c(10,-4,-5,-6,-20,-8,-8,-8,-5,-3,20,20,20))
r <- ranef(cmod1,condVar=TRUE)
## set up plots
## FIXME:: replace with broom!
library(coefplot2)
f1 <- coeftab(cmod1)
f3 <- coeftab(cmod3)
f4 <- coeftab(cmod4)[1:4,]
r1 <- coeftab(cmod1,ptype="ranef")
r2 <- f3[1,]  ## intercept
r4 <- coeftab(cmod4)[-(1:4),]
v <- data.frame(predict(cmod4,
             newdata=data.frame(block=factor(1:10),ttt="none"),
             se.fit=TRUE))
v$fit <- v$fit-mean(v$fit)
v <- setNames(v[,1:2],c("Estimate","Std..Error"))
w <- data.frame(Estimate=rep(0,10),
                Std..Error=NA,row.names=1:10)
## less clunky way to do this?
dd <- function(X,method,type,
               strip="^ttt") {
  data.frame(method,type,
             p=gsub(strip,"",rownames(X)),X[,1:2],
             stringsAsFactors=FALSE)
}
allEsts <- rbind(dd(f1,"mixed","ttt"),
      dd(f3,"pooled","ttt"),
      dd(f4,"fixed","ttt"),
      dd(r1,"mixed","blocks",
         strip="block\\.\\(Intercept\\)"),
      dd(v,"fixed","blocks"),
      dd(w,"pooled","blocks"))
allEsts <- transform(allEsts,
   p=factor(p,levels=gtools::mixedsort(unique(p))))
```

```{r bplot1,fig.width=8,fig.height=6}
op <- par(las=1,cex=1.5,bty="l")
bb <- barplot(p,ylab="Number of blocks",xlab="Symbionts",
               main="Number of predation events")
bloc <- apply(p,2,cumsum)
midb <- rbind(bloc[1,]/2,
             (bloc[1,]+bloc[2,])/2,
             (bloc[2,]+bloc[3,])/2)
text(bb[col(p)][p>0],midb[p>0],c(0,1,2)[row(p)][p>0])
par(op)
```

## Environmental stress: *Glycera* cell survival (D. Julian unpubl.)

```{r glycera}
x <- read.csv("../../data/Live-Dead Tabulated Counts.csv")
## utility function for factoring/renaming variables before
##  lattice plot
rnfac <- function(dat,vars) {
  if (!all(vars %in% names(dat))) stop("unknown variable")
  for (v in vars) {
    dat[[v]] <- factor(dat[[v]])
    levels(dat[[v]]) <- paste(v,"=",round(as.numeric(levels(dat[[v]])),2),sep="")
  }
  dat
}
sc <- function(x) { (x-min(x))/diff(range(x))}
xsc <- x
predvars <- c("Osm","Cu","H2S","Anoxia")
for (i in predvars) {
  xsc[[i]] <- sc(xsc[[i]])
}
xsc$Osm <- xsc$Osm-0.5
## xsc$O2 <- 1-xsc$O2
## names(xsc)[names(xsc)=="O2"] <- "anox"
xr0 <- within(x,FractionAlive <- Alive/(Alive+Dead))
xr <- melt(subset(xr0,select=-c(Alive,Dead)),id.vars=1:5)

x4 <- dcast(xr,H2S+Anoxia+Cu+Osm~.,fun.aggregate=mean)
names(x4)[5] <- "all"
x5 <- rnfac(x4,c("Anoxia","Osm"))

## FIXME: replace with ColorBrewer colours?
cmgen.colors  <- function (n,h1=6/12,h2=10/12,maxs=0.5)  {
    if ((n <- as.integer(n[1])) > 0) {
        even.n <- n%%2 == 0
        k <- n%/%2
        l1 <- k + 1 - even.n
        l2 <- n - k + even.n
        c(if (l1 > 0) hsv(h = h1,
                          s = seq(maxs, ifelse(even.n, 0.5/k, 0), length.out = l1),
                          v = 1),
          if (l2 > 1) hsv(h = h2,
                          s = seq(0, maxs, length.out = l2)[-1],
                          v = 1))
    }
    else character(0)
}
rb.colors <- function(n) {
  cmgen.colors(n,h1=0,h2=0.7,maxs=1)
}
```
 

```{r glycplot1,out.width="\\textwidth",fig.width=8,fig.height=5}
orig <- trellis.par.get()
pad <- 0 ## 15 for regular layout
trellis.par.set(layout.widths=list(right.padding=pad,left.padding=pad),
                regions=list(col=rb.colors(100)),
##                regions=list(col=brewer.pal(11,"RdBu")),
## leave text alone for regular layout
                add.text=list(cex=0.8),axis.text=list(cex=0.5))
levels(x5$Anoxia) <- c("Normoxia","Anoxia")
## print(levelplot(`(all)`~factor(H2S)*factor(Cu)|Anoxia*Osm,
##          col.region=rb.colors(100),
##          data=x5,
##          xlab=expression(H[2]*S),
##          ylab="Copper"))
levelplot(all~factor(H2S)*factor(Cu)|Osm*Anoxia,
                col.region=rb.colors(100), ## brewer.pal(11,"RdBu"), ## rb.colors(100),
                data=x5,
                xlab=expression(H[2]*S),
                ylab="Copper")
trellis.par.set(theme=orig) ## restore settings
## FIXME: redo in ggplot2?  LOW PRIORITY
```

## *Arabidopsis* response to fertilization & herbivory [@banta_comprehensive_2010]

```{r arabplot1,fig.width=7,fig.height=6}
panel.stripplot2 <-
function (x, y, jitter.data = FALSE, factor = 0.5, amount = NULL, 
    horizontal = TRUE, groups = NULL, ...) 
{
    if (!any(is.finite(x) & is.finite(y))) 
        return()
    panel.sizeplot(x = x, y = y, jitter.x = jitter.data && !horizontal, 
        jitter.y = jitter.data && horizontal, factor = factor, 
        amount = amount, groups = groups, horizontal = horizontal, 
        ...)
}
load("../../data/Banta.RData")
trellis.par.set(list(fontsize=list(text=20)))
stripplot(jltf ~ amd|nutrient, 
                data=within(dat.tf,jltf <-jitter(log(total.fruits+1),
                  amount=0.05)),
                strip=strip.custom(strip.names=c(TRUE,TRUE)),
                groups=gen, type=c('p','a'),
                ylab="Log(1+fruit set)")
##                main="panel: nutrient, color: genotype")
## trellis.par.set(theme=orig) ## restore settings
```

## Coral demography (J.-S. White unpubl.)

```{r coral_demog,warning=FALSE}
L <- load("../../data/m.acr.jagsout.RData")
L2 <- load("../../data/m.acr.lme4out.RData")
L3 <- load("../../data/demog.mort.18apr.RData")
source("../../R/demog_mort_funs.R")
plotfun(j.red2,m.acr.nofr,drop.cols=c(4,8))
```

## Simple mixed models<br>(scalar, intercept-only)

$$
\begin{split}
y_{ij} & = \beta_0 + \beta_1 x_{ij} + \epsilon_{0,ij} + \epsilon_{1,j} \\
& = (\beta_0 + \epsilon_{1,j}) + \beta_1 x_{ij} + \epsilon_{1,j} \\
\epsilon_{0,ij} & \sim \textrm{Normal}(0,\sigma_0^2) \\
\epsilon_{1,j} & \sim \textrm{Normal}(0,\sigma_1^2)
\end{split}
$$

- Could have multiple, nested levels of random effects  
(genotype within population within region ...), or *crossed* REs

## Random-slopes model

$$
\begin{split}
y_{ij} & = \beta_0 + \beta_1 x_{ij} + \epsilon_{0,ij} + \epsilon_{1,j} +
\epsilon_{2,j} x_{ij} \\
& = (\beta_0 + \epsilon_{1,j}) + (\beta_1 + \epsilon_{2,j}) x_{ij})  \\
\epsilon_{0,ij} & \sim \textrm{Normal}(0,\sigma_0^2) \\
\{\epsilon_{1,j}, \epsilon_{2,j}\} & \sim \textrm{MVN}(0,\Sigma)
\end{split}
$$

- variation in the *effect* of a treatment or covariate across groups
- variation in intercept, slope allowed to be correlated

## General definition

$$
\begin{split}
\underbrace{Y_i}_{\text{response}} & \sim \overbrace{\text{Distr}}^{\substack{\text{conditional} \\ \text{distribution}}}(\underbrace{g^{-1}(\eta_i)}_{\substack{\text{inverse} \\ \text{link} \\ \text{function}}},\underbrace{\phi}_{\substack{\text{scale} \\ \text{parameter}}}) \\
\underbrace{\boldsymbol \eta}_{\substack{\text{linear} \\ \text{predictor}}} & 
 = 
\underbrace{\boldsymbol X \boldsymbol \beta}_{\substack{\text{fixed} \\ \text{effects}}} + 
\underbrace{\boldsymbol Z \boldsymbol b}_{\substack{\text{random} \\ \text{effects}}}
\\
\underbrace{\boldsymbol b}_{\substack{\text{conditional} \\ \text{modes}}}  & 
\sim \text{MVN}(\boldsymbol 0, \underbrace{\Sigma(\boldsymbol \theta)}_{\substack{\text{variance-} \\ \text{covariance} \\ \text{matrix}}})
\end{split}
$$

## What are random effects?

A method for …

-   accounting for among-individual, within-block correlation
-   compromising between\
    *complete pooling* (no among-block variance)\
     and *fixed effects* (large among-block variance)
-   handling levels selected at random from a larger population
-   sharing information among levels (*shrinkage
    estimation*)
-   estimating variability among levels
-   allowing predictions for unmeasured levels

## Random-effect myths

-   levels of random effects must always be sampled at random
-   a complete sample cannot be treated as a random effect
-   random effects are always a *nuisance variable*
-   nothing can be said about the predictions of a random effect
-   you should always use a random effect no matter how few levels you
have

## Reasons for random effects (inferential/philosophical)

- **don't** want to
     - test hypotheses about differences between responses at particular levels of
the grouping variable;
- **do** want to
     - quantify variation among groups
     - make predictions about unobserved groups
- have levels that are randomly sampled from/representative of a larger population (*exchangeable*)


## Reasons for random effects (practical)

- want to combine information across groups
- have variation in information per level (number of samples or noisiness);
- have a categorical predictor that is a nuisance variable (i.e., it is not of direct interest, but should be controlled for).
- have more than 5-6 groups, or regularizing/using priors (otherwise, use fixed)

See also @Crawley2002; @gelman_analysis_2005

## Avoiding MM

* average: for *nested* LMMs (when uninterested in within-group variation) [@murtaugh_simplicity_2007] 
* use fixed effects (or two-stage models) instead of random effects when
     * many samples per block (shrinkage unimportant)
     * few blocks (little advantage; bad variance estimates)

Estimation
==========

Overview
--------

## Maximum likelihood estimation

-   Best fit is a compromise between two components\
    (consistency of data with fixed effects and conditional modes;
    consistency of random effect with RE distribution)
	-   Goodness-of-fit *integrates* over conditional modes

## ...

```{r plotex,message=FALSE}
set.seed(101)
dd <- data.frame(f=gl(5,5))
dd$y <- simulate(~1+(1|f),newdata=dd,
                 family=gaussian,seed=101,
                 newparams=list(theta=1,beta=0,sigma=1))[[1]]
ggplot(dd,aes(x=f,y=y))+geom_point()+
    stat_summary(fun.y=mean,geom="point",size=3,colour="blue",
                 pch=3)+
     geom_point(data=subset(dd,y<(-2)),colour="red",size=2)+
         theme_update(panel.grid.major=element_blank(),
                      panel.grid.minor=element_blank())
```

## Shrinkage: *Arabidopsis* conditional modes

```{r arabshrink,fig.height=6,fig.width=8}
z<- subset(dat.tf,amd=="clipped" & nutrient=="1")
m1 <- glm(total.fruits~gen-1,data=z,family="poisson")
m2 <- glmer(total.fruits~1+(1|gen),data=z,family="poisson")
tt <- table(z$gen)
rr <- unlist(ranef(m2)$gen)[order(coef(m1))]+fixef(m2)
m1s <- sort(coef(m1))
m1s[1:2] <- rep(-5,2)
gsd <- attr(VarCorr(m2)$gen,"stddev")
gm <- fixef(m2)
nseq <- seq(-3,6,length.out=50)
sizefun <- function(x,smin=0.5,smax=3,pow=2) {
    smin+(smax-smin)*((x-min(x))/diff(range(x)))^pow
}
nv <- dnorm(nseq,mean=gm,sd=gsd)
##
op <- par(las=1,cex=1.5,bty="l")
plot(exp(m1s),xlab="Genotype",ylab="Mean fruit set",
     axes=FALSE,xlim=c(-0.5,25),log="y",yaxs="i",xpd=NA,
     pch=16,cex=0.5)
axis(side=1)
axis(side=2,at=c(exp(-5),0.1,1,10,20),
     labels=c(0,0.1,1,10,20),cex=0.8)
##     ylim=c(-3,5))
polygon(c(rep(0,50),nv*10),exp(c(rev(nseq),nseq)),col="gray",xpd=NA)
n <- tt[order(coef(m1))]
points(exp(rr),pch=16,col=adjustcolor("red",alpha=0.5),
       cex=sizefun(n),xpd=NA)
## text(seq_along(rr),rr,n,pos=3,xpd=NA,cex=0.6)
box()
axis.break(axis=2,breakpos=exp(-4))
legend("bottomright",
       c("group mean","shrinkage est."),
       pch=16,pt.cex=c(1,2),
       col=c("black",adjustcolor("red",alpha=0.5)),
       bty="n")
par(op)
```

## Shrinkage in a random-slopes model

From Christophe Lalanne, see [here](https://stats.stackexchange.com/questions/51186/what-would-be-an-illustrative-picture-for-linear-mixed-models):

```{r sleepstudy_shrinkage, echo=FALSE}
library(lme4)
data(sleepstudy)

## Fit individual regression lines for each subject
dfrm <- coef(lmList(Reaction ~ Days | Subject, sleepstudy))

## Estimate parameters of a random intercept and random intercept and slope model
m1 <- lmer(Reaction ~ Days + (1 | Subject), data=sleepstudy)
m2 <- lmer(Reaction ~ Days + (Days | Subject), data=sleepstudy)

## Put all estimates (intercept + slope for each model) into the same data.frame
dfrm <- cbind.data.frame(dfrm,
                         as.data.frame(coef(m1)[["Subject"]]),
                         as.data.frame(coef(m2)[["Subject"]]))

## Kernel density estimates for the distribution of individual intercepts
intcpt.dens <- list()
idx <- seq(1, ncol(dfrm), by=2)
for (i in seq_along(idx))
  intcpt.dens[[i]] <- density(as.numeric(dfrm[,idx[i]]), adj=1.4)
len <- length(intcpt.dens[[1]]$x)

## Show all
cols <- c("grey30", "#D95F02", "#669999")
xyplot(Reaction ~ Days, data=sleepstudy,
       xlim=c(0, 8), ylim=c(150, 450), ylab="Fitted reaction time",
       scales=list(x=list(at=seq(0, 8, by=1))),
       key=list(corner=c(0,1), text=list(c("within-group",
                                 "random intercept",
                                 "random intercept and slope"),
                                 col=cols, cex=0.8)),
       panel=function(...) {
         apply(dfrm[,1:2], 1, panel.abline, col=cols[1], alpha=.5, lwd=1.2)
         apply(dfrm[,3:4], 1, panel.abline, col=cols[2], alpha=.5, lwd=1.2)
         for (i in seq_along(idx))
         panel.lines(x=c(intcpt.dens[[i]]$y*100, rep(0, len)),
                     y=c(intcpt.dens[[i]]$x, rev(intcpt.dens[[i]]$x)), col=cols[i], lwd=1.8)
})
```

Methods
-------

## Estimation methods

- deterministic
    - various approximate integrals [@breslow_whither_2004]
    - penalized quasi-likelihood, Laplace, Gauss-Hermite quadrature, … [@biswas2015];  
best methods needed for large variance, small clusters
    -   flexibility and speed vs. accuracy
- stochastic
- stochastic (Monte Carlo): frequentist and Bayesian
    - [@booth_maximizing_1999; @sung_monte_2007; @ponciano_hierarchical_2009]
    -   usually slower but flexible and accurate

## Estimation: *Culcita* [@mckeon_multiple_2012]

```{r stierfig2,warning=FALSE,fig.width=8,fig.height=8}
library(coefplot2)
cvec <- c("purple","magenta","blue","black","darkgreen")
mnames <- c("GLM (fixed)","GLM (pooled)","PQL","Laplace","AGQ")
coefplot2(list(mod2B,mod2,mod5,mod1,
               coeftab(mod0)[1:4,]),
          col.pts=cvec,
          varnames=c("Symbiont",
            "Crab vs. Shrimp",
            "Added symbiont"),xlim=c(-7,3),
          ylim=c(0.9,3.7),
          spacing=0.14,
          main="Log-odds of predation")
par(xpd=NA)
text(rep(1,5),seq(1.6,by=-0.18,length.out=5),
     mnames,col=rev(cvec),cex=0.7,adj=0)
```

# Inference

## Wald tests

-   typical results of `summary()`
-   exact for ANOVA, regression:  
approximation for GLM(M)s
-   fast
-   approximation is sometimes awful (Hauck-Donner effect)

## Likelihood ratio tests

-   better than Wald, but still have two problems:
-   “denominator degrees of freedom”  
(when estimating scale)
    -   finite-size issues for GLM(M)s: Bartlett corrections
    -   Kenward-Roger correction? [@stroup_rethinking_2014]
-   Profile confidence intervals: expensive/fragile

## finite-size p-values??

- hope that min grps > 42
- guess denominator DF from classic design ([R code](../../R/calcDenDF.R))
- conservative: take minimum number of groups - 1
- Satterthwaite/Kenward-Roger (`lmerTest`, LMMs only)
- parametric bootstrap (`pbkrtest`)

## Parametric bootstrapping

-   fit null model to data
-   simulate “data” from null model
-   fit null and working model, compute likelihood difference
-   repeat to estimate null distribution
-   should be OK but ??? not well tested  
(assumes estimated parameters are "sufficiently" good)

## Bayesian inference

-   If we have a good sample from the posterior distribution (Markov chains have converged etc. etc.) we get most of the inferences we want for free by summarizing the marginal posteriors
-   *post hoc* Bayesian methods: use deterministic/frequentist methods to find the maximum, then sample around it

# Challenges & open questions

## On beyond `lme4`

-  `glmmTMB`: zero-inflated and other distributions
-  `brms`,`rstanarm`: interfaces to Stan
-  `INLA`: spatial and temporal correlations
- [rethinking package](https://github.com/rmcelreath/rethinking)

## Toolboxes

- JAGS (R: `rjags`, `r2jags`)
- TensorFlow (R: `greta`)
- NIMBLE (R: `nimble` package)
- Stan (R: `rstan`)
- TMB (R: `TMB`)

## On beyond R

-   Julia: `MixedModels.jl` package
-   SAS: PROC MIXED, NLMIXED
-   AS-REML
-   Stata (GLLAMM, xtmelogit)
-   HLM, MLWiN

![image](pix/300px-On_Beyond_Zebra.jpg){width="3cm"}

## Challenges

-   Small clusters: need AGQ/MCMC
-   Small numbers of clusters: need finite-size corrections (KR/PB/MCMC)
-   Small data sets: issues with *singular* fits  
[@barr_random_2013] vs. [@bates_parsimonious_2015]
-   Big data: speed!
-   Model diagnosis
-   Confidence intervals accounting for uncertainty in variances

## The need for speed [@brooks_modeling_2017]

<!-- https://tex.stackexchange.com/questions/250623/how-to-generate-a-4-up-table-of-images-with-pandoc -->


![](pix/bigfit2.png)

## more speed

![](pix/bigfit1.png)
 

## Spatial and temporal correlations

-   Sometimes blocking takes care of non-independence ...
-   but sometimes there is temporal or spatial correlation *within* blocks
-    ... also phylogenetic ... [@ives_statistics_2006]
-   "G-side" vs. "R-side" effects
-   tricky to implement for GLMMs, but new possibilities on the horizon
[@Rue+2009; @rousset_testing_2014];
[https://github.com/stevencarlislewalker/lme4ord](lme4ord package)

## Next steps

-   Complex random effects:   
    regularization, model selection, penalized methods (lasso/fence)
-   Flexible correlation and variance structures
-   Flexible/nonparametric random effects distributions
-   hybrid & improved MCMC methods
-   *Reliable* assessment of out-of-sample performance

## resources

- <http://ms.mcmaster.ca/~bolker/misc/private/14-Fox-Chap13.pdf>
- <https://bbolker.github.io/mixedmodels-misc/ecostats_chap.html>
-  @bolker_glmm_2014

<img src="pix/foxbook.png" height="100">

(code ASPROMP8)

## references  {.refs}
