R/methods.r

In b-steve/ascr: Acoustic spatial capture-recapture

#' Extract ascr model coefficients
#'
#' Extracts estimated and derived parameters from a model fitted using
#' \link{fit.ascr}.
#'
#' @param object A fitted model from \link[ascr]{fit.ascr}.
#' @param pars A character vector containing either parameter names,
#'     or a subset of \code{"all"}, \code{"derived"}, \code{"fitted"},
#'     and \code{"linked"}; \code{"fitted"} corresponds to the
#'     parameters of interest, \code{"derived"} corresponds to
#'     quantities that are functions of these parameters (e.g., the
#'     effective survey area or animal density from an acoustic
#'     survey), and \code{"linked"} corresponds to the parameters AD
#'     Model Builder has maximised the likelihood over.
#' @param ... Other parameters (for S3 generic compatibility).
#'
#' @examples
#' coef(example.data$fits$simple.hn)
#' coef(example.data$fits$simple.hn, pars = "all")
#' coef(example.data$fits$simple.hn, pars = "derived")
#'
#' @export
coef.ascr <- function(object, pars = "fitted", ...){
    if ("all" %in% pars){
        pars <- c("fitted", "derived", "linked")
    }
    if (any(pars == "esa")){
        pars <- pars[-which(pars == "esa")]
        pars <- c(pars, paste("esa", 1:object$n.sessions, sep = "."))
    }
    par.names <- names(object$coefficients)
    if (!all(pars %in% c("fitted", "derived", "linked", "esa", par.names))){
        stop("Argument 'pars' must either contain a vector of parameter names, or a subset of \"fitted\", \"derived\", \"linked\", and \"all\".")
    }
    if (any(c("fitted", "derived", "linked") %in% pars)){
        which.linked <- grep("_link", par.names)
        linked <- object$coefficients[which.linked]
        which.derived <- which(substr(par.names, 1, 3) == "esa" | par.names == "Da")
        derived <- object$coefficients[which.derived]
        fitted <- object$coefficients[-c(which.linked, which.derived)]
        out <- mget(pars)
        names(out) <- NULL
        out <- c(out, recursive = TRUE)
        if (!object$fit.ihd){
            out <- out[c("D", names(out)[!(names(out) %in% c("D.(Intercept)", "D"))])]
        }
    } else {
        out <- object$coefficients[pars]
    }
    out
}

#' @rdname coef.ascr
#'
#' @param correct.bias Logical, if \code{TRUE}, estimated biases are
#'     subtracted from estimated parameter values.
#'
#' @export
coef.ascr.boot <- function(object, pars = "fitted",
                           correct.bias = FALSE, ...){
    out <- coef.ascr(object, pars)
    if (correct.bias){
        out <- out - get.bias(object, pars)
    }
    out
}
#' Extract the variance-covariance matrix from an ascr model
#' object
#'
#' Extracts the variance-covariance matrix for parameters in a model
#' fitted using \link[ascr]{fit.ascr}.
#'
#' @inheritParams coef.ascr
#'
#' @examples
#' vcov(example.data$fits$simple.hn)
#' vcov(example.data$fits$simple.hn, pars = "all")
#' vcov(example.data$fits$simple.hn, pars = "derived")
#'
#' @export
vcov.ascr <- function(object, pars = "fitted", ...){
    if ("all" %in% pars){
        pars <- c("fitted", "derived", "linked")
    }
    if (any(pars == "esa")){
        pars <- pars[-which(pars == "esa")]
        pars <- c(pars, paste("esa", 1:object$n.sessions, sep = "."))
    }
    par.names <- names(object$coefficients)
    if (!all(pars %in% c("fitted", "derived", "linked", "esa", par.names))){
        stop("Argument 'pars' must either contain a vector of parameter names, or a subset of \"fitted\", \"derived\", \"linked\", and \"all\".")
    }
    if (any(c("fitted", "derived", "linked") %in% pars)){
        which.linked <- grep("_link", par.names)
        which.derived <- which(substr(par.names, 1, 3) == "esa" | par.names == "Da")
        which.fitted <- (1:length(par.names))[-c(which.linked, which.derived)]
        keep <- NULL
        if ("fitted" %in% pars){
            keep <- c(keep, which.fitted)
        }
        if ("derived" %in% pars){
            keep <- c(keep, which.derived)
        }
        if ("linked" %in% pars){
            keep <- c(keep, which.linked)
        }
    } else {
        keep <- pars
    }
    if (!object$fit.ihd){
        keep <- keep[par.names[keep] != "D.(Intercept)"]
        keep <- keep[c(which(par.names[keep] == "D"), which(par.names[keep] != "D"))]
    }
    object$vcov[keep, keep, drop = FALSE]
}

