Dear Steve,
I am exploring the trajectory analysis functions of the new
adehabitatLT package and am coming across an error while I am trying
to partition the trajectory using the function partmod.ltraj.
Unfortunately I do not get the error on all data sets that I am
working on, which is making it very difficult for me to track down
exactly what the cause is. The error is "Error in if (prot > 1e-08) {
: missing value where TRUE/FALSE needed"
As an example, the code below shows the process for one particular
data set, the first completes successfully, while the second results
in the above error:
# Example 1
tested.means <- round(seq(0,20000, length=10),0)
limod <- as.list(paste("dnorm(dist, mean=", tested.means, ", sd=750)"))
(mod <- modpartltraj(tr3, limod))
bestpartmod(mod)
(pm <- partmod.ltraj(tr3,14, mod))
# Example 2
tested.means <- round(seq(0,15000, length=10),0)
limod <- as.list(paste("dnorm(dist, mean=", tested.means, ", sd=750)"))
(mod <- modpartltraj(tr3, limod))
bestpartmod(mod)
(pm <- partmod.ltraj(tr3,8, mod))
In the second example, I have only changed the upper bound of the set
of models. As a result, the recommended number of partitions
(bestpartmod) changes to 8, which I use to create 'pm' but the above
error is the result. If I continue to reduce the number of
partitions, it's finally successful at 3. Alternatively, if I reduce
the upper bound of the model from 20000 to 16000 and also reduce the
number of models from 10 to 8, it will also run successfully.
This is indeed a bug that occurs in partmod.ltraj when missing values
occur in the trajectory and na.manage="prop.move". It is now corrected.
Until I submit a revised version (probably by the end of the week), you
can use this corrected version of partmod.ltraj:
partmod.ltraj <- function (tr, npart, mods, na.manage = c("prop.move",
"locf"))
{
if (!inherits(tr, "ltraj"))
stop("tr should be of class \"ltraj\"")
if (length(tr) > 1)
stop("only one traject can be passed")
if (!inherits(mods, "modpartltraj"))
stop("mods should be of class modpartltraj")
na.manage <- match.arg(na.manage)
cor <- tr[[1]]
indiceNA <- attr(mods, "nna.places")
if (npart > nrow(mods))
stop("too large number of segments required")
toto <- .C("partrajr", as.double(t(as.matrix(mods))), double(npart),
integer(npart), integer(npart + 1), as.integer(nrow(mods)),
as.integer(ncol(mods)), as.integer(npart), PACKAGE =
"adehabitatLT")
curloc <- rev(toto[[4]])
curloc[2:length(curloc)] <- curloc[2:length(curloc)] + 1
curmod <- rev(toto[[3]])
curma <- rev(toto[[2]])
filo <- curloc[-length(curloc)]
lalo <- curloc[-1]
lalo[length(lalo)] <- nrow(cor)
resltr <- lapply(1:length(lalo), function(i) {
if (i == 1) {
xyt <- cor[1:indiceNA[lalo[i]], c("x", "y", "date")]
}
else {
if (i == length(lalo)) {
xyt <- cor[indiceNA[filo[i]]:nrow(cor), c("x",
"y", "date")]
}
else {
xyt <- cor[indiceNA[filo[i]]:indiceNA[lalo[i]],
c("x", "y", "date")]
}
}
return(as.ltraj(xyt[, c("x", "y")], xyt[, c("date")],
id = id(tr), burst = i))
})
cseq <- function(x) {
id <- diff(c(1, c(1:length(x))[abs(c(0, diff(x))) > 0],
length(x) + 1))
split(x, unlist(sapply(1:length(id), function(i) rep(i,
id[i]))))
}
if (na.manage == "prop.move") {
nadf <- do.call("rbind", lapply(1:length(resltr), function(i) {
nas <- is.na(resltr[[i]][[1]]$dist[-nrow(resltr[[i]][[1]])])
vec <- cseq(nas)
beg <- sum(vec[[length(vec)]])
intern <- sum(nas) - beg
return(c(beg, intern))
}))
nadf <- as.data.frame(nadf)
typmod <- tapply(nadf[, 2], factor(curmod), sum)
typmod <- typmod/sum(typmod)
for (i in 2:length(resltr)) {
gg <- resltr[[i - 1]][[1]]
gg2 <- resltr[[i]][[1]]
gg <- gg[-nrow(gg), ]
ff <- cseq(is.na(gg$dist))
nna <- sum(ff[[length(ff)]])
if (nna > 1) {
prot <- sum(typmod[names(typmod)%in%c(curmod[i - 1],
curmod[i])])
if (prot > 1e-08) {
nna1 <- floor(nna * typmod[names(typmod)==curmod[i -
1]]/prot)
}
else {
nna1 <- floor(nna/2)
}
nna2 <- nna - nna1
gg2 <- rbind(gg[(nrow(gg) - nna2):nrow(gg), ],
gg2)
gg <- gg[1:(nrow(gg) - nna2), ]
resltr[[i - 1]] <- as.ltraj(gg[, c("x", "y")],
gg[, c("date")], id = id(tr), burst = i - 1)
resltr[[i]] <- as.ltraj(gg2[, c("x", "y")], gg2[,
c("date")], id = id(tr), burst = i)
}
}
}
resltr <- do.call("c.ltraj", resltr)
resu <- list(ltraj = resltr, stats = list(locs = curloc,
Mk = curma, mod = curmod, which.mod = colnames(mods)[curmod]))
attr(resu, "nna.places") <- indiceNA
class(resu) <- "partltraj"
return(resu)
}
Thanks for reporting,
Best regards,
Clément Calenge.
If anyone has any experience with this or thoughts that might help, I
would appreciate hearing from you.
Thanks in advance
Steve
_______________________________________________
AniMov mailing list
AniMov@faunalia.it
http://lists.faunalia.it/cgi-bin/mailman/listinfo/animov
--
Clément CALENGE
Cellule d'appui à l'analyse de données
Direction des Etudes et de la Recherche
Office national de la chasse et de la faune sauvage
Saint Benoist - 78610 Auffargis
tel. (33) 01.30.46.54.14
_______________________________________________
AniMov mailing list
AniMov@faunalia.it
http://lists.faunalia.it/cgi-bin/mailman/listinfo/animov
