Covariance structures with glmmTMB

The AR(1) covariance structure

n <- 25                                              ## Number of time points
x <- MASS::mvrnorm(mu = rep(0,n),
             Sigma = .7 ^ as.matrix(dist(1:n)) )    ## Simulate the process using the MASS package
y <- x + rnorm(n)                                   ## Add measurement noise

times <- factor(1:n, levels=1:n)
head(levels(times))

group <- factor(rep(1,n))

dat0 <- data.frame(y, times, group)

glmmTMB(y ~ ar1(times + 0 | group), data=dat0)

simGroup <- function(g, n=6, phi=0.7) {
    x <- MASS::mvrnorm(mu = rep(0,n),
             Sigma = phi ^ as.matrix(dist(1:n)) )   ## Simulate the process
    y <- x + rnorm(n)                               ## Add measurement noise
    times <- factor(1:n)
    group <- factor(rep(g,n))
    data.frame(y, times, group)
}
simGroup(1)

dat1 <- do.call("rbind", lapply(1:1000, simGroup) )

(fit.ar1 <- glmmTMB(y ~ ar1(times + 0 | group), data=dat1))

The unstructured covariance

We can try to fit an unstructured covariance to the previous dataset dat. For this case an unstructured covariance has 300 correlation parameters and 25 variance parameters. Adding \(\sigma_0^2 I\) on top would cause a strict overparameterization, as these would be redundant with the diagonal elements in the covariance matrix. Hence, when fitting the model with glmmTMB, we have to disable the \(\varepsilon\) term (the dispersion) by setting dispformula=~0:

fit.us <- glmmTMB(y ~ us(times + 0 | group), data=dat1, dispformula=~0)
fit.us$sdr$pdHess ## Converged ?

The estimated variance and correlation parameters are:

VarCorr(fit.us)

