nlme: spatial autocorrelation on a sphere

This message from Malcolm Fairbrother at R-sig-ME, who has a potential solution...

[Begin forwarded message]
I struggled with this a couple years ago, and the best solution I could find (or at least that I think I found) was to modify corStruct.R, gls.R, and lme.R, adding a new function "dist2" which can calculate distances using "rdist.earth" from the "fields" package. I re-compiled nlme with these modifications, to make dist2 load automatically with nlme (and "fields" loads when required). This makes rdist.earth a valid choice of metrics, for corSpatial and similar functions.
The only change to gls.R and lme.R is to replace:
       distance <- lapply(covar,
                          function(el, metric) dist(as.matrix(el), metric),
                          metric = metric)
with:
       distance <- lapply(covar,
                          function(el, metric) dist2(as.matrix(el), metric),
                          metric = metric)
In corStruct.R you just need to replace calls to "dist" with calls to "dist2", and add the code below to create the function dist2. Then the call looks like, for example:
lme(y ~ x, random = ~ 1 | group, correlation = corExp(form = ~ LONG + LAT | group, metric="rdist.earth"), data=dat)
There may be more elegant ways of doing this (and obviously so for more professional coders), but I believe this works. Maybe this will save you some time, Ben.
Cheers,
Malcolm

dist2 <- function (x, method = "euclidean", diag = FALSE, upper = FALSE,
     p = 2)
{
     if (!is.na(pmatch(method, "euclidian")))
         method <- "euclidean"
     METHODS <- c("euclidean", "maximum", "manhattan", "canberra",
         "binary", "minkowski", "rdist.earth")
     method <- pmatch(method, METHODS)
     if (is.na(method))
         stop("invalid distance method")
     if (method == -1)
         stop("ambiguous distance method")
     N <- nrow(x <- as.matrix(x))
     if (method == 7) {
         require(fields)
         d <- rdist.earth(x, miles=F, R=6371)[lower.tri(rdist.earth(x, miles=F, R=6371))]
         } else {
         d <- .C("R_distance", x = as.double(x), nr = N, nc = ncol(x),
           d = double(N * (N - 1)/2), diag = as.integer(FALSE),
           method = as.integer(method), p = as.double(p), DUP = FALSE,
           NAOK = TRUE, PACKAGE = "stats")$d
         }
     attr(d, "Size") <- N
     attr(d, "Labels") <- dimnames(x)[[1L]]
     attr(d, "Diag") <- diag
     attr(d, "Upper") <- upper
     attr(d, "method") <- METHODS[method]
     if (method == 6)
         attr(d, "p") <- p
     attr(d, "call") <- match.call()
     class(d) <- "dist"
     return(d)
}
[End forwarded message]

Date: Tue, 2 Oct 2012 23:40:24 +0000 (UTC)
From: Ben Bolker <bbolker at gmail.com>
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] nlme: spatial autocorrelation on a sphere

Dan Bebber <dbebber at ...> writes:

I would like to write a spatial correlation structure function for
nlme, that replicates the ability of such functions in the 'ramps'
package to calculate great circle distances, e.g. corRSpher(form = ~
lat + lon, metric = "haversine").
The authors of the ramps package, Smith et al. (2008) write "The
spatial correlation structures in nlme are not directly used because
they do not allow great circle distance, which is very commonly
needed for spatial data" so there may be a case for adding this
functionality to nlme.
Unfortunately, I have been unable to locate the code to do this
within nlme.  My question is: is this code accessible, and if so,
how do I access it to modify it and then use it in functions like
corSpher?

[snip snip snip]

  In my opinion, corRSpher *ought* to work -- this is exactly the
sort of situation that user-defined (or other-package-defined)
corStructs are supposed to be designed for.  I confirmed, however,
that it doesn't work for gls() [I was feeling a little bit lazy
making up my example data, so I used gls() instead of [n]lme()].
I know that other 'pseudo-spatial' corStructs such as corBrownian etc.
(in the ape package) *do* work, so I'm not sure exactly why
corRSpher doesn't -- a quick look suggests that the problem *might*
be that it inherits from 'corSpatial', and there's some hard-coded
stuff in there.
  If I get a chance I will look further, but I don't think you're
missing something obvious (unfortunately).

  Ben Bolker