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

[Mathieu Basille]

Dear all,

The answer Tom gave was actually incorrect. Let me first explain why,
and then give a solution to the problem.

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. So that if you copy the $li of the first ENFA
object to the second, considering that the vector of presence $pr
remains the same, then you will have exactly the same predictions, which
is not exactly what we want.

Now, let's try to give a solution to this problem. 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.



Tomislav Hengl a écrit :
> Dear Alexandre,
> 
> I think that this should be doable.
> 
> Once you run the ENFA to estimate the habitat model, copy the values of new rasters (future habitat)
> and replace the original values in the "" object. Then simply re-run the predict.enfa method e.g.:
> 
> # prepare the dataset for ENFA:
>> beidata <- data2enfa(as.kasc(list(dem=import.asc("dem.asc"), grad=import.asc("grad.asc"),
> twi=import.asc("twi.asc"), achan=import.asc("achan.asc"))), bei.sub.pnt at coords)
> # run ENFA: 
>> enfa.bei <- enfa(dudi.pca(beidata$tab, scannf=FALSE), beidata$pr, scannf=FALSE, nf=2)
> # same area, new environmental conditions:
>> beidata.new <- data2enfa(as.kasc(list(dem=import.asc("dem_new.asc"),
> grad=import.asc("grad_new.asc"), twi=import.asc("twi_new.asc"), achan=import.asc("achan_new.asc"))),
> bei.sub.pnt at coords)  # needs to have the same mask!
>> enfa.bei.new <- enfa(dudi.pca(beidata.new$tab, scannf=FALSE), beidata.new$pr, scannf=FALSE, nf=2)
> # copy the model parameters fitted previously:
>> enfa.bei.new$m <- enfa.bei$m
>> enfa.bei.new$s <- enfa.bei$s
>> enfa.bei.new$lw <- enfa.bei$lw
>> enfa.bei.new$li <- enfa.bei$li
>> enfa.bei.new$co <- enfa.bei$co
>> enfa.bei.new$mar <- enfa.bei$mar
>> bei.dist.new <- predict(enfa.bei.new, beidata.new$index, beidata.new$attr)
> 
> ...something like this (I am not really sure which parameters do you have to copy, but the $tab
> data.frame certainly needs to be replaced with new predictors).
> 
> Let me know if it works.
> 
> cheers,
> 
> T. Hengl
> http://spatial-analyst.net/wiki/index.php?title=Species_Distribution_Modelling 
> 
> 
>> -----Original Message-----
>> From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
>> Of Alexandre VILLERS
>> Sent: Monday, September 07, 2009 9:31 AM
>> To: Aide R SIG GEO; r-sig-ecology at r-project.org
>> Subject: [R-sig-Geo] Producing ENFA derived suitability maps with adehabitat
>>
>> 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
>>
>> --
>> Alexandre Villers
>> PhD Student
>> Team "Biodiversity"
>> Centre d'Etudes Biologiques de Chizé-CNRS UPR1934
>> 79360 Beauvoir sur Niort
>>
>> Phone +33 (0)5 49 09 96 13
>> Fax   +33 (0)5 49 09 65 26


-- 

~$ 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.