This vignette demonstrates some of the covariance structures available in the glmmTMB package. Currently the available covariance structures are:

Covariance Notation Parameter count Requirement
Heterogeneous unstructured us \(n(n+1)/2\)
Heterogeneous Toeplitz toep \(2n-1\)
Heterogeneous compound symmetry cs \(n+1\)
Heterogeneous diagonal diag \(n\)
AR(1) ar1 \(2\)
Ornstein–Uhlenbeck ou \(2\) Coordinates
Spatial exponential exp \(2\) Coordinates
Spatial Gaussian gau \(2\) Coordinates
Spatial Matern mat \(3\) Coordinates

The word ‘heterogeneous’ refers to the marginal variances of the model. Beyond correlation parameters, a heteorogenous structure uses \(n\) additional variance parameters where \(n\) is the dimension.

Some of the structures require temporal or spatial coordinates. We will show examples of this in a later section.

The AR(1) covariance structure

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

Description Parameter Value
Mean \(\mu\) 0
Process variance \(\sigma^2\) 1
Measurement variance \(\sigma_0^2\) 1
One-step correlation \(e^{-\theta}\) 0.7

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

n <- 6                                              ## Number of time points
x <- 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.

times <- factor(1:n)
levels(times)
## [1] "1" "2" "3" "4" "5" "6"

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

Now fit the model using

glmmTMB(y ~ ar1(times + 0 | group))

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 \(e^{-\theta}\) (First off-diagonal of “Corr”) in the output:

FIXME: Try a longer time series when the print.VarCorr is fixed.

## Formula:          y ~ ar1(times + 0 | group)
##      AIC      BIC   logLik df.resid 
## 17.07962 16.24666 -4.53981        2 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups   Name   Std.Dev.  Corr                         
##  group    times1 0.5753756                              
##           times2 0.5753756 -0.48                        
##           times3 0.5753756  0.23 -0.48                  
##           times4 0.5753756 -0.11  0.23 -0.48            
##           times5 0.5753756  0.05 -0.11  0.23 -0.48      
##           times6 0.5753756 -0.03  0.05 -0.11  0.23 -0.48
##  Residual        0.0000137                              
## 
## Number of obs: 6 / Conditional model: group, 1
## 
## Dispersion estimate for gaussian family (sigma^2): 1.88e-10 
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)  
##       1.064

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) {
    x <- mvrnorm(mu = rep(0,n),
             Sigma = .7 ^ 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)
##             y times group
## 1  0.43831650     1     1
## 2 -0.15292422     2     1
## 3  0.82537310     3     1
## 4  1.13161026     4     1
## 5  1.91410931     5     1
## 6  0.06970856     6     1

A dataset with 1000 groups is generated:

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

And fitting the model on this larger dataset gives estimates close to the true values:

fit.ar1 <- glmmTMB(y ~ ar1(times + 0 | group), data=dat)
fit.ar1
## Formula:          y ~ ar1(times + 0 | group)
## Data: dat
##       AIC       BIC    logLik  df.resid 
##  20554.47  20581.27 -10273.24      5996 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups   Name   Std.Dev. Corr                    
##  group    times1 1.002                            
##           times2 1.002    0.72                    
##           times3 1.002    0.52 0.72               
##           times4 1.002    0.38 0.52 0.72          
##           times5 1.002    0.27 0.38 0.52 0.72     
##           times6 1.002    0.20 0.27 0.38 0.52 0.72
##  Residual        1.020                            
## 
## Number of obs: 6000 / Conditional model: group, 1000
## 
## Dispersion estimate for gaussian family (sigma^2): 1.04 
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)  
##     0.04378

The unstructured covariance

We can try to fit an unstructured covariance to the previous dataset dat. For this case an unstructered covariance has 15 correlation parameters and 6 variance parameters. Adding \(\sigma_0^2 I\) on top would cause a strict overparameterization. Hence, when fitting the model with glmmTMB, we have to disable the \(\varepsilon\) term (the dispersion):

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

The estimated variance and correlation parameters are:

VarCorr(fit.us)
## 
## Conditional model:
##  Groups   Name   Std.Dev.   Corr                         
##  group    times1 1.41052242                              
##           times2 1.39875744 0.354                        
##           times3 1.48794614 0.267 0.384                  
##           times4 1.42176303 0.166 0.230 0.349            
##           times5 1.42250955 0.154 0.193 0.227 0.353      
##           times6 1.44075053 0.107 0.143 0.192 0.280 0.344
##  Residual        0.00012207

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 disperison is effectively included in the estimated covariance matrix.

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 5 correlation parameters and 6 variance parameters. This time we do not have to disable the dispersion parameter.

fit.toep <- glmmTMB(y ~ toep(times + 0 | group), data=dat)
fit.toep$sdr$pdHess ## Converged ?
## [1] TRUE

The estimated variance and correlation parameters are:

VarCorr(fit.toep)
## 
## Conditional model:
##  Groups   Name   Std.Dev. Corr                         
##  group    times1 1.10964                               
##           times2 1.09543  0.567                        
##           times3 1.20845  0.396 0.567                  
##           times4 1.12680  0.299 0.396 0.567            
##           times5 1.12423  0.236 0.299 0.396 0.567      
##           times6 1.14961  0.168 0.236 0.299 0.396 0.567
##  Residual        0.86955

