R help - Oct 2012 - nlme: spatial autocorrelation on a sphere

I have spatial data on a sphere (the Earth) for which I would like to run an gls
model assuming that the errors are autcorrelated, i.e. including a corSpatial
correlation in the model specification.

In this case the distance metric should be calculated on the sphere, therefore
metric =  "euclidean" in (for example) corSpher would be incorrect.

I would be grateful for help on how to write a new distance metric for the
corSpatial function.
I believe there are several ways that distances on a sphere can be calculated in
R, for example the "distMeeus" function in the geosphere library.
However, I have no idea how to write this into a corSpatial function.

The aim is to end up with a metric  = "sphere" option that calculates
great circle distances between points using latitude and longitude.

Many thanks,
Dan

On Sep 30, 2012, at 6:48 PM, Dan Bebber wrote:
> I have spatial data on a sphere (the Earth) for which I would like to run
an gls model assuming that the errors are autcorrelated, i.e. including a
corSpatial correlation in the model specification.
> 
> In this case the distance metric should be calculated on the sphere,
therefore metric =  "euclidean" in (for example) corSpher would be
incorrect.
> 
> I would be grateful for help on how to write a new distance metric for the
corSpatial function.
> I believe there are several ways that distances on a sphere can be
calculated in R, for example the "distMeeus" function in the geosphere
library. However, I have no idea how to write this into a corSpatial function.
> 
> The aim is to end up with a metric  = "sphere" option that
calculates great circle distances between points using latitude and longitude.
LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html

-- 

David Winsemius, MD
Alameda, CA, USA

Thanks, but the problem is quite specific and not addressed on the Spatial Data
taskview page.
Quite specifically, I would like to know how to edit corSpatial functions to
calculate great circle distances.
The Bayesian equivalent, georamps in the ramps package, is able to do this,
therefore I imagine it must be possible for nlme.

Dan
________________________________________
From: David Winsemius [dwinsemius at comcast.net]
Sent: 01 October 2012 08:38
To: Dan Bebber
Cc: r-help at r-project.org
Subject: Re: [R] nlme: spatial autocorrelation on a sphere

On Sep 30, 2012, at 6:48 PM, Dan Bebber wrote:
> I have spatial data on a sphere (the Earth) for which I would like to run
an gls model assuming that the errors are autcorrelated, i.e. including a
corSpatial correlation in the model specification.
>
> In this case the distance metric should be calculated on the sphere,
therefore metric =  "euclidean" in (for example) corSpher would be
incorrect.
>
> I would be grateful for help on how to write a new distance metric for the
corSpatial function.
> I believe there are several ways that distances on a sphere can be
calculated in R, for example the "distMeeus" function in the geosphere
library. However, I have no idea how to write this into a corSpatial function.
>
> The aim is to end up with a metric  = "sphere" option that
calculates great circle distances between points using latitude and longitude.
LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html

--

David Winsemius, MD
Alameda, CA, USA

1.  Have you consulted Pinheiro and Bates (2000) Mixed-Effects 
Modeling in S and S-Plus (Springer)?  It does not has much information 
on the cor* constructor functions, but it does have some.  "corExp" is
discussed on p. 238, and information on other cor* functions should be 
in pages near 238.


       2.  Have you tried subscribing and posting to 
"r-sig-mixed-models at r-project.org"?  That email list is devoted 
primarily to "lmer", the replacement for nlme.  However, with luck,
you
will get an answer there, including, (a) whether what you want is 
available for "lmer", and (b) how to use it either with lmer or nlme.


       Best Wishes,
       Spencer


On 10/1/2012 12:59 AM, Dan Bebber wrote:
> Thanks, but the problem is quite specific and not addressed on the Spatial
Data taskview page.
> Quite specifically, I would like to know how to edit corSpatial functions
to calculate great circle distances.
> The Bayesian equivalent, georamps in the ramps package, is able to do this,
therefore I imagine it must be possible for nlme.
>
> Dan
> ________________________________________
> From: David Winsemius [dwinsemius at comcast.net]
> Sent: 01 October 2012 08:38
> To: Dan Bebber
> Cc: r-help at r-project.org
> Subject: Re: [R] nlme: spatial autocorrelation on a sphere
>
> On Sep 30, 2012, at 6:48 PM, Dan Bebber wrote:
>
>> I have spatial data on a sphere (the Earth) for which I would like to
run an gls model assuming that the errors are autcorrelated, i.e. including a
corSpatial correlation in the model specification.
>>
>> In this case the distance metric should be calculated on the sphere,
therefore metric =  "euclidean" in (for example) corSpher would be
incorrect.
>>
>> I would be grateful for help on how to write a new distance metric for
the corSpatial function.
>> I believe there are several ways that distances on a sphere can be
calculated in R, for example the "distMeeus" function in the geosphere
library. However, I have no idea how to write this into a corSpatial function.
>>
>> The aim is to end up with a metric  = "sphere" option that
calculates great circle distances between points using latitude and longitude.
> LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html
>
> --
>
> David Winsemius, MD
> Alameda, CA, USA
>
>
>
> ______________________________________________
> 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


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