#' Extract the variance-covariance matrix from a bootstrapped ascr
#' model object
#'
#' Extracts the variance-covariance matrix for parameters in a model
#' fitted using \link[ascr]{fit.ascr}, with a bootstrap procedure
#' carried out using \link[ascr]{boot.ascr}.
#'
#' @inheritParams coef.ascr
#'
#' @export
vcov.ascr.boot <- function(object, pars = "fitted", ...){
    if ("all" %in% pars){
        pars <- c("fitted", "derived", "linked")
    }
    if (any(pars == "esa")){
        pars <- pars[-which(pars == "esa")]
        pars <- c(pars, paste("esa", 1:object$n.sessions, sep = "."))
    }
    par.names <- names(object$coefficients)
    if (!all(pars %in% c("fitted", "derived", "linked", "esa", par.names))){
        stop("Argument 'pars' must either contain a vector of parameter names, or a subset of \"fitted\", \"derived\", \"linked\", and \"all\".")
    }
    if (any(c("fitted", "derived", "linked") %in% pars)){
        which.linked <- grep("_link", par.names)
        which.derived <- which(substr(par.names, 1, 3) == "esa" | par.names == "Da")
        which.fitted <- (1:length(par.names))[-c(which.linked, which.derived)]
        keep <- NULL
        if ("fitted" %in% pars){
            keep <- c(keep, which.fitted)
        }
        if ("derived" %in% pars){
            keep <- c(keep, which.derived)
        }
        if ("linked" %in% pars){
            keep <- c(keep, which.linked)
        }
    } else {
        keep <- pars
    }
    if (!object$fit.ihd){
        keep <- keep[par.names[keep] != "D.(Intercept)"]
        keep <- keep[c(which(par.names[keep] == "D"), which(par.names[keep] != "D"))]
    }
    object$boot$vcov[keep, keep, drop = FALSE]
}

#' Extract standard errors from an ascr model fit
#'
#' Extracts standard errors for estimated and derived parameters from
#' a model fitted using \link[ascr]{fit.ascr}.
#'
#' @inheritParams coef.ascr
#'
#' @examples
#' stdEr(example.data$fits$simple.hn)
#' stdEr(example.data$fits$simple.hn, pars = "all")
#' stdEr(example.data$fits$simple.hn, pars = "derived")
#'
#' @export
stdEr.ascr <- function(object, pars = "fitted", ...){
    if ("all" %in% pars){
        pars <- c("fitted", "derived", "linked")
    }
    if (any(pars == "esa")){
        pars <- pars[-which(pars == "esa")]
        pars <- c(pars, paste("esa", 1:object$n.sessions, sep = "."))
    }
    par.names <- names(object$coefficients)
    if (!all(pars %in% c("fitted", "derived", "linked", "esa", par.names))){
        stop("Argument 'pars' must either contain a vector of parameter names, or a subset of \"fitted\", \"derived\", \"linked\", and \"all\".")
    }
    if (any(c("fitted", "derived", "linked") %in% pars)){
        which.linked <- grep("_link", par.names)
        linked <- object$se[which.linked]
        which.derived <- which(substr(par.names, 1, 3) == "esa" | par.names == "Da")
        derived <- object$se[which.derived]
        fitted <- object$se[-c(which.linked, which.derived)]
        out <- mget(pars)
        names(out) <- NULL
        out <- c(out, recursive = TRUE)
        if (!object$fit.ihd){
            out <- out[c("D", names(out)[!(names(out) %in% c("D.(Intercept)", "D"))])]
        }
    } else {
        out <- object$se[pars]
    }
    out
}

