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
fit.mat <- glmmTMB(y ~ mat(times + 0 | group), data=dat1, dispformula=~0)
fit.mat$sdr$pdHess ## Converged ?## [1] TRUEVarCorr(fit.mat)##
## Conditional model:
## Groups Name Std.Dev. Corr
## group times(1) 1.42996896
## times(2) 1.42996896 0.357
## times(3) 1.42996896 0.250 0.357
## times(4) 1.42996896 0.187 0.250 0.357
## times(5) 1.42996896 0.143 0.187 0.250 0.357
## times(6) 1.42996896 0.111 0.143 0.187 0.250 0.357
## Residual 0.00012207Gaussian
“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 ?## [1] TRUEVarCorr(fit.gau)##
## Conditional model:
## Groups Name Std.Dev. Corr
## group times(1) 1.41749843
## times(2) 1.41749843 0.280
## times(3) 1.41749843 0.006 0.280
## times(4) 1.41749843 0.000 0.006 0.280
## times(5) 1.41749843 0.000 0.000 0.006 0.280
## times(6) 1.41749843 0.000 0.000 0.000 0.006 0.280
## Residual 0.00012207Exponential
fit.exp <- glmmTMB(y ~ exp(times + 0 | group), data=dat1)
fit.exp$sdr$pdHess ## Converged ?## [1] TRUEVarCorr(fit.exp)##
## Conditional model:
## Groups Name Std.Dev. Corr
## group times(1) 1.0034
## times(2) 1.0034 0.721
## times(3) 1.0034 0.520 0.721
## times(4) 1.0034 0.375 0.520 0.721
## times(5) 1.0034 0.270 0.375 0.520 0.721
## times(6) 1.0034 0.195 0.270 0.375 0.520 0.721
## Residual 1.0187A 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")
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")## 2.5 % 97.5 % Estimate
## sigma 9.492257 16.46098 12.50007The 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## pos group
## 1 (1,1) 1
## 2 (2,1) 1
## 3 (3,1) 1
## 4 (1,2) 1
## 5 (2,2) 1
## 6 (3,2) 1
## 7 (1,3) 1
## 8 (2,3) 1
## 9 (3,3) 1and predict using
predict(f, newdata, type="response", allow.new.levels=TRUE)## [1] 104.4634 104.8042 105.1246 104.6228 104.9824 105.3195 104.7828 105.1666
## [9] 105.5366A 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")