On Oct 1, 2012, at 8:34 AM, Spencer Graves wrote:
> On 10/1/2012 12:38 AM, David Winsemius wrote:
>> <snip>
> 
>> LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html
>> 
> 
> 
>      What's LMCTVIFY?
Please accept my apologies for attempting to be cute and I also apologize to Mr
Bebber for reading his request far too quickly. I was trying to use (C)RAN
(T)ask (V)iews as a verb.
>      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).
I don't think it warrants being enshrined. I might also  have written: 
require(sos); findFn("metric spherical latitude longitude")  might
produce. In this instance that approach found the ramps::corRSpher function,
which has a "haversine" metric, much  more quickly than did my
subsequent efforts with RSeek, which I conducted after I realized that Bebber
had already made a good faith effort at identifying resources.

-- 


David Winsemius, MD
Alameda, CA, USA

Hi Spencer,

thanks, I'll try R-sig-mixed-models.

I tried using the corRSpher function from ramps in gls but received the
following error:
> Error in recalc.corSpatial(object[[i]], conLin) : 
>   Unknown spatial correlation class
Will post back here if I am successful.

Dan
________________________________________
From: Spencer Graves [spencer.graves at prodsyse.com]
Sent: 02 October 2012 04:29
To: Dan Bebber
Cc: David Winsemius; r-help at r-project.org
Subject: Re: [R] nlme: spatial autocorrelation on a sphere

Hi, Dan:


On 10/1/2012 8:07 PM, Dan Bebber wrote:
> Please be assured, I really did "RTFM" (i.e. Pinheiro &
Bates) as well as all the online sources, and have accessed the codes of the
corStruct functions, but they do not reveal how the "metric" argument
is used. Therefore, I can't reprogram to include "haversine"
distance, as seen in 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.
>
> If anyone knows a way in to the relevant code, please let me know.

       Thanks for the clarification and for your persistence in this
issue.  Others clearly have been asking for an enhancement of this nature.


       I assume you've checked and confirmed that
"ramps::corRSpher", as
suggested by David Winsemius, will NOT do what you want.  If so, I
suggest you subscribe to "R-sig-mixed-models" (per r-project.org ->
"Mailing Lists"), and post your question, citing Smith et al., to
there.  When you do, I suggest you include cc:  Douglas Bates
<bates at stat.wisc.edu>.


       I'm sorry I can't help more, but I hope you will soon find a
solution to this problem.


       Best Wishes,
       Spencer
>
> Thanks
> Dan
>
> Reference
> Smith BJ, Yan J & Cowles MK 2008. Unified Geostatistical Modeling for
Data Fusion and Spatial Heteroskedasticity with R Package ramps. Journal of
Statistical Software 25(10)
>
>
> ________________________________________
> From: David Winsemius [dwinsemius at comcast.net]
> Sent: 01 October 2012 17:32
> To: Spencer Graves
> Cc: Dan Bebber; r-help at r-project.org
> Subject: Re: [R] nlme: spatial autocorrelation on a sphere
>
> On Oct 1, 2012, at 8:34 AM, Spencer Graves wrote:
>
>> On 10/1/2012 12:38 AM, David Winsemius wrote:
>>> <snip>
>>> LMCTVTFY: http://cran.r-project.org/web/views/Spatial.html
>>>
>>
>>       What's LMCTVIFY?
> Please accept my apologies for attempting to be cute and I also apologize
to Mr Bebber for reading his request far too quickly. I was trying to use (C)RAN
(T)ask (V)iews as a verb.
>
>>       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).
> I don't think it warrants being enshrined. I might also  have written: 
require(sos); findFn("metric spherical latitude longitude")  might
produce. In this instance that approach found the ramps::corRSpher function,
which has a "haversine" metric, much  more quickly than did my
subsequent efforts with RSeek, which I conducted after I realized that Bebber
had already made a good faith effort at identifying resources.
>
> --
>
>
> David Winsemius, MD
> Alameda, CA, USA
>

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


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