[R-sig-ME] nlme: spatial autocorrelation on a sphere
Malcolm Fairbrother
m.fairbrother at bristol.ac.uk
Wed Oct 3 12:49:55 CEST 2012
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
More information about the R-sig-mixed-models
mailing list