if (FALSE) {
########################################
##
## example using NMs.randomShiftRotation
## first load the data:
data(puechcirc)
data(puechabonsp)
map <- puechabonsp$map
## Consider the first animal
## on an elevation map
anim1 <- puechcirc[1]
plot(anim1, spixdf=map[,1])
## We define a very simple treatment function
## for a NMs model: it just plots the randomized trajectory
## over the study area
## As required, the function takes two arguments:
## x is a data.frame storing a randomized trajectory (three
## columns: the x, y coordinates and the date)
## par contains the map of the study area
myfunc <- function(x, par)
{
par(mar = c(0,0,0,0))
## first plot the map
image(par)
## then add the trajectory
lines(x[,1], x[,2], lwd=2)
}
## Then we define the null model
##
## We define the range of the study area where the trajectory
## will be shifted:
rxy <- apply(coordinates(map),2,range)
rxy
## We define the null model with 9 repetitions
nmo <- NMs.randomShiftRotation(na.omit(anim1), rshift = TRUE, rrot = TRUE,
rx = rxy[,1], ry = rxy[,2], treatment.func = myfunc,
treatment.par = map, nrep=9)
## Then apply the null model
par(mfrow = c(3,3))
tmp <- testNM(nmo)
## You may try variations, by setting rshift or rrot to FALSE, to see
## the differences
## Note that some of the randomized trajectories are located outside the
## study area, although all barycentres are located within the X and Y
## limits of this study area.
## We may define a constraint function returning TRUE only if all
## relocations are located within the study area
## again, note that the two parameters are x and par
consfun <- function(x, par)
{
## first convert x to the class SpatialPointsDataFrame
coordinates(x) <- x[,1:2]
## then use the function over from the package sp
## to check whether all points in x are located inside
## the study area
ov <- over(x, geometry(map))
return(all(!is.na(ov)))
}
## Now fit again the null model under these constraints:
nmo2 <- NMs.randomShiftRotation(na.omit(anim1), rshift = TRUE, rrot = TRUE,
rx = rxy[,1], ry = rxy[,2], treatment.func = myfunc,
treatment.par = map,
constraint.func = consfun,
constraint.par = map, nrep=9)
## Then apply the null model
par(mfrow = c(3,3))
tmp <- testNM(nmo2)
## all the relocations are now inside the study area.
########################################
##
## example using NMs.randomCRW
## We generate correlated random walks with the same starting
## point as the original trajectory, the same turning angle
## distribution, and the same distance between relocation
## distribution. We use the same constraint function as previously
## (all relocations falling within the study area), and we
## use the same treatment function as previously (just plot
## the result).
mo <- NMs.randomCRW(na.omit(anim1), rangles=TRUE, rdist=TRUE,
treatment.func = myfunc,
treatment.par = map, constraint.func=consfun,
constraint.par = map, nrep=9)
par(mfrow = c(3,3))
tmp <- testNM(mo)
## Now, try a different treatment function: e.g. measure
## the distance between the first and last relocation,
## to test whether the animal is performing a return trip
myfunc2 <- function(x, par)
{
sqrt(sum(unlist(x[1,1:2] - x[nrow(x),1:2])^2))
}
## Now fit again the null model with this new treatment and 499 repetitions:
mo2 <- NMs.randomCRW(na.omit(anim1), rangles=TRUE, rdist=TRUE,
treatment.func = myfunc2,
treatment.par = map, constraint.func=consfun,
constraint.par = map, nrep=499)
## Then apply the null model
suppressWarnings(RNGversion("3.5.0"))
set.seed(298) ## to make the calculation reproducible
rand <- testNM(mo2)
## rand is a list with one element (there is one trajectory in anim1).
length(rand[[1]])
## The first element of rand is a list of length 499 (there are 499
## randomizations).
head(rand[[1]])
## unlist this list:
rand2 <- unlist(rand[[1]])
## calculate the observed average elevation:
obs <- myfunc2(na.omit(anim1)[[1]][,1:3], map)
## and performs a randomization test:
(rt <- as.randtest(rand2, obs, alter="less"))
plot(rt)
## Comparing to a model where the animal is moving randomly, and based
## on the chosen criterion (distance between the first and last
## relocation), we can see that the distance between the first and last
## relocation is rarely observed. It seems to indicate that the animal
## tends to perform a loop.
########################################
##
## example using NMs2NMm
## Given the previous results, we may try to see if all
## the trajectories in puechcirc are characterized by return
## trips
## We need a NMm approach. Because we have 3 burst in puechcirc
## we need a summary criterion. For example, the mean
## distance between the first and last relocation.
## We program a treatment function: it also takes two arguments, but x
## is now an object of class "ltraj" !
## par is needed, but will not be used in the function
myfunm <- function(x, par)
{
di <- unlist(lapply(x, function(y) {
sqrt(sum(unlist(y[1,1:2] - y[nrow(y),1:2])^2))
}))
return(mean(di))
}
## Now, prepare the NMs object: we do not indicate any treatment
## function (it would not be taken into account when NMs would be
## transformed to NMm). However, we keep the constraint function
## the simulated trajectories should fall within the study area
mo2s <- NMs.randomCRW(na.omit(puechcirc), constraint.func=consfun,
constraint.par = map)
## We convert this object to NMm, and we pass the treatment function
mo2m <- NMs2NMm(mo2s, treatment.func = myfunm, nrep=499)
## and we fit the model
suppressWarnings(RNGversion("3.5.0"))
set.seed(908)
resu <- testNM(mo2m)
## We calculate the observed mean distance between the
## first and last relocation
obs <- myfunm(na.omit(puechcirc))
## and performs a randomization test:
(rt <- as.randtest(unlist(resu), obs, alter="less"))
plot(rt)
## The test is no longer significant
########################################
##
## example using NMs.randomCs
## Consider this sample of 5 capture sites:
cs <- list(c(701184, 3161020), c(700164, 3160473),
c(698797, 3159908), c(699034, 3158559),
c(701020, 3159489))
image(map)
lapply(cs, function(x) points(x[1], x[2], pch=16))
## Consider this sample of distances:
dist <- c(100, 200, 150)
## change the treatment function so that the capture sites are showed as
## well. Now, par is a list with two elements: the first one is the map
## and the second one is the list of capture sites
myfunc <- function(x, par)
{
par(mar = c(0,0,0,0))
## first plot the map
image(par[[1]])
lapply(par[[2]], function(x) points(x[1], x[2], pch=16))
## then add the trajectory
lines(x[,1], x[,2], lwd=2)
}
## Now define the null model, with the same constraints
## and treatment as before
mod <- NMs.randomCs(na.omit(anim1), newCs=cs, newDistances=dist,
treatment.func=myfunc, treatment.par=list(map, cs),
constraint.func=consfun, constraint.par=map,
nrep=9)
## apply the null model
par(mfrow = c(3,3))
tmp <- testNM(mod)
}
Run the code above in your browser using DataLab