[R] nlme: spatial autocorrelation on a sphere

[Spencer Graves]

Hi, Dan:


       That's great.


       Do you have a package that would be a reasonable home for this 
solution?  If yes, I would encourage you to add functions similar to 
those described by Malcolm Fairbrother.  If it were mine, I think I 
might prefer to call rdist.earth directly rather than make a local copy, 
and name the new functions something like corStruct.earth, gls.earth, 
and lme.earth.  (If you haven't done this before but are willing to try, 
you may want to look at the two "contributed" documents on "Creating R 
Packages" on CRAN.)


       I've also had moderate luck writing to the maintainer of another 
contributed package that might logically accept things like this and 
offering to do the work if they would consider such a friendly addition 
to their package.  Some have said, "yes". When that happened, I felt the 
addition would likely be more likely to reach its target audience -- and 
not contribute another package to be kept separate on CRAN.  When 
they've said, "no", I've created another package.


       Thanks for pursuing this issue and publicizing the solution you 
found.
       Spencer


On 10/3/2012 4:34 PM, Dan Bebber wrote:
> 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
> ________________________________________
> From: Spencer Graves [spencer.graves at prodsyse.com]
> Sent: 01 October 2012 16:34
> To: David Winsemius
> Cc: Dan Bebber; r-help at r-project.org
> Subject: Re: [R] nlme: spatial autocorrelation on a sphere
>
> On 10/1/2012 12:38 AM, David Winsemius wrote:
>> <snip>
>> LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html
>>
>
>         What's LMCTVIFY?
>
>
>         LMGTFY for LMCTVIY produced only LMGTFY, at least for me. There's
> space for it on Wikipedia ("RTFM") after LMGTFY (which I found using
> LMGTFY Wikipedia).
>
>
>         Spencer
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.


-- 
Spencer Graves, PE, PhD
President and Chief Technology Officer
Structure Inspection and Monitoring, Inc.
751 Emerson Ct.
San José, CA 95126
ph:  408-655-4567
web:  www.structuremonitoring.com