#' Extract standard errors from a bootstrapped ascr model object
#'
#' Extracts standard errors for parameters of a model fitted using
#' \link[ascr]{fit.ascr}, with a bootstrap procedure carried out
#' using \link[ascr]{boot.ascr}.
#'
#' @param mce Logical, if \code{TRUE} Monte Carlo error for the
#'     standard errors is also returned.
#' @inheritParams coef.ascr
#'
#' @export
stdEr.ascr.boot <- function(object, pars = "fitted", mce = FALSE, ...){
    if ("all" %in% pars){
        pars <- c("fitted", "derived", "linked")
    }
    if (any(pars == "esa")){
        pars <- pars[-which(pars == "esa")]
        pars <- c(pars, paste("esa", 1:object$n.sessions, sep = "."))
    }
    par.names <- names(object$coefficients)
    if (!all(pars %in% c("fitted", "derived", "linked", "esa", par.names))){
        stop("Argument 'pars' must either contain a vector of parameter names, or a subset of \"fitted\", \"derived\", \"linked\", and \"all\".")
    }
    mces <- get.mce(object, estimate = "se")
    if (any(c("fitted", "derived", "linked") %in% pars)){
        which.linked <- grep("_link", par.names)
        linked <- object$boot$se[which.linked]
        which.derived <- which(substr(par.names, 1, 3) == "esa" | par.names == "Da")
        derived <- object$boot$se[which.derived]
        fitted <- object$boot$se[-c(which.linked, which.derived)]
        out <- mget(pars)
        names(out) <- NULL
        out <- c(out, recursive = TRUE)
        if (!object$fit.ihd){
            out <- out[c("D", names(out)[!(names(out) %in% c("D.(Intercept)", "D"))])]
        }
    } else {
        out <- object$boot$se[pars]
    }
    if (mce){
        out.vec <- out
        out <- cbind(out.vec, mces[names(out)])
        rownames(out) <- names(out.vec)
        colnames(out) <- c("Std. Error", "MCE")
    }
    out
}

#' Extract AIC from an ascr model object
#'
#' Extracts the AIC from an ascr model object.
#'
#' If the model is based on an acoustic survey where there are
#' multiple calls per individual, then AIC should not be used for
#' model selection. This function therefore returns NA in this case.
#'
#' @inheritParams coef.ascr
#' @inheritParams stats::AIC
#'
#' @export
AIC.ascr <- function(object, ..., k = 2){
    if (object$fit.freqs){
        message("NOTE: Use of AIC for this model relies on independence between locations of calls from the same animal, which may not be appropriate")
    }
    deviance(object) + k*length(coef(object))
}

#' Summarising ascr model fits
#'
#' Provides a useful summary of the model fit.
#'
#' @inheritParams coef.ascr
#'
#' @export
summary.ascr <- function(object, ...){
    coefs <- coef(object, "fitted")
    derived <- coef(object, "derived")
    coefs.se <- stdEr(object, "fitted")
    derived.se <- stdEr(object, "derived")
    infotypes <- object$infotypes
    detfn <- object$args$detfn
    n.sessions <- object$n.sessions
    out <- list(coefs = coefs, derived = derived, coefs.se = coefs.se,
                derived.se = derived.se, infotypes = infotypes,
                detfn = detfn, n.sessions = n.sessions)
    class(out) <- c("summary.ascr", class(out))
    out
}

#' @export
print.summary.ascr <- function(x, ...){
    n.coefs <- length(x$coefs)
    n.derived <- length(x$derived)
    mat <- matrix(0, nrow = n.coefs + n.derived + 1, ncol = 2)
    mat[1:n.coefs, 1] <- c(x$coefs)
    mat[1:n.coefs, 2] <- c(x$coefs.se)
    mat[n.coefs + 1, ] <- NA
    mat[(n.coefs + 2):(n.coefs + n.derived + 1), ] <- c(x$derived, x$derived.se)
    rownames(mat) <- c(names(x$coefs), "---", names(x$derived))
    colnames(mat) <- c("Estimate", "Std. Error")
    detfn <- c(hn = "Halfnormal", hhn = "Hazard halfnormal", hr = "Hazard rate", th = "Threshold",
               lth = "Log-link threshold", ss = "Signal strength")[x$detfn]
    infotypes <- c(bearing = "Bearings", dist = "Distances", ss = "Signal strengths",
                   toa = "Times of arrival", mrds = "Exact locations")[x$infotypes]
    cat("Detection function:", detfn, "\n")
    cat("Information types: ")
    cat(infotypes, sep = ", ")
    cat("\n", "\n", "Parameters:", "\n")
    printCoefmat(mat, na.print = "")
}

