Example of making a spatial kernel and using it as the inverse G matrix in a mixed-effects model with MCMCglmm.

geoK

# Spatial kernel for mixed-effects model

library(geosphere) # for distances between countries
library(ggplot2)

# Life-expectancy (le) and center coordinates
geo <- structure(list(ISO = structure(1:10, .Label = c("ALB", "BIH", "GRC", "HRV", "ITA", "MKD", "MLT", "MNE", "SRB", "TUN"), class = "factor"), Country = structure(c(1L, 2L, 4L, 3L, 5L, 6L, 7L, 8L, 9L, 10L), .Label = c("Albania", "Bosnia and Herzegovina", "Croatia", "Greece", "Italy", "Macedonia, FYR", "Malta", "Montenegro", "Serbia", "Tunisia"), class = "factor"), le = c(70.4, 70.7, 75.2, 70.5, 74.5, 69.6, 73.6, 74.1, 70.2, 64.1), longitude = c(20, 18, 22, 15.5, 12.8333, 22, 14.5833, 19, 21, 9), latitude = c(41, 44, 39, 45.1667, 42.8333, 41.8333, 35.8333, 42, 44, 34)), .Names = c("ISO", "Country", "le", "longitude", "latitude"), row.names = 1:10, class = "data.frame")


latlong <- geo[c('longitude', 'latitude')]

# Calculate great circle distances
latlong_dist <- distm(latlong)

# Tuning parameter (this was determined by cross-validation)
h <- 1 / 1e12
# distance kernel. 1 = close, 0 = far away
K <- exp(-h * latlong_dist^2)
# name rows and columns for data points
# they represent
rownames(K) <- colnames(K) <- geo$ISO

# Plot example distances for Malta
malta <- data.frame(long=geo$longitude, lat=geo$latitude, K=K['MLT',], group='target')

world <- map_data('world')
ggplot(world, aes(x=long, y=lat, group=group)) + geom_path() + geom_point(aes(col=K), data=malta, group=NULL, size=3)+ scale_y_continuous(breaks=(-2:2) * 30) + scale_x_continuous(breaks=(-4:4) * 45) + coord_map('ortho', orientation=c(35, 15, 0), xlim=c(0, 30), ylim=c(20, 60))

# Invert the matrix
Kinv <- chol2inv(chol(K))
# if matrix does not invert because it is not
# positive definite (check for negative values
# in eigen(K)$values) you can smooth it by
# adding small values to the diagonals e.g.,
#   chol2inv(chol(K + diag(.01, nrow(K))))

rownames(Kinv) <- colnames(Kinv) <- rownames(K)

library(MCMCglmm)

# for MCMCglmm, the inverse matrix needs to be
# a sparse matrix
Kinv_sparse <- as(Kinv, 'sparseMatrix')

# simple spatial correlation
# ginverse is a list of with names matching
# a matrix to a random effect 
m <- MCMCglmm(le ~ 1, random=~ISO, data=geo, ginverse=list(ISO=Kinv_sparse))