Covariance structures with glmmTMB
This vignette demonstrates some of
the covariance structures available in the glmmTMB package.
Currently the available covariance structures are:
| Covariance | Notation | no. parameters | Requirement | Parameters |
|---|---|---|---|---|
| Unstructured (general positive definite) | us |
n(n + 1)/2 | See Mappings | |
| Heterogeneous Toeplitz | toep |
2n − 1 | log-SDs (θ1 − θn); correlations $\rho_k = \theta_{n+k}/\sqrt{1+\theta_{n+k}^2}$, k = abs(i − j + 1) | |
| Het. compound symmetry | cs |
n + 1 | log-SDs (θ1 − θn); correlation $\rho = \theta_{n+1}/\sqrt{1+\theta_{n+1}^2}$ | |
| Homogenous diagonal | homdiag |
1 | log-SD | |
| 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 (ρ = exp (−exp (θ2)dij)) |
| Spatial exponential | exp |
2 | Coordinates | log-SD; log-scale (ρ = exp (−exp (−θ2)dij)) |
| Spatial Gaussian | gau |
2 | Coordinates | log-SD; log-scale (ρ = exp (−exp (−2θ2)dij2) |
| Spatial Matèrn | mat |
3 | Coordinates | log-SD, log-range, log-shape (power) |
| Reduced-rank | rr |
nd − d(d − 1)/2 | rank (d) | |
| Proptional | propto |
1 | Variance-covariance matrix |
The word ‘heterogeneous’ refers to the marginal variances of the model.
Homogenous versions of some structures (e.g. Toeplitz, compound
symmetric) can be implemented by using the map argument to
set all log-SD parameters equal to each other.
Some of the structures require temporal or spatial coordinates. We will show examples 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 ∈ {1, ..., n} by
Y(t) = μ + X(t) + ε(t)
where
- μ is the mean value parameter.
- X(t) is a stationary AR(1) process, i.e. has covariance cov(X(s), X(t)) = σ2exp (−θ|t − s|).
- ε(t) is iid. N(0, σ02) measurement error.
A simulation experiment is set up using the parameters
| Description | Parameter | Value |
|---|---|---|
| Mean | μ | 0 |
| Process variance | σ2 | 1 |
| Measurement variance | σ02 | 1 |
| One-step correlation | ϕ | 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 noiseIn 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:
We also need a grouping variable. In the current case there is only one time-series so the grouping is:
We combine the data into a single data frame (not absolutely required, but good practice):
Now fit the model using
This formula notation follows that of the lme4
package.
- The left hand side of the bar
times + 0corresponds 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) - The distribution of u is
ar1(this is the onlyglmmTMBspecific 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 μ (intercept), σ02 (dispersion), σ (Std. Dev.) and ϕ (First off-diagonal of “Corr”) in the output:
For those trying to make sense of the internal parameterization, the
internal transformation from the parameter (θ2) 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, 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:
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):
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 σ02I
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
ε 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:
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. ρ′ = ρσAR2/(σAR2 + σmeas2)).
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())1
fit.toep <- glmmTMB(y ~ toep(times + 0 | group), data=dat1,
dispformula=~0)
fit.toep$sdr$pdHess ## Converged ?The estimated variance and correlation parameters are:
The diagonal elements are all approximately equal to the true total variance (σAR2 + σmeas2=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):
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:
Anova tables
The models ar1, toep, and us
are nested so we can use:
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:
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
we can generate a factor representing (x, y) coordinates by
Numeric coordinates can be recovered from the factor levels:
In order to try the remaining structures on our test data we
re-interpret the time factor using numFactor:
Ornstein–Uhlenbeck
Having the numeric times encoded in the factor levels we can now try the Ornstein–Uhlenbeck covariance structure.
It should give the exact same results as ar1 in this
case since the times are equidistant:
However, note the differences between ou and
ar1:
oucan handle irregular time points.ouonly 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.
Matern
Gaussian
“Gaussian” refers here to a Gaussian decay in correlation with distance, i.e. ρ = exp (−dx2), not to the conditional distribution (“family”).
Exponential
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:
Next, add random normal noise to the pixel intensities and extract a small subset of 100 pixels. This is our spatial dataset:
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
zi = μ + u(xi, yi) + εi
To fit the model, a numFactor and a dummy grouping
variable must be added to the dataset:
The model is fit by
Recall that a standard deviation sd=15 was used to
distort the image. A confidence interval for this parameter is
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)))
newdataand predict using
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…):
Finally plot the re-constructed image by
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.
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 Σ = D−1/2LL⊤D−1/2,
where D = diag(LL⊤).
For a single correlation parameter θ0, this works out to
$\rho = \theta_0/\sqrt{1+\theta_0^2}$
(with inverse $\theta_0
= \rho/\sqrt(1-\rho^2)$. You can use the utility functions
get_cor() (transform a theta vector into the
upper triangular [rowwise] elements of a correlation matrix, or the full
correlation matrix) and put_cor() (translate a correlation
matrix, or the values from the lower triangle, into a theta
vector) to perform these transformations.
(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)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)
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)
Generalized latent variable model
Consider a generalized linear mixed model
where g(.) is the link function; β is a p-dimensional vector of regression coefficients related to the covariates; X is an n × p model matrix; and Z is the n × q model matrix for the q-dimensional vector-valued random effects variable U which is multivariate normal with mean zero and a parameterized q × q variance-covariance matrix, Σ, i.e., U ∼ N(0, Σ).
A general latent variable model (GLVM) requires many fewer parameters for the variance-covariance matrix, Σ. To a fit a GLVM we add a reduced-rank (rr) covariance structure, so the model becomes where ⊗ is the Kronecker product and Λ = (λ1, …, λd)′ is the q × d matrix of factor loadings (with d ≪ q). The upper triangular elements of Λ are set to be zero to ensure parameter identifiability. Here we assume that the latent variables follow a multivariate standard normal distribution, b ∼ N(0, 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.
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
## 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)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.sdadds variation to the starting values of latent variables whenmethod = "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)
(Hui et al. 2015).
We can extract the random effects (predicted values for each site by
species combination) with ranef();
`as.data.frame(ranef()) (or
broom.mixed::tidy(..., effects = "ran_vals")) gives the
results in a more convenient format. Based on this information, we can
plot the predictions for species (ordered by their predicted presence at
site 1). (We’ve arbitrarily chosen d=3 here.)
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")
}If we instead want to get the factor loadings by Species and
latent variables by site, we can use a (so far experimental)
function to get a list with components $fl (factor
loadings) and $b (latent variables by site)
source(system.file("misc", "extract_rr.R", package = "glmmTMB"))
rr_info <- extract_rr(spider_rr)
lapply(rr_info, dim)We can use this information to create an (ugly) biplot. (Improvements welcome!)
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)Proportional
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 × 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 × 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")2 the corresponding factor will have a
baseline level of "c" by default, and the model matrix will
be:
## (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 × 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:
## 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 dij (in space or time), we typically want their correlation to be σ2ρ(dij), where ρ() is a temporal or spatial autocorrelation function (e.g. in the AR1 model, ρ(dij) = ϕdij). 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 × 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.
ar1orou) 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 ith column of the
corresponding model matrix to correspond to an indicator variable for
whether an observation is at the ith location —
not to a contrast between the ith 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)).
References
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 σt2, the off-diagonals to σt2 ⋅ ρ(|i − j|). If we add a residual variance to the model then the diagonal of the combined covariance matrix becomes σt2 + σr2 and the off-diagonals become (σt2 + σr2)ρ(|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.↩︎
chocolate, strawberry, vanilla↩︎