ase-research.org : Mathieu Basille -- R functions for SSF

R functions for SSF

by Mathieu Basille.

Here is a set of R functions to analyse movement data. They are directed towards the computation of Step Selection Functions (SSF). Reference: Fortin D., Beyer H., Boyce M., Smith D., Duchesne T. & Mao J. 2005 Wolves influence elk movements: Behavior shapes a trophic cascade in Yellowstone National Park. Ecology 86: 1320-1330.

distLocs

distLocs <- function(from, to, dint = 7, x = "x", y = "y", date = "date")
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-03
{
    distLoc <- function(fromi, to, dint) {
        cond <- to[, date] > fromi[, date] - dint * 24 * 3600 &
            to[, date] <= fromi[, date]
        return(min(sqrt((fromi[, x] - to[cond, x])^2 + (fromi[,
            y] - to[cond, y])^2), na.rm = TRUE))
    }
    dl <- unlist(lapply(1:nrow(from), function(i) distLoc(fromi = from[i,
        ], to = to, dint = dint)))
    is.na(dl) <- which(dl == Inf)
    return(dl)
}

lincirc.cor.test computes a linear-circular correlation, and tests it if appropriate (See Mardia, K. V. & Jupp, P. E. 2000 Directional statistics. Wiley, 429 pp. for more details).

x is a linear variable
theta is a circular variable, should be in radians (rad = deg*pi/180)
test allows to turn off or on the F-test, if x, and theta are independent and x is normally distributed.

lincirc.cor.test <- function(x, theta, test = TRUE)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-26
{
    if (length(x) != length(theta))
        stop("x and theta should be vector of same length.")
    if (max(theta, na.rm = TRUE) > 2 * pi | min(theta, na.rm = TRUE) <
        -2 * pi)
        warning("theta should be in radians. Calculations go on.")
    nna <- !is.na(x) & !is.na(theta)
    x <- x[nna]
    theta <- theta[nna]
    rxc <- cor(x, cos(theta))
    rxs <- cor(x, sin(theta))
    rcs <- cor(cos(theta), sin(theta))
    R2xt <- (rxc^2 + rxs^2 - 2 * rxc * rxs * rcs)/(1 - rcs^2)
    if (test) {
        stat <- (length(x) - 3) * R2xt/(1 - R2xt)
        pv <- 1 - pf(stat, 2, length(x) - 3)
        return(data.frame(cor = R2xt, n = length(x), `p-value` = pv))
    }
    else return(data.frame(cor = R2xt, n = length(x)))
}

rdSteps computes random steps. It takes a data frame as input, mandatory fields are id, x/ystart, rel./abs.angle, and dist.

distMax allows to use a maximum distance in the empirical distributions of distances used for random steps
other allows to draw random steps from distributions of all other individuals. If they don't exist, two new columns xend and yend are created with the end of the random step coordinates computed based on dx/dy.
A new column case is created (or overwritten) which takes the value 1 for observed steps and 0 for random ones.
The protect argument is to copy the value of the set of protected variables from the observed step to the random ones. For example, use protect = c("Area", "Sex") to copy the value of Area and Sex.

