[R-sig-ME] nlme: spatial autocorrelation on a sphere

[Malcolm Fairbrother]

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)
}




> 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