library(knitr) library(glmmTMB) library(MASS) ## for mvrnorm() library(TMB) ## for tmbprofile() library(mvabund) ## for spider data ## devtools::install_github("kaskr/adcomp/TMB") ## get development version knitr::opts_chunk$set(echo = TRUE, eval=if (exists("params")) params$EVAL else FALSE) do_image <- exists("params") && params$EVAL ## want to *store* images within package save_vig_dir <- file.path("inst","vignette_data") pkg_dir <- "glmmTMB" ## guess where we are ... if (grepl("/vignettes$",getwd())) { ## in vignettes dir save_vig_dir <- file.path("../",save_vig_dir) } else if (grepl(paste0("/",pkg_dir,"$"),getwd())) { ## in repo head save_vig_dir <- file.path(pkg_dir,save_vig_dir) } ## want to *retrieve* images from system files use_vig_dir <- system.file("vignette_data",package="glmmTMB") mkfig <- function(expr,fn) { png(normalizePath(file.path(save_vig_dir,fn))) eval(substitute(expr)) invisible(dev.off()) } usefig <- function(fn) { knitr::include_graphics(file.path(use_vig_dir,fn)) } ## turned off caching for now: got error in chunk 'fit.us.2' ## Error in retape() : ## Error when reading the variable: 'thetaf'. Please check data and parameters. ## In addition: Warning message: ## In retape() : Expected object. Got NULL. set.seed(1) ## run this in interactive session if you actually want to evaluate chunks ... ## Sys.setenv(NOT_CRAN="true")
This vignette demonstrates some of the covariance structures available in the glmmTMB
package.
Currently the available covariance structures are:
ctab <- read.delim(sep = "#", comment = "", header = TRUE, check.names = FALSE, text = " Covariance # Notation # no. parameters # Requirement # Parameters Unstructured (general positive definite) # `us` # $n(n+1)/2$ # # See [Mappings] Heterogeneous Toeplitz # `toep` # $2n-1$ # # log-SDs ($\\theta_1-\\theta_n$); correlations $\\rho_k = \\theta_{n+k}/\\sqrt{1+\\theta_{n+k}^2}$, $k = \\textrm{abs}(i-j+1)$ Het. compound symmetry # `cs` # $n+1$ # # log-SDs ($\\theta_1-\\theta_n$); correlation $\\rho = \\theta_{n+1}/\\sqrt{1+\\theta_{n+1}^2}$ Het. diagonal # `diag` # $n$ # # log-SDs AR(1) # `ar1` # $2$ # Unit spaced levels # log-SD; $\\rho = \\left(\\theta_2/\\sqrt{1+\\theta_2^2}\\right)^{d_{ij}}$ Ornstein-Uhlenbeck # `ou` # $2$ # Coordinates # log-SD; log-OU rate ($\\rho = \\exp(-\\exp(\\theta_2) d_{ij})$) Spatial exponential # `exp` # $2$ # Coordinates # log-SD; log-scale ($\\rho = \\exp(-\\exp(-\\theta_2) d_{ij})$) Spatial Gaussian # `gau` # $2$ # Coordinates # log-SD; log-scale ($\\rho = \\exp(-\\exp(-2\\theta_2) d_{ij}^2$) Spatial Matèrn # `mat` # $3$ # Coordinates # log-SD, log-range, log-shape (power) Reduced rank # `rr` # $nd-d(d-1)/2$ # rank (d) " ) knitr::kable(ctab)
The word 'heterogeneous' refers to the marginal variances of the model.
Some of the structures require temporal or spatial coordinates. We will show examples of this in a later section.
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
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 | $\phi$ | 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 <- 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.
times + 0
corresponds to a design
matrix $Z$ linking observation vector $y$ (rows) with a random
effects vector $u$ (columns) (see [Construction of structured covariance matrices] for why we need the + 0
)ar1
(this is the only glmmTMB
specific part of the formula).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:
glmmTMB(y ~ ar1(times + 0 | group), data=dat0)
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.)
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, 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)
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))
We can try to fit an unstructured covariance to the previous dataset
dat
. For this case an unstructured covariance has r (n*n-n)/2
correlation parameters and r n
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)
\newcommand{\textsub}[2]{#1_{{\text {#2}}}} 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 \textsub{\sigma^2}{AR}/(\textsub{\sigma^2}{AR} + \textsub{\sigma^2}{meas})$).
The next natural step would be to reduce the number of parameters by
collecting correlation parameters within the same off-diagonal. This
amounts to r (n-1)
correlation parameters and r n
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()
)^[Why do we do this? Consider the slightly simplified case of a homogeneous Toeplitz structure where all of the variance parameters are identical. The diagonal elements of the covariance matrix are equal to $\sigma_t^2$, the off-diagonals to $\sigma_t^2 \cdot \rho(|i-j|)$. If we add a residual variance to the model then the diagonal of the combined covariance matrix becomes $\sigma_t^2 + \sigma_r^2$ and the off-diagonals become $(\sigma_t^2 + \sigma_r^2) \rho(|i-j|)$. However, by reparameterizing the Toeplitz model to ${{\sigma_t^2}' = \sigma_t^2 + \sigma_r^2, \rho'(|i-j|) = \rho(|i-j|) \cdot \frac{\sigma_t^2}{\sigma_t^2 + \sigma_r^2}}$ --- that is, by inflating the variance and deflating the correlation parameters --- we can get back to an equivalent Toeplitz model. This implies that the residual variance and the Toeplitz covariance parameters are jointly unidentifiable, which is likely to make problems for the fitting procedure.]
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 ($\textsub{\sigma^2}{AR} + \textsub{\sigma^2}{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)])
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)
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).
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)
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.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)
"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)
fit.exp <- glmmTMB(y ~ exp(times + 0 | group), data=dat1) fit.exp$sdr$pdHess ## Converged ?
VarCorr(fit.exp)
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)
if (do_image) { volcano.data <- array(NA, dim(volcano)) volcano.data[cbind(d$x, d$y)] <- d$z mkfig(image(volcano.data, main="Spatial data"),"volcano_data.png") }
usefig("volcano_data.png")
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")
if (do_image) { mkfig(image(pred, main="Reconstruction", useRaster=TRUE), "volcano_results.png") }
usefig("volcano_results.png")
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.
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)
mkfig(plot(tt),"us_profile_plot.png")
usefig("us_profile_plot.png")
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)
mkfig(plot(tt2),"cs_profile_plot.png")
usefig("cs_profile_plot.png")
Consider a generalized linear mixed model
\begin{equation} g(\boldsymbol{\mu}) = \boldsymbol{X\beta} + \boldsymbol{Zb} \end{equation}
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})$.
A general latent variable model (GLVM) requires many fewer parameters for the variance-covariance matrix, $\boldsymbol{\Sigma}$. To a fit a GLVM we add a reduced-rank (rr) covariance structure, so the model becomes \begin{align} g(\boldsymbol{\mu}) &= \boldsymbol{X\beta} + \boldsymbol{Z(I_n \otimes \Lambda)b} \ &= \boldsymbol{X\beta} + \boldsymbol{Zb_{new}} \end{align} 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. @niku2019gllvm 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.
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) [@hui2015model].
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")
^[chocolate, strawberry, vanilla] 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"))))
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"))))
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
numFactor()
to compute the corresponding pairwise distancesar1
or ou
) and the autocorrelation parameters (drawn from the parameter vector) to specify the autocorrelation functionIn 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)
).
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.