rdSteps <- function(df, emp = df, nr = 10, distMax = Inf, other = TRUE,
    id = "id", xstart = "x", ystart = "y", date = "date", dx = "dx",
    dy = "dy", dist = "dist", dt = "dt", abs.angle = "abs.angle",
    rel.angle = "rel.angle", xend = "xend", yend = "yend", case = "case",
    protect = NULL)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-10-26
{
    if (!exists(xend, df))
        df[, xend] <- df[, xstart] + df[, dx]
    if (!exists(yend, df))
        df[, yend] <- df[, ystart] + df[, dy]
    df[, case] <- 1
    angles <- na.omit(emp[, rel.angle])
    idA <- emp[!is.na(emp[, rel.angle]), id]
    dists <- na.omit(emp[emp[, dist] <= distMax, dist])
    idD <- emp[!is.na(emp[, dist]) & emp[, dist] <= distMax,
        id]
    rdStepsId <- function(ldfk) {
        idk <- as.character(ldfk[1, id])
        if (other)
            return(do.call("rbind", lapply(1:nrow(ldfk), function(i) rdStep(ldfk[i,
                ], anglesk = angles[idA != idk], distsk = dists[idD !=
                idk]))))
        else return(do.call("rbind", lapply(1:nrow(ldfk), function(i) rdStep(ldfk[i,
            ], anglesk = angles, distsk = dists))))
    }
    rdStep <- function(pt, anglesk = anglesk, distsk = distsk) {
        if (is.na(pt[, xstart]) | is.na(pt[, rel.angle]))
            return()
        else {
            rhord <- sample(anglesk, nr, replace = TRUE)
            alphard <- pt[, abs.angle] - pt[, rel.angle] + rhord
            distrd <- sample(distsk, nr, replace = TRUE)
            rd <- pt
            rd[1, ] <- NA
            rd[, id] <- pt[1, id]
            rd <- rd[rep(1, nr), ]
            rd[, xstart] <- pt[1, xstart]
            rd[, ystart] <- pt[1, ystart]
            rd[, date] <- pt[1, date]
            rd[, dx] <- cos(alphard * 180/pi) * distrd
            rd[, dy] <- sin(alphard * 180/pi) * distrd
            rd[, dist] <- distrd
            rd[, dt] <- pt[1, dt]
            rd[, abs.angle] <- alphard
            rd[, rel.angle] <- rhord
            rd[, xend] <- rd[, xstart] + rd[, dx]
            rd[, yend] <- rd[, ystart] + rd[, dy]
            rd[, case] <- 0
            if (!is.null(protect))
                rd[, protect] <- pt[1, protect]
            return(rbind(pt, rd))
        }
    }
    ldf <- split(df, f = df[, id])
    return(do.call("rbind", lapply(ldf, rdStepsId)))
}

steps2sldf converts a data frame of steps into a SpatialLinesDataFrame, with one step per Lines. It takes a data frame as input, mandatory fields are x/ystart. If x/yend are not provided (set the arguments to NULL), dx/y are used to generate x/yend. To export directly into a shapefile, use (this requires maptools):

writeLinesShape(locs2sldf(my_df, ...), fn = "my_steps")

Note: POSIXt format are not allowed by writeLinesShape and need to be converted to characters.

steps2sldf <- function(df, xstart = "x", ystart = "y", xend = "xend",
    yend = "yend", dx = "dx", dy = "dy")
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-09-08
{
    require(sp)
    if (!exists(xend, df))
        df[, xend] <- df[, xstart] + df[, dx]
    if (!exists(yend, df))
        df[, yend] <- df[, ystart] + df[, dy]
    df <- df[!is.na(df[, xend]), ]
    row.names(df) <- 1:nrow(df)
    return(SpatialLinesDataFrame(SpatialLines(lapply(1:nrow(df),
        function(i) Lines(Line(matrix(unlist(df[i, c(xstart,
            ystart, xend, yend)]), ncol = 2, byrow = TRUE)),
            ID = i))), df))
}

acftest

acftest <- function(residuals, id, type = c("correlation", "covariance",
    "partial"), ci = 0.95)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-12-19
{
    type <- match.arg(type)
    acfk <- lapply(levels(id), function(x) acf(residuals[id ==
        x], type = type, plot = FALSE))
    threshold <- lapply(acfk, function(x) qnorm((1 + ci)/2)/sqrt(x$n.used))
    lag <- unlist(lapply(1:length(acfk), function(i) which(acfk[[i]]$acf <
        threshold[[i]])[1]))
    return(list(acfk = acfk, threshold = threshold, lag = lag))
}

makeClusterSeq

