Bootstrap for R/qtl2 scan1()

scan1boot.R

# bootstrap of scan1()

scan1boot <-
    function(genoprobs, pheno, kinship=NULL, addcovar=NULL, Xcovar=NULL, intcovar=NULL,
             weights=NULL, n_boot=1, ...)
{
    # to contain the results
    result <- vector("list", n_boot)

    # drop individuals that are not in common across all data
    ind2keep <- get_common_ids(genoprobs, kinship, addcovar, Xcovar, intcovar, weights,
                               complete.cases=TRUE)
    ind2keep <- get_common_ids(ind2keep, rownames(pheno)[rowSums(is.finite(pheno))>0])
    if(length(ind2keep)<=2) stop("Need > 2 individuals with data")

    # subset to the common individuals
    genoprobs <- subset(genoprobs, ind=ind2keep)
    pheno <- pheno[ind2keep,,drop=FALSE]
    if(!is.null(addcovar)) addcovar <- addcovar[ind2keep,,drop=FALSE]
    if(!is.null(Xcovar)) Xcovar <- Xcovar[ind2keep,,drop=FALSE]
    if(!is.null(intcovar)) intcovar <- intcovar[ind2keep,,drop=FALSE]
    if(!is.null(weights)) weights <- weights[ind2keep]

    n <- length(ind2keep)

    for(i in seq_len(n_boot)) {

        # sample which individuals to use
        boot_ind <- sample(n, replace=TRUE)
        # create individual ids for them
        ids <- as.character(seq_len(n))

        # go through each piece and
        gp <- genoprobs
        for(chr in seq_along(gp)) {
            gp[[chr]] <- gp[[chr]][boot_ind,,]
            rownames(gp[[chr]]) <- ids
        }
        phe <- pheno[boot_ind,];rownames(phe) <- ids
        ac <- xc <- ic <- w <- k <- NULL
        if(!is.null(addcovar)) {
            ac <- addcovar[boot_ind,]
            rownames(ac) <- ids
        }
        if(!is.null(Xcovar)) {
            xc <- Xcovar[boot_ind,]
            rownames(xc) <- ids
        }
        if(!is.null(intcovar)) {
            ic <- intcovar[boot_ind,]
            rownames(ic) <- ids
        }
        if(!is.null(weights)) {
            w <- weights[boot_ind]
            names(w) <- ids
        }
        if(!is.null(kinship)) {
            if(is.matrix(kinship)) {
                k <- kinship[boot_ind,boot_ind]
                dimnames(k) <- list(ids, ids)
            } else { # loco-type kinship
                k <- kinship
                for(chr in seq_along(k)) {
                    k[[chr]] <- k[[chr]][boot_ind, boot_ind]
                    dimnames(k[[chr]]) <- list(ids, ids)
                }
            }
        }

        result[[i]] <- scan1(genoprobs=gp, pheno=phe, kinship=k,
                             addcovar=ac, Xcovar=xc, intcovar=ic, weights=w, ...)

    }

    result
}

### example of use:
# iron <- read_cross2(system.file("extdata", "iron.zip", package="qtl2"))
# map <- insert_pseudomarkers(iron$gmap, step=1)
# pr <- calc_genoprob(iron, map, error_prob=0.002, map_function="c-f")
# Xcovar <- get_x_covar(iron)
# out <- scan1(pr, iron$pheno, Xcovar=Xcovar)
# z <- scan1boot(pr, iron$pheno, Xcovar=Xcovar, n_boot=10)