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

[Tomislav Hengl]

Dear Mathieu,

Thanks for the clarification. I have underestimated the complexity, but I
admit - habitat analysis is not really 'my cup of tee' (maybe this
question was more suited for the Environmentrics R-sig).

As least I was correct about the whole concept - it is doable  ;)

I would be interested to learn from Alexandre how did the prediction of
future habitat work and does it satisfies their objectives.

T. Hengl
http://home.medewerker.uva.nl/t.hengl/


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