makeClusterSeq <- function(seq, strata, case, nloc, na)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-12-19
{
    pat <- rep(c(0, NA), times = c(nloc, na))
    compt <- 0
    clust <- numeric()
    for (i in unique(seq)) {
        tabi <- table(strata[seq == i])
        maxi <- ceiling(sum(seq[case == 1] == i)/(nloc + na))
        if (maxi > 1)
            clusti <- (rep(pat, times = maxi) + rep(1:maxi +
                compt, each = (nloc + na)))[1:sum(seq[case ==
                1] == i)]
        else clusti <- rep(compt + 1, sum(seq[case == 1] == i))
        clusti <- rep(clusti, times = tabi)
        clust <- c(clust, clusti)
        compt <- compt + maxi
    }
    return(clust)
}

QIC

QIC <- function(mod, ...)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2010-11-10
### With the help of Thierry Duchesne and Marc Mazerolle
    UseMethod("QIC")

QIC.coxph <- function(mod, details = FALSE)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2010-11-10
### With the help of Thierry Duchesne and Marc Mazerolle
{
    trace <- sum(diag(solve(mod$naive.var) %*% mod$var))
    quasi <- mod$loglik[2]
    if (details)
        return(data.frame(QICR = -2 * quasi + 2 * trace, QICI = -2 *
            quasi + 2 * length(mod$coefficients), QuasiLL = quasi,
            K = length(mod$coefficients), Trace = trace))
    else return(-2 * quasi + 2 * trace)
}

modWeight

modWeight <- function(listMod, criterion = QIC)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-12-19
{
    CritMod <- unlist(lapply(listMod, criterion))
    dCritMod <- CritMod - min(CritMod)
    return(exp(-0.5 * dCritMod)/sum(exp(-0.5 * dCritMod)))
}

confintMod

confintMod <- function(object, parm, whichfac = NULL, level = 0.95,
    plot = FALSE, mar = c(5, 7, 2, 2) + 0.1, col, ...)
### Mathieu Basille, basille@ase-research.org
### Last modified: 2011-12-19
{
    cf <- coef(object)
    pnames <- names(cf)
    if (missing(parm))
        parm <- pnames
    else if (is.numeric(parm))
        parm <- pnames[parm]
    parmt <- names(cf) %in% parm
    a <- (1 - level)/2
    a <- c(a, 1 - a)
    pct <- stats:::format.perc(a, 3)
    fac <- qnorm(a)
    ci <- array(NA, dim = c(length(parm), 2L), dimnames = list(parm,
        pct))
    ses <- sqrt(diag(vcov(object)))[parm]
    ci[] <- cf[parm] + ses %o% fac
    ci <- data.frame(ci, cf[parm])
    names(ci) <- c(pct, "coef")
    if (missing(col))
        col <- c(rep("black", length(cf[parmt]) - length(whichfac)),
            rep("blue", length(whichfac)))
    if (!is.null(whichfac)) {
        ordfac <- whichfac[order(cf[whichfac])]
        whichfact <- (1:nrow(ci)) %in% whichfac
        ordcont <- (1:nrow(ci))[!whichfact][order(cf[!whichfact])]
        ci <- ci[c(ordcont, ordfac), ]
    }
    else ci <- ci[order(cf[parmt]), ]
    if (plot) {
        par(mar = mar)
        plot(ci[, 3], 1:length(cf[parmt]), xlim = range(ci),
            axes = FALSE, xlab = paste("Confidence intervals, level = ",
                level * 100, "%", sep = ""), ylab = NA, col = col,
            ...)
        grid()
        for (i in 1:nrow(ci)) lines(x = c(ci[i, 1], ci[i, 2]),
            y = c(i, i))
        abline(v = 0)
        axis(2, at = 1:nrow(ci), labels = FALSE)
        mtext(side = 2, text = row.names(ci), at = 1:nrow(ci),
            col = col, line = 1, las = 1)
        axis(1)
        box()
    }
    return(ci[nrow(ci):1, ])
}

[last update: 2011/12/22]