Covariance structures with glmmTMB

Demonstration on simulated data

First, let's consider a simple time series model. Assume that our measurements \(Y(t)\) are given at discrete times \(t \in \{1,...,n\}\) by

\[Y(t) = \mu + X(t) + \varepsilon(t)\]

where

• \(\mu\) is the mean value parameter.
• \(X(t)\) is a stationary AR(1) process, i.e. has covariance \(cov(X(s), X(t)) = \sigma^2\exp(-\theta |t-s|)\).
• \(\varepsilon(t)\) is iid. \(N(0,\sigma_0^2)\) measurement error.

A simulation experiment is set up using the parameters

The following R-code draws a simulation based on these parameter values. For illustration purposes we consider a very short time series.

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

In order to fit the model with glmmTMB we must first specify a time variable as a factor. The factor levels correspond to unit spaced time points. It is a common mistake to forget some factor levels due to missing data or to order the levels incorrectly. We therefore recommend to construct factors with explicit levels, using the levels argument to the factor function:

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

We also need a grouping variable. In the current case there is only one time-series so the grouping is:

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

We combine the data into a single data frame (not absolutely required, but good practice):

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

Now fit the model using

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

This formula notation follows that of the lme4 package.

• The left hand side of the bar times + 0 corresponds to a design matrix \(Z\) linking observation vector \(y\) (rows) with a random effects vector \(u\) (columns).
• The distribution of \(u\) is ar1 (this is the only glmmTMB specific part of the formula).
• The right hand side of the bar splits the above specification independently among groups. Each group has its own separate \(u\) vector but shares the same parameters for the covariance structure.

After running the model, we find the parameter estimates \(\mu\) (intercept), \(\sigma_0^2\) (dispersion), \(\sigma\) (Std. Dev.) and \(\phi\) (First off-diagonal of "Corr") in the output:

For those trying to make sense of the internal parameterization, the internal transformation from the parameter (\(\phi\)) to the AR1 coefficient is \(\phi = theta_2/sqrt(1+\theta_2^2)\); the inverse transformation is \(\theta_2 = \phi/sqrt(1-\phi^2)\). (The first element of the theta vector is the log-standard-deviation.)

Increasing the sample size

A single time series of 6 time points is not sufficient to identify the parameters. We could either increase the length of the time series or increase the number of groups. We'll try the latter:

simGroup <- function(g, n=6, rho=0.7) {
    x <- MASS::mvrnorm(mu = rep(0,n),
             Sigma = rho ^ 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)

Generate a dataset with 1000 groups:

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

And fitting the model on this larger dataset gives estimates close to the true values (AR standard deviation=1, residual (measurement) standard deviation=1, autocorrelation=0.7):

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

Matern

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

VarCorr(fit.mat)

Gaussian

"Gaussian" refers here to a Gaussian decay in correlation with distance, i.e. \(\rho = \exp(-d x^2)\), not to the conditional distribution ("family").

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

VarCorr(fit.gau)

Exponential

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

VarCorr(fit.exp)

A spatial covariance example

Starting out with the built in volcano dataset we reshape it to a data.frame with pixel intensity z and pixel position x and y:

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

Next, add random normal noise to the pixel intensities and extract a small subset of 100 pixels. This is our spatial dataset:

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

Display sampled noisy volcano data:

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

Based on this data, we'll attempt to re-construct the original image.

As model, it is assumed that the original image image(volcano) is a realization of a random field with correlation decaying exponentially with distance between pixels.

Denoting by \(u(x,y)\) this random field the model for the observations is

\[ z_{i} = \mu + u(x_i,y_i) + \varepsilon_i \]

To fit the model, a numFactor and a dummy grouping variable must be added to the dataset:

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

The model is fit by

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

Recall that a standard deviation sd=15 was used to distort the image. A confidence interval for this parameter is

confint(f, "sigma")

The glmmTMB predict method can predict unseen levels of the random effects. For instance to predict a 3-by-3 corner of the image one could construct the new data:

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

and predict using

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

A specific image column can thus be predicted using the function

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)
}

Prediction of the entire image is carried out by (this takes a while...):

pred <- sapply(1:61, predict_col)

Finally plot the re-constructed image by

image(pred, main="Reconstruction")

Unstructured

For an unstructured matrix of size n, parameters 1:n represent the log-standard deviations while the remaining n(n-1)/2 (i.e. (n+1):(n:(n*(n+1)/2))) are the elements of the scaled Cholesky factor of the correlation matrix, filled in row-wise order (see TMB documentation). In particular, if \(L\) is the lower-triangular matrix with 1 on the diagonal and the correlation parameters in the lower triangle, then the correlation matrix is defined as \(\Sigma = D^{-1/2} L L^\top D^{-1/2}\), where \(D = \textrm{diag}(L L^\top)\). For a single correlation parameter \(\theta_0\), this works out to \(\rho = \theta_0/\sqrt{1+\theta_0^2}\).

(See calculations here.)

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]]))

Profiling (experimental/exploratory):

## 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)

Reduced rank

One of our main motivations for adding this variance-covariance structure is to enable the analysis of multivariate abundance data, for example to model the abundance of different taxa across multiple sites. Typically an unstructured random effect is assumed to account for correlation between taxa; however the number of parameters required quickly becomes large with increasing numbers of taxa. A GLVM is a flexible and more parsimonious way to account for correlation so that one can fit a joint model across many taxa.

A GLVM can be fit by specifying a reduced rank (rr) covariance structure. For example, the code for modeling the mean abundance against taxa and to account for the correlation between taxa using two latent variables is as follows

if (require(mvabund)) {
    data(spider)
    ## organize data into long format
    sppTot <- sort(colSums(spider$abund), decreasing = TRUE)
    tmp <- cbind(spider$abund, spider$x)
    tmp$id <- 1:nrow(tmp)
    spiderDat <- reshape(tmp,
                         idvar = "id",
                         timevar = "Species",
                         times =  colnames(spider$abund),
                         varying = list(colnames(spider$abund)),
                         v.names = "abund",
                         direction = "long")
    ## fit rank-reduced models with varying dimension
    fit_list <- lapply(2:10,
                       function(d) {
                           fit.rr <- glmmTMB(abund ~ Species + rr(Species + 0|id, d = d),
                                             data = spiderDat)
                       })
    ## compare fits via AIC
    aic_vec <- sapply(fit_list, AIC)
    aic_vec - min(aic_vec, na.rm = TRUE)

The left hand side of the bar taxa + 0 corresponds to a factor loading matrix that accounts for the correlations among taxa. The right hand side of the bar splits the above specification independently among sites. The d is a non-negative integer (which defaults to 2).

An option in glmmTMBControl() has been included to initialize the starting values for the parameters based on the approach mentioned above with the default set at glmmTMBControl(start_method = list(method = NULL, jitter.sd = 0):

• method = "res" initializes starting values from the results of fitting a GLM, and fitting a reduced rank model to the residuals to obtain starting values for the fixed coefficients, the latent variables and the factor loadings.
• jitter.sd adds variation to the starting values of latent variables when method = "res" (default 0).

For a reduced rank matrix of rank d, parameters 1:d represent the diagonal factor loadings while the remaining \(nd-d(d-3)/2\), (i.e. parameters (d+1):(nd-d(d-1)/2) are the lower diagonal factor loadings filled in column-wise order. The factor loadings from a model can be obtained by fit.rr$obj$env$report(fit.rr$fit$parfull)$fact_load[[1]]. An appropriate rank for the model can be determined by standard model selection approaches such as information criteria (e.g. AIC or BIC) (F. K. Hui et al. 2015).