#' Confidence intervals for ascr model parameters
#'
#' Computes confidence intervals for one or more parameters estimated
#' in an ascr model object.
#'
#' Options for the argument \code{method} are as follows:
#' \code{"default"} for intervals based on a normal approximation
#' using the calculated standard errors (for objects of class
#' \code{ascr.boot}, these standard errors are calculated from the
#' bootstrap procedure); \code{"default.bc"} is a bias-corrected
#' version of \code{default}, whereby the estimated bias is subtracted
#' from each confidence limit; \code{"basic"} for the so-called
#' "basic" bootstrap method; and \code{"percentile"} for intervals
#' calculated using the bootstrap percentile method (the latter three
#' are only available for objects of class \code{ascr.boot}; see
#' Davison & Hinkley, 1997, for details).
#'
#' For method \code{"default"} with objects of class
#' \code{ascr.boot}, the appropriateness of the normal
#' approximation can be evaluated by setting \code{qqnorm} to
#' \code{TRUE}. If this indicates a poor fit, set \code{linked} to
#' \code{TRUE} and evaluate the QQ plot to see if this yields an
#' improvement (see Davison & Hinkley, 1997, pp. 194, for details).
#'
#' @references Davison, A. C., and Hinkley, D. V. (1997)
#' \emph{Bootstrap methods and their application}. Cambridge:
#' Cambridge University Press.
#'
#' @param parm A character vector containing either parameter names,
#' or a subset of \code{"all"}, \code{"derived"}, \code{"fitted"}, and
#' \code{"linked"}; \code{"fitted"} corresponds to the parameters of
#' interest, \code{"derived"} corresponds to quantities that are
#' functions of these parameters (e.g., the effective survey area or
#' animal density from an acoustic survey), and \code{"linked"}
#' corresponds to the parameters AD Model Builder has maximised the
#' likelihood over.
#' @param linked Logical, if \code{TRUE}, intervals for fitted
#' parameters are calculated on their link scales, then transformed
#' back onto their "real" scales.
#' @inheritParams coef.ascr
#' @inheritParams stats::confint
#'
#' @export
confint.ascr <- function(object, parm = "fitted", level = 0.95, linked = FALSE, ...){
    if (!object$args$hess){
        stop("Standard errors not calculated; use boot.ascr() or refit with 'hess = TRUE', if appropriate.")
    }
    calc.cis(object, parm, level, method = "default", linked, qqplot = FALSE,
             boot = FALSE, ask = FALSE, ...)
}

#' @param method A character string specifying the method used to
#' calculate the confidence intervals. See 'Details' below.
#' @param qqplot Logical, if \code{TRUE} and \code{method} is
#' \code{"default"} then a normal QQ plot is plotted. The default
#' method is based on a normal approximation; this plot tests its
#' validity.
#' @param ask Logical, if \code{TRUE}, hitting return will show the
#' next plot.
#'
#' @rdname confint.ascr
#'
#' @export
confint.ascr.boot <- function(object, parm = "fitted", level = 0.95, method = "default",
                                  linked = FALSE, qqplot = FALSE, ask = TRUE, ...){
    calc.cis(object, parm, level, method, linked, qqplot, boot = TRUE, ask, ...)
}

calc.cis <- function(object, parm, level, method, linked, qqplot, boot, ask, ...){
    if (any(c("all", "derived", "fitted") %in% parm)){
        parm <- names(coef(object, pars = parm))
    }
    if (linked){
        fitted.names <- names(coef(object, "fitted")[parm])
        fitted.names <- fitted.names[fitted.names != "mu.rates"]
        linked.names <- paste(fitted.names, "_link", sep = "")
        linked.names[linked.names == "D_link"] <- "D.(Intercept)_link"
        link.parm <- linked.names[!(linked.names %in% parm)]
        all.parm <- c(parm, link.parm)
    } else {
        all.parm <- parm
    }
    if (any(all.parm == "esa")){
        all.parm <- all.parm[-which(all.parm == "esa")]
        all.parm <- c(all.parm, paste("esa", 1:object$n.sessions, sep = "."))
    }
    if (method == "default" | method == "default.bc"){
        mat <- cbind(coef(object, pars = "all")[all.parm],
                     stdEr(object, pars = "all")[all.parm])
        FUN.default <- function(x, level){
            x[1] + qnorm((1 - level)/2)*c(1, -1)*x[2]
        }
        out <- t(apply(mat, 1, FUN.default, level = level))
        if (method == "default.bc"){
            out <- out - get.bias(object, parm)
        }
        if (qqplot & boot){
            opar <- par(ask = ask)
            for (i in parm){
                if (linked){
                    if (i %in% fitted.names){
                        j <- linked.names[fitted.names == i]
                    }
                } else {
                    j <- i
                }
                qqnorm(object$boot$boots[, j], main = i)
                abline(mean(object$boot$boots[, j], na.rm = TRUE),
                       sd(object$boot$boots[, j], na.rm = TRUE))
            }
            par(opar)
        }
    } else if (method == "basic"){
        qs <- t(apply(object$boot$boots[, all.parm, drop = FALSE], 2, quantile,
                      probs = c((1 - level)/2, 1 - (1 - level)/2),
                      na.rm = TRUE))
        mat <- cbind(coef(object, pars = "all")[all.parm], qs)
        FUN.basic <- function(x){
            2*x[1] - c(x[3], x[2])
        }
        out <- t(apply(mat, 1, FUN.basic))
    } else if (method == "percentile"){
        out <- t(apply(object$boot$boots[, all.parm, drop = FALSE], 2, quantile,
                       probs = c((1 - level)/2, 1 - (1 - level)/2),
                       na.rm = TRUE))
    }
    if (linked){
        for (i in fitted.names){
            linked.name <- paste(i, "_link", sep = "")
            if (linked.name == "D_link"){
                linked.name <- "D.(Intercept)_link"
                object$par.unlinks[[i]] <- exp
            }
            out[i, ] <- object$par.unlinks[[i]](out[linked.name, ])
        }
        out <- out[parm, , drop = FALSE]
    }
    percs <- c(100*(1 - level)/2, 100*(1 - (1 - level)/2))
    colnames(out) <- paste(round(percs, 2), "%")
    out
}

