Efficient algorithms for CAR models have been widely discussed in particular in the econometric literature (e.g., LeSage and Pace 2009), but these models are not necessarily recommended for irregular lattices (see Wall, 2004 and Martellosio, 2012 for some insights on the implications of autoregressive models).
In CAR models, the covariance matrix of random effects u can be described as \(\lambda\)(I\(-\rho\) W\()^{-1}\) where W is the (symmetric) adjacency matrix. HLCor uses the spectral decomposition of the adjacency matrix, written as W=VDV' where D is a diagonal matrix of eigenvalues \(d_i\). The covariance of V'u is
\(\lambda\)(I\(-\rho\) D\()^{-1}\), which is a diagonal matrix with elements
\(\lambda_i\)=\(\lambda\)/(1\(-\rho d_i\)). Hence \(1/\lambda_i\) is in the linear predictor form \(\alpha\)+\(\beta d_i\) This can be used to fit \(\lambda\) and \(\rho\) efficiently. A call to corrHLfit with the additional argument
init.HLfit=list(rho=0) should be equivalent in speed and result to the HLCor call.
This is fast for small datasets (as in the example below) but more generic maximization algorithms may be preferable for large ones. It is suggested to use fitme generally unless one has a large number of small data sets to analyze. A call to fitme or corrHLfit without that initial value does not use the spectral decomposition. It performs numerical maximization of the likelihood (or restricted likelihood) as function of the correlation parameter \(\rho\). The choice of fitting function may slightly impact the results. The ML fits by corrHLfit and HLCor should be practically equivalent. The REML fits should slightly differ from each other, due to the fact that the REML approximation for GLMMs does not maximize a single likelihood function.
If HLCor is used, the results are reported as the coefficients \(\alpha\) ((Intercept)) and \(\beta\) (adjd) of the predictor for \(1/\lambda_i\), in addition to the resulting values of \(\rho\) and of the common \(\lambda\) factor.
Different fits may also differ in using or not algorithms that exploit the sparsity of the precision matrix of the autoregressive random effect. By default, spaMM tends to select sparse-precision algorithms for large datasets and large (i.e. many-level) random effects (details are complex).
However, for AR1 models, the dimension of the implied precision matrix is determined by the extreme values of grouping factor (typically interpreted as a time index), as all intermediate values must be considered. Then, the correlation-based algorithms may be more efficient if only a few levels are present in the data, as only a small correlation matrix is required in that case.