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