#' Extract density estimates given a set of covariate values
#'
#' Extracts density estimates from a set of supplied parameter values.
#'
#' @param object A fitted model from \link[ascr]{fit.ascr}.
#' @param newdata An optional data frame in which to look for
#'     variables with which to estimate density. If omitted, the
#'     function will return the estimated densities at the mask
#'     points.
#' @param se.fit A switch indicating if standard errors are
#'     required. At present, this will only work if \code{newdata} is
#'     also provided.
#' @param use.log If \code{TRUE}, density estimates and standard
#'     errors (if calculated) are provided on the log scale.
#' @param set.zero Indices for effects to ignore. For example,
#'     \code{set.zero = c(1, 3)} will set the first (probably an
#'     intercept) and third terms in the linear predictor to zero when
#'     calculating the predictions.
#' @param ... Other parameters (for S3 generic compatibility).
#'
#' @export
predict.ascr <- function(object, newdata = NULL, se.fit = FALSE,
                         use.log = FALSE, set.zero = NULL, ...){
    if (is.null(newdata)){
        out <- object$D.mask
    } else {
        ## Filling in x and y if required.
        if (!any(colnames(newdata) == "y")){
            newdata <- data.frame(y = rep(NA, nrow(newdata)), newdata)
        }
        if (!any(colnames(newdata) == "x")){
            newdata <- data.frame(x = rep(NA, nrow(newdata)), newdata)
        }
        ## Scaling data.
        newdata.scaled <- object$scale.covs(newdata)
        ## Creating model matrix.
        mm <- predict(gam(G = object$fgam), newdata = newdata.scaled, type = "lpmatrix")
        if (!is.null(set.zero)){
            mm[, set.zero] <- 0
        }
        ## Calculated estimated density.
        out <- as.vector(exp(mm %*% get.par(object, object$D.betapars)))
        if (use.log){
            out <- log(out)
        }
        if (se.fit){
            ## Gotta implement the delta method.
            n.predict <- length(out)
            n.betapars <- length(object$D.betapars)
            est.vcov <- vcov(object)
            est.vcov <- est.vcov[rownames(est.vcov) != "mu.rates", colnames(est.vcov) != "mu.rates"]
            n.par <- nrow(est.vcov)
            jacobian <- matrix(0, nrow = n.predict, ncol = n.par)
            colnames(jacobian) <- colnames(est.vcov)
            for (i in 1:ncol(jacobian)){
                if (colnames(jacobian)[i] %in% object$D.betapars){
                    if (use.log){
                        jacobian[, i] <- mm[, which(colnames(jacobian)[i] == object$D.betapars)]
                    } else {
                        jacobian[, i] <- mm[, which(colnames(jacobian)[i] == object$D.betapars)]*out
                    }
                } else {
                    jacobian[, i] <- 0
                }
            }
            se <- sqrt(diag(jacobian %*% est.vcov %*% t(jacobian)))
            out <- cbind(out, se)
            colnames(out) <- c("predict", "se")
        }
    }
    out
}