The residual variance appears downward biased. REML estimation (currently not part of glmmTMB) would probably give a better estimate of the variance and thereby the correlation parameters.

FIXME: Add REML argument to glmmTMB

Compound symmetry

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

fit.cs <- glmmTMB(y ~ cs(times + 0 | group), data=dat)
## Warning in glmmTMB(y ~ cs(times + 0 | group), data = dat): Model
## convergence problem. Hessian is not positive definite. This may indicate
## that a model is overparameterized.
fit.cs$sdr$pdHess ## Converged ?
## [1] FALSE

The estimated variance and correlation parameters are:

VarCorr(fit.cs)
## 
## Conditional model:
##  Groups   Name   Std.Dev. Corr                         
##  group    times1 0.60040                               
##           times2 0.71307  1.000                        
##           times3 0.85664  1.000 1.000                  
##           times4 0.77512  1.000 1.000 1.000            
##           times5 0.71430  1.000 1.000 1.000 1.000      
##           times6 0.63714  1.000 1.000 1.000 1.000 1.000
##  Residual        1.23565

Anova tables

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

anova(fit.ar1, fit.toep, fit.us)
## Data: dat
## Models:
## fit.ar1: y ~ ar1(times + 0 | group)
## fit.toep: y ~ toep(times + 0 | group)
## fit.us: y ~ us(times + 0 | group)
##          Df   AIC   BIC logLik deviance  Chisq Chi Df Pr(>Chisq)
## fit.ar1   4 20554 20581 -10273    20546                         
## fit.toep 13 20566 20653 -10270    20540 6.6986      9     0.6685
## fit.us   22 20578 20725 -10267    20534 5.9611      9     0.7438

The models cs is a sub-model of toep:

anova(fit.cs, fit.toep)
## Data: dat
## Models:
## fit.cs: y ~ cs(times + 0 | group)
## fit.toep: y ~ toep(times + 0 | group)
##          Df   AIC   BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## fit.cs    9                                                    
## fit.toep 13 20566 20653 -10270    20540            4

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)
pos
##  [1] (2,1) (2,2) (2,1) (1,2) (2,1) (1,1) (1,2) (2,2) (1,1) (1,2)
## Levels: (1,1) (2,1) (1,2) (2,2)

Numeric coordinates can be recovered from the factor levels

parseNumLevels(levels(pos))
##      [,1] [,2]
## [1,]    1    1
## [2,]    2    1
## [3,]    1    2
## [4,]    2    2

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

dat$times <- numFactor(dat$times)
levels(dat$times)
## [1] "(1)" "(2)" "(3)" "(4)" "(5)" "(6)"

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=dat)
fit.ou$sdr$pdHess ## Converged ?
## [1] TRUE

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

VarCorr(fit.ou)
## 
## Conditional model:
##  Groups   Name     Std.Dev. Corr                         
##  group    times(1) 1.0021                                
##           times(2) 1.0021   0.722                        
##           times(3) 1.0021   0.521 0.722                  
##           times(4) 1.0021   0.376 0.521 0.722            
##           times(5) 1.0021   0.272 0.376 0.521 0.722      
##           times(6) 1.0021   0.196 0.272 0.376 0.521 0.722
##  Residual          1.0205

However, note the differences between ou and ar1:

Spatial correlations

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

FIXME: Maybe try some spatial data instead ?

Matern

fit.mat <- glmmTMB(y ~ mat(times + 0 | group), data=dat, dispformula=~0)
fit.mat$sdr$pdHess ## Converged ?
## [1] TRUE
VarCorr(fit.mat)
## 
## Conditional model:
##  Groups   Name     Std.Dev.   Corr                         
##  group    times(1) 1.43028671                              
##           times(2) 1.43028671 0.356                        
##           times(3) 1.43028671 0.250 0.356                  
##           times(4) 1.43028671 0.187 0.250 0.356            
##           times(5) 1.43028671 0.143 0.187 0.250 0.356      
##           times(6) 1.43028671 0.112 0.143 0.187 0.250 0.356
##  Residual          0.00012207

Gaussian

fit.gau <- glmmTMB(y ~ gau(times + 0 | group), data=dat, dispformula=~0)
fit.gau$sdr$pdHess ## Converged ?
## [1] TRUE
VarCorr(fit.gau)
## 
## Conditional model:
##  Groups   Name     Std.Dev.   Corr                         
##  group    times(1) 1.41782818                              
##           times(2) 1.41782818 0.279                        
##           times(3) 1.41782818 0.006 0.279                  
##           times(4) 1.41782818 0.000 0.006 0.279            
##           times(5) 1.41782818 0.000 0.000 0.006 0.279      
##           times(6) 1.41782818 0.000 0.000 0.000 0.006 0.279
##  Residual          0.00012207

Exponential

fit.exp <- glmmTMB(y ~ exp(times + 0 | group), data=dat)
fit.exp$sdr$pdHess ## Converged ?
## [1] TRUE
VarCorr(fit.exp)
## 
## Conditional model:
##  Groups   Name     Std.Dev. Corr                         
##  group    times(1) 1.0021                                
##           times(2) 1.0021   0.722                        
##           times(3) 1.0021   0.521 0.722                  
##           times(4) 1.0021   0.376 0.521 0.722            
##           times(5) 1.0021   0.272 0.376 0.521 0.722      
##           times(6) 1.0021   0.196 0.272 0.376 0.521 0.722
##  Residual          1.0205