[R-sig-eco] Producing ENFA derived suitability maps with adehabitat

[Mathieu Basille]

Dear all,

[This is a slightly modified version of a message posted on R-SIG-Geo. I
post here too for the sake of consistency]

The so-called predictions of the ENFA (function 'predict.enfa') are
based on the row coordinates ($li) and the vector of presence ($pr).
>From the row coordinates, the function computes Mahalanobis distances
from the center of the niche to every pixel (row), given the niche
covariance structure.

The row coordinates gives you the coordinates of every pixel projected
into the new space created by the ENFA. You can compute them by hand by
multiplying the (scaled) table of environmental variables with the
columns coordinates. For example, following the example of the function
'enfa':

## We load the data
data(lynxjura)
map <- lynxjura$map
tmp <- lynxjura$locs[,4]!="D"
locs <- lynxjura$locs[tmp, c("X","Y")]
map[,4] <- sqrt(map[,4])

## We perform the ENFA
dataenfa1 <- data2enfa(map, locs)
pc <- dudi.pca(dataenfa1$tab, scannf = FALSE)
(enfa1 <- enfa(pc, dataenfa1$pr, scannf = FALSE))

## We compute the row coordinates by hand and check that this is
## correct
bla <- as.matrix(enfa1$tab) %*% as.matrix(enfa1$co)
all.equal(bla, as.matrix(enfa1$li))

This means that we can then reproject any new environmental variable
into the space created by the ENFA. Let's say that a second set of
environmental maps is called 'map2'. It can be then projected as follows:

li2 <- as.matrix(scale(kasc2df(map2)$tab)) %*% as.matrix(enfa1$co)

Note that:
1) the table must be scaled, as is the $tab component of an ENFA object
("modified array").
2) for the sake of simplicity, I used here the function 'scale', which
used a N-1 denominator. In the ENFA example, a PCA ('dudi.pca') is ran
first on the data, before we use the scale table. In the 'dudi.pca'
function, the denominator is N, so that it might give slightly different
results (negligible for large N).

Then, to "predict", you must use an adapted 'predict.enfa' function,
with a 'new' argument in which we feed the new maps:

predict.enfa <- function (object, index, attr, nf, new, ...)
{
    if (!inherits(object, "enfa"))
        stop("should be an object of class \"enfa\"")
    warning("the enfa is not mathematically optimal for prediction:\n
please consider the madifa instead")
    if ((missing(nf)) || (nf > object$nf))
        nf <- object$nf
    Zli <- object$li[, 1:(nf + 1)]
    f1 <- function(x) rep(x, object$pr)
    Sli <- apply(Zli, 2, f1)
    m <- apply(Sli, 2, mean)
    cov <- t(as.matrix(Sli)) %*% as.matrix(Sli)/nrow(Sli)
    if (!missing("new"))
        Zli <- as.matrix(dudi.pca(kasc2df(new)$tab, scannf = FALSE,
               nf = 1)$tab) %*% as.matrix(object$co)
    maha <- mahalanobis(Zli, center = m, cov = cov)
    map <- getkasc(df2kasc(data.frame(toto = maha, tutu = maha),
        index, attr), "toto")
    return(invisible(map))
}

And this should work like this:

pred2 <- predict(enfa1, dataenfa1$index, dataenfa1$attr, new = map2)

Note that the two maps have here the same attributes. It still works if
it's not the case, simply by adjusting the index and attr in the predict
call (given by a kasc2df on map2).

Hope this helps,
Mathieu.



Alexandre VILLERS a écrit :
> Dear list(s) 'members,
> 
> I have a dataset representing the position of a species over a large
> region for 3 different years (2000, 2004 and 2008) and seven
> ecogeographical variables. Given that the study takes place on an
> agricultural region, the landscape changes every year at fine scale and
> there is a long term trend in crops sowed, leading me to account for
> between years variability.
> I would like to use the niche determined by the landscape and location
> of birds in 2000 and then, appply this niche over the landscape in 2004
> and 2008 (this would, I believe, give me a first answer on whether
> changes in birds' location are due to a decrease in habitat suitability
> or birds).
> I have already computed ENFA with adehabitat (using the doc provided
> with the package "adehabitat" and the www.spatial-analyst.net of T.
> Hengl) but I don't see how exactly using the result of an ENFA on a
> "new" landscape...
> 
> Any link or help would be welcome.
> 
> Alex
> 

-- 

~$ whoami
Mathieu Basille, PhD

~$ locate
Laboratoire de Biométrie et Biologie Évolutive
Université Claude Bernard Lyon 1 - France
http://lbbe.univ-lyon1.fr/

~$ info
http://ase-research.org/basille

~$ fortune
``If you can't win by reason, go for volume.''
Calvin, by Bill Watterson.