The estimated correlation is approximately constant along diagonals (apparent Toeplitz structure) and we note that the first off-diagonal is now ca. half the true value (0.7) because the dispersion is effectively included in the estimated covariance matrix (i.e. \(\rho' = \rho {\sigma^2_{{\text {AR}}}}/({\sigma^2_{{\text {AR}}}} + {\sigma^2_{{\text {meas}}}})\)).

The Toeplitz structure

The next natural step would be to reduce the number of parameters by collecting correlation parameters within the same off-diagonal. This amounts to 24 correlation parameters and 25 variance parameters.

We use dispformula = ~0 to suppress the residual variance (it actually gets set to a small value controlled by the zerodisp_val argument of glmmTMBControl())¹

fit.toep <- glmmTMB(y ~ toep(times + 0 | group), data=dat1,
                    dispformula=~0)
fit.toep$sdr$pdHess ## Converged ?

The estimated variance and correlation parameters are:

(vc.toep <- VarCorr(fit.toep))

The diagonal elements are all approximately equal to the true total variance (\({\sigma^2_{{\text {AR}}}} + {\sigma^2_{{\text {meas}}}}\)=2), and the off-diagonal elements are approximately equal to the expected value of 0.7/2=0.35.

vc1 <- vc.toep$cond[[1]] ## first term of var-cov for RE of conditional model
summary(diag(vc1))
summary(vc1[row(vc1)!=col(vc1)])

We can get a slightly better estimate of the variance by using REML estimation (however, the estimate of the correlations seems to have gotten slightly worse):

fit.toep.reml <- update(fit.toep, REML=TRUE)
vc1R <- VarCorr(fit.toep.reml)$cond[[1]]
summary(diag(vc1R))
summary(vc1R[row(vc1R)!=col(vc1R)])

Compound symmetry

The compound symmetry structure collects all off-diagonal elements of the correlation matrix to one common value.

We again use dispformula = ~0 to make the model parameters identifiable (see the footnote in The Toeplitz structure; a similar, although slightly simpler, argument applies here).

fit.cs <- glmmTMB(y ~ cs(times + 0 | group), data=dat1, dispformula=~0)
fit.cs$sdr$pdHess ## Converged ?

The estimated variance and correlation parameters are:

VarCorr(fit.cs)

Anova tables

The models ar1, toep, and us are nested so we can use:

anova(fit.ar1, fit.toep, fit.us)

ar1 has the lowest AIC (it's the simplest model, and fits the data adequately); we can't reject the (true in this case!) null model that an AR1 structure is adequate to describe the data.

The model cs is a sub-model of toep:

anova(fit.cs, fit.toep)

Here we can reject the null hypothesis of compound symmetry (i.e., that all the pairwise correlations are the same).

Adding coordinate information

Coordinate information can be added to a variable using the glmmTMB function numFactor. This is necessary in order to use those covariance structures that require coordinates. For example, if we have the numeric coordinates

x <- sample(1:2, 10, replace=TRUE)
y <- sample(1:2, 10, replace=TRUE)

we can generate a factor representing \((x,y)\) coordinates by

(pos <- numFactor(x,y))

Numeric coordinates can be recovered from the factor levels:

parseNumLevels(levels(pos))

In order to try the remaining structures on our test data we re-interpret the time factor using numFactor:

dat1$times <- numFactor(dat1$times)
levels(dat1$times)

Ornstein–Uhlenbeck

Having the numeric times encoded in the factor levels we can now try the Ornstein–Uhlenbeck covariance structure.

fit.ou <- glmmTMB(y ~ ou(times + 0 | group), data=dat1)
fit.ou$sdr$pdHess ## Converged ?

It should give the exact same results as ar1 in this case since the times are equidistant:

VarCorr(fit.ou)

However, note the differences between ou and ar1:

• ou can handle irregular time points.
• ou only allows positive correlation between neighboring time points.

Spatial correlations

The structures exp, gau and mat are meant to used for spatial data. They all require a Euclidean distance matrix which is calculated internally based on the coordinates. Here, we will try these models on the simulated time series data.

An example with spatial data is presented in a later section.

fit.mat <- glmmTMB(y ~ mat(times + 0 | group), data=dat1, dispformula=~0)
fit.mat$sdr$pdHess ## Converged ?

VarCorr(fit.mat)

fit.gau <- glmmTMB(y ~ gau(times + 0 | group), data=dat1, dispformula=~0)
fit.gau$sdr$pdHess ## Converged ?

VarCorr(fit.gau)

fit.exp <- glmmTMB(y ~ exp(times + 0 | group), data=dat1)
fit.exp$sdr$pdHess ## Converged ?

VarCorr(fit.exp)

d <- data.frame(z = as.vector(volcano),
                x = as.vector(row(volcano)),
                y = as.vector(col(volcano)))

set.seed(1)
d$z <- d$z + rnorm(length(volcano), sd=15)
d <- d[sample(nrow(d), 100), ]

volcano.data <- array(NA, dim(volcano))
volcano.data[cbind(d$x, d$y)] <- d$z
image(volcano.data, main="Spatial data", useRaster=TRUE)

d$pos <- numFactor(d$x, d$y)
d$group <- factor(rep(1, nrow(d)))

f <- glmmTMB(z ~ 1 + exp(pos + 0 | group), data=d)

confint(f, "sigma")

newdata <- data.frame( pos=numFactor(expand.grid(x=1:3,y=1:3)) )
newdata$group <- factor(rep(1, nrow(newdata)))
newdata

predict(f, newdata, type="response", allow.new.levels=TRUE)

predict_col <- function(i) {
    newdata <- data.frame( pos = numFactor(expand.grid(1:87,i)))
    newdata$group <- factor(rep(1,nrow(newdata)))
    predict(f, newdata=newdata, type="response", allow.new.levels=TRUE)
}

pred <- sapply(1:61, predict_col)

image(pred, main="Reconstruction")

Mappings

For various advanced purposes, such as computing likelihood profiles, it is useful to know the details of the parameterization of the models - the scale on which the parameters are defined (e.g. standard deviation, variance, or log-standard deviation for variance parameters) and their order.

vv0 <- VarCorr(fit.us)
vv1 <- vv0$cond$group          ## extract 'naked' V-C matrix
n <- nrow(vv1)
rpars <- getME(fit.us,"theta") ## extract V-C parameters
## first n parameters are log-std devs:
all.equal(unname(diag(vv1)),exp(rpars[1:n])^2)
## now try correlation parameters:
cpars <- rpars[-(1:n)]
length(cpars)==n*(n-1)/2      ## the expected number
cc <- diag(n)
cc[upper.tri(cc)] <- cpars
L <- crossprod(cc)
D <- diag(1/sqrt(diag(L)))
round(D %*% L %*% D,3)
round(unname(attr(vv1,"correlation")),3)

all.equal(c(cov2cor(vv1)),c(fit.us$obj$env$report(fit.us$fit$parfull)$corr[[1]]))

## want $par, not $parfull: do NOT include conditional modes/'b' parameters
ppar <- fit.us$fit$par
length(ppar)
range(which(names(ppar)=="theta")) ## the last n*(n+1)/2 parameters
## only 1 fixed effect parameter
tt <- tmbprofile(fit.us$obj,2,trace=FALSE)

confint(tt)
plot(tt)

ppar <- fit.cs$fit$par
length(ppar)
range(which(names(ppar)=="theta")) ## the last n*(n+1)/2 parameters
## only 1 fixed effect parameter, 1 dispersion parameter
tt2 <- tmbprofile(fit.cs$obj,3,trace=FALSE)

plot(tt2)

Generalized latent variable model

Consider a generalized linear mixed model

where \(g(.)\) is the link function; \(\boldsymbol{\beta}\) is a p-dimensional vector of regression coefficients related to the covariates; \(\boldsymbol{X}\) is an \(n \times p\) model matrix; and \(\boldsymbol{Z}\) is the \(n\times q\) model matrix for the \(q\)-dimensional vector-valued random effects variable \(\boldsymbol{U}\) which is multivariate normal with mean zero and a parameterized \(q \times q\) variance-covariance matrix, \(\boldsymbol{\Sigma}\), i.e., \(\boldsymbol{U} \sim N(\boldsymbol{0}, \boldsymbol{\Sigma})\).

where \(\otimes\) is the Kronecker product and \(\boldsymbol{\Lambda} = (\boldsymbol{\lambda_1}, \ldots, \boldsymbol{\lambda_d})'\) is the \(q \times d\) matrix of factor loadings (with \(d \ll q\)). The upper triangular elements of \(\boldsymbol{\Lambda}\) are set to be zero to ensure parameter identifiability. Here we assume that the latent variables follow a multivariate standard normal distribution, \(\boldsymbol{b} \sim N(\boldsymbol{0}, \boldsymbol{I})\).

For GLVMs it is important to select initial starting values for the parameters because the observed likelihood may be multimodal, and maximization algorithms can end up in local maxima. Niku et al. (2019) describe methods to enable faster and more reliable fits of latent variable models by carefully choosing starting values of the parameters.

A similar method has been implemented in glmmTMB. A generalized linear model is fitted to the data to obtain initial starting values for the fixed parameters in the model. Residuals from the fitted GLM are calculated; Dunn-Smyth residuals are calculated for common families while residuals from the dev.resids() function are used otherwise. Initial starting values for the latent variables and their loadings are obtained by fitting a reduced rank model to the residuals.

## fit rank-reduced models with varying dimension
dvec <- 2:10
fit_list <- lapply(dvec,
                   function(d) {
                       glmmTMB(abund ~ Species + rr(Species + 0|id, d = d),
                               data = spider_long)
                   })
names(fit_list) <- dvec
## compare fits via AIC
aic_vec <- sapply(fit_list, AIC)
delta_aic  <- aic_vec - min(aic_vec, na.rm = TRUE)

spider_rr <- glmmTMB(abund ~ Species + rr(Species + 0|id, d = 3),
                     data = spider_long)
re <- as.data.frame(ranef(spider_rr))
re <- within(re, {
    ## sites in numeric order
    grp <- factor(grp, levels = unique(grp))
    ## species in site-1-predicted-abundance order
    term <- reorder(term, condval, function(x) x[1])
    lwr <- condval - 2*condsd
    upr <- condval + 2*condsd
})
if (require("ggplot2")) {
    ggplot(re, aes(grp, condval)) +
        geom_pointrange(aes(ymin=lwr, ymax = upr)) +
        facet_wrap(~term, scale = "free")
}

source(system.file("misc", "extract_rr.R", package = "glmmTMB"))
rr_info <- extract_rr(spider_rr)
lapply(rr_info, dim)

par(las = 1)
afac <- 4
sp_names <- abbreviate(gsub("Species", "", rownames(rr_info$fl)))
plot(rr_info$fl[,1], rr_info$fl[,2], xlab = "factor 1", ylab = "factor 2", pch = 16, cex = 2)
text(rr_info$b[,1]*afac*1.05, rr_info$b[,2]*afac*1.05, rownames(rr_info$b))
arrows(0, 0, rr_info$b[,1]*afac, rr_info$b[,2]*afac)
text(rr_info$fl[,1], rr_info$fl[,2], sp_names, pos = 3, col = 2)

Proptional

The random effect structure propto fits multivariate random effects proportional to a known variance-covariance matrix. One way the propto structure can be used is in phylogenetic analysis where a random effect proportional to a phylogenetic variance-covariance matrix aims to account for the correlation across species due to their shared ancestry. For example, the carni70 data set from the ade4 package describes the phylogeny along with the geographic range and body size of 70 carnivora. To account for the dependence among species due to shared ancestral history we can include a phylogenetically structured error term in the model via the propto structure as follows:

require(ade4)
require(ape)
data(carni70)
carnidat <- data.frame(species = rownames(carni70$tab), carni70$tab)
tree <- read.tree(text=carni70$tre)
phylo_varcov <- vcv(tree)# phylogenetic variance-covariance matrix
# row/column names of phylo_varcov must match factor levels in data
rownames(phylo_varcov) <- colnames(rownames(phylo_varcov)) <- gsub(".", "_", rownames(phylo_varcov)) 
carnidat$dummy <- factor(1) # a dummy grouping variable must be added to the dataset

fit_phylo <- glmmTMB(log(range) ~ log(size) + propto(0 + species | dummy, phylo_varcov),
                     data = carnidat)

dummy is a dummy variable equal to one for all observations to specify that all observations belong to the same cluster. The intercept term is excluded from the proportional random effect -- this is to ensure that each random effect corresponds to the effect for its corresponding species. It is important that the row/columns names of the matrix match the terms in the random effect (see Construction of structured covariance matrices for how the terms are constructed).

Construction of structured covariance matrices

This section will explain how covariance matrices are constructed "under the hood", and in particular why the 0+ term is generally required in models for temporal and spatial covariances.

Probably the key insight here is that the terms in a random effect (the f formula in a random-effects term (f|g) are expanded using the base-R machinery for regression model formulas. In the case of an intercept-only random effect (1|g), the model matrix is a column of ones, so we have a \(1 \times 1\) covariance matrix - a single variance. For a random-slopes model (x|g) or (1+x|g), where x is a numeric variable, the model matrix has two columns, a column of ones and column of observed values of x, and the covariance matrix is \(2 \times 2\) (intercept variance, slope variance, intercept-slope covariance).

Things start to get weird when we have (f|g) (or (1+f|g)) where f is a factor (representing a categorical variable). R uses treatment contrasts by default; if the observed values of f are c("c", "s", "v")² the corresponding factor will have a baseline level of "c" by default, and the model matrix will be:

model.matrix(~f, data.frame(f=factor(c("c", "s", "v"))))

##   (Intercept) fs fv
## 1           1  0  0
## 2           1  1  0
## 3           1  0  1
## attr(,"assign")
## [1] 0 1 1
## attr(,"contrasts")
## attr(,"contrasts")$f
## [1] "contr.treatment"

i.e., an intercept (which corresponds to the predicted mean value for observations in group c) followed by dummy variables that describe contrasts between the predicted mean values for s and c (fs) and between v and c (fv). The covariance matrix is \(3 \times 3\) and looks like this:

\[ \newcommand{\ssub}[1]{\sigma^2_{\textrm{#1}}} \newcommand{\csub}[2]{\sigma^2_{\textrm{#1}, \textrm{#2}}} \left( \begin{array}{ccc} \ssub{c} & \csub{c}{s-c} & \csub{c}{v-c} \\ \csub{c}{s-c} & \ssub{s-c} & \csub{s-c}{v-c} \\ \csub{c}{v-c} & \csub{s-c}{v-c} & \ssub{v-c} \end{array} \right) \]

This might be OK for some problems, but the parameters of the model will often be more interpretable if we remove the intercept from the formula:

model.matrix(~0+f, data.frame(f=factor(c("c", "s", "v"))))

##   fc fs fv
## 1  1  0  0
## 2  0  1  0
## 3  0  0  1
## attr(,"assign")
## [1] 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$f
## [1] "contr.treatment"

The corresponding covariance matrix is

\[ \left( \begin{array}{ccc} \ssub{c} & \csub{c}{s} & \csub{c}{v} \\ \csub{c}{s} & \ssub{s} & \csub{s}{v} \\ \csub{c}{v} & \csub{s}{v} & \ssub{v} \end{array} \right) \]

This is easier to understand (the elements are the variances of the intercepts for each group, and the covariances between intercepts of different groups). If we use an 'unstructured' model (us(f|g), or just plain (f|g)), then this reparameterization won't make any difference in the overall model fit. However, if we use a structured covariance model, then the choice matters: for example, the two models diag(f|g) and diag(0+f|g) give rise to the covariance matrices

\[ \left( \begin{array}{ccc} \ssub{c} & 0 & 0 \\ 0 & \ssub{s-c} & 0 \\ 0 & 0 & \ssub{v-c} \end{array} \right) \;\; \textrm{vs} \;\; \left( \begin{array}{ccc} \ssub{c} & 0 & 0 \\ 0 & \ssub{s} & 0 \\ 0 & 0 & \ssub{v} \end{array} \right) \]

which cannot be made equivalent by changing parameters.

What does this have to do with temporally/spatially structured covariance matrices? In this case, if two points are separated by a distance \(d_{ij}\) (in space or time), we typically want their correlation to be \(\sigma^2 \rho(d_{ij})\), where \(\rho()\) is a temporal or spatial autocorrelation function (e.g. in the AR1 model, \(\rho(d_{ij}) = \phi^{d_{ij}}\)). So we want to set up a covariance matrix

\[ \sigma^2 \left( \begin{array}{cccc} 1 & \rho(d_{12}) & \rho(d_{13}) & \ldots \\ \rho(d_{12}) & 1 & \rho(d_{23}) & \ldots \\ \rho(d_{13}) & \rho(d_{23}) & 1 & \ldots \\ \vdots & \vdots & \vdots & \ddots \end{array} \right) \]

How glmmTMB actually does this internally is to

• treat the temporal/spatial locations as a factor to construct the data structure for a \(n \times n\) covariance matrix (where \(n\) is the number of unique locations)
• use the information encoded in the levels by numFactor() to compute the corresponding pairwise distances
• use the prefix of the random term (e.g. ar1 or ou) and the autocorrelation parameters (drawn from the parameter vector) to specify the autocorrelation function
• use the autocorrelation function and the distances to fill in the values in the correlation matrix
• multiply the correlation matrix by the variance (also drawn from the parameter vector) to compute the covariance matrix

In order for this to work, we need the \(i^\textrm{th}\) column of the corresponding model matrix to correspond to an indicator variable for whether an observation is at the \(i^\textrm{th}\) location --- not to a contrast between the \(i\textrm{th}\) level and the first level! So, we want to use e.g. ar1(0 + time|g), not ar1(time|g) (which is equivalent to ar1(1+time|g)).