inst/doc/permutations.R

In permute: Functions for Generating Restricted Permutations of Data

## ----setup, echo = FALSE, results = "hide", message = FALSE, cache = FALSE----
knitr::opts_chunk$set(message = FALSE, warning = FALSE)
library("permute")

## ----load_jackal--------------------------------------------------------------
library("permute")
data(jackal)
jackal

## ----ttest_jackal-------------------------------------------------------------
jack.t <- t.test(Length ~ Sex, data = jackal, var.equal = TRUE,
                 alternative = "greater")
jack.t

## ----meanFun------------------------------------------------------------------
meanDif <- function(x, grp) {
 mean(x[grp == "Male"]) - mean(x[grp == "Female"])
}

## ----randJackal---------------------------------------------------------------
Djackal <- numeric(length = 5000)
N <- nrow(jackal)
set.seed(42)
for(i in seq_len(length(Djackal) - 1)) {
    perm <- shuffle(N)
    Djackal[i] <- with(jackal, meanDif(Length, Sex[perm]))
}
Djackal[5000] <- with(jackal, meanDif(Length, Sex))

## ----hist_jackal, fig=FALSE, echo=TRUE, eval=FALSE----------------------------
# hist(Djackal, main = "",
#      xlab = expression("Mean difference (Male - Female) in mm"))
# rug(Djackal[5000], col = "red", lwd = 2)

## -----------------------------------------------------------------------------
(Dbig <- sum(Djackal >= Djackal[5000]))

## -----------------------------------------------------------------------------
Dbig / length(Djackal)

## ----draw_hist_jackal, fig=TRUE, echo=FALSE, fig.cap="Distribution of the difference of mean mandible length in random allocations, ten to each sex."----
hist(Djackal, main = "",
     xlab = expression("Mean difference (Male - Female) in mm"))
rug(Djackal[5000], col = "red", lwd = 2)

## -----------------------------------------------------------------------------
choose(20, 10)

## ----show_args----------------------------------------------------------------
args(shuffle)

## ----show_str-----------------------------------------------------------------
str(how())

## ----compare_shuffle_sample---------------------------------------------------
set.seed(2)
(r1 <- shuffle(10))
set.seed(2)
(r2 <- sample(1:10, 10, replace = FALSE))
all.equal(r1, r2)

## ----series1------------------------------------------------------------------
set.seed(4)
x <- 1:10
CTRL <- how(within = Within(type = "series"))
perm <- shuffle(10, control = CTRL)
perm
x[perm] ## equivalent

## ----grid1--------------------------------------------------------------------
set.seed(4)
plt <- gl(3, 9)
CTRL <- how(within = Within(type = "grid", ncol = 3, nrow = 3),
            plots = Plots(strata = plt))
perm <- shuffle(length(plt), control = CTRL)
perm

## ----vis_grid1, keep.source=TRUE----------------------------------------------
## Original
lapply(split(seq_along(plt), plt), matrix, ncol = 3)
## Shuffled
lapply(split(perm, plt), matrix, ncol = 3)

## ----grid_2, keep.source=TRUE-------------------------------------------------
set.seed(4)
CTRL <- how(within = Within(type = "grid", ncol = 3, nrow = 3,
                            constant = TRUE),
            plots = Plots(strata = plt))
perm2 <- shuffle(length(plt), control = CTRL)
lapply(split(perm2, plt), matrix, ncol = 3)

## ----series_2, results="hide"-------------------------------------------------
how(nperm = 10, within = Within(type = "series"))

## ----shuffleSet_1-------------------------------------------------------------
set.seed(4)
CTRL <- how(within = Within(type = "series"))
pset <- shuffleSet(10, nset = 5, control = CTRL)
pset

## ----results="hide"-----------------------------------------------------------
how(nperm = 999)

## ----withinArgs, echo=FALSE---------------------------------------------------
args(Within)

## ----ptest-fun----------------------------------------------------------------
pt.test <- function(x, group, nperm = 199) {
    ## mean difference function
    meanDif <- function(i, x, grp) {
        grp <- grp[i]
        mean(x[grp == "Male"]) - mean(x[grp == "Female"])
    }
    ## check x and group are of same length
    stopifnot(all.equal(length(x), length(group)))
    ## number of observations
    N <- nobs(x)
    ## generate the required set of permutations
    pset <- shuffleSet(N, nset = nperm)
    ## iterate over the set of permutations applying meanDif
    D <- apply(pset, 1, meanDif, x = x, grp = group)
    ## add on the observed mean difference
    D <- c(meanDif(seq_len(N), x, group), D)
    ## compute & return the p-value
    Ds <- sum(D >= D[1]) # how many >= to the observed diff?
    Ds / (nperm + 1)     # what proportion of perms is this (the pval)?
}

## ----run-ptest----------------------------------------------------------------
set.seed(42) ## same seed as earlier
pval <- with(jackal, pt.test(Length, Sex, nperm = 4999))
pval

## ----parallel-ptest-fun-------------------------------------------------------
ppt.test <- function(x, group, nperm = 199, cores = 2) {
    ## mean difference function
    meanDif <- function(i, .x, .grp) {
        .grp <- .grp[i]
        mean(.x[.grp == "Male"]) - mean(.x[.grp == "Female"])
    }
    ## check x and group are of same length
    stopifnot(all.equal(length(x), length(group)))
    ## number of observations
    N <- nobs(x)
    ## generate the required set of permutations
    pset <- shuffleSet(N, nset = nperm)
    if (cores > 1) {
        ## initiate a cluster
        cl <- makeCluster(cores)
        on.exit(stopCluster(cl = cl))
        ## iterate over the set of permutations applying meanDif
        D <- parRapply(cl, pset, meanDif, .x = x, .grp = group)
    } else {
        D <- apply(pset, 1, meanDif, .x = x, .grp = group)
    }
    ## add on the observed mean difference
    D <- c(meanDif(seq_len(N), x, group), D)
    ## compute & return the p-value
    Ds <- sum(D >= D[1]) # how many >= to the observed diff?
    Ds / (nperm + 1)     # what proportion of perms is this (the pval)?
}

## ----run-pptest---------------------------------------------------------------
require("parallel")
set.seed(42)
system.time(ppval <- ppt.test(jackal$Length, jackal$Sex, nperm = 9999,
                              cores = 2))
ppval

## ----run-pptest2--------------------------------------------------------------
set.seed(42)
system.time(ppval2 <- ppt.test(jackal$Length, jackal$Sex, nperm = 9999,
                               cores = 1))
ppval2

## ----get-set-eg0--------------------------------------------------------------
hh <- how()

## ----get-set-eg1--------------------------------------------------------------
getNperm(hh)

## ----<get-set-eg2-------------------------------------------------------------
getCall(hh)
setNperm(hh) <- 999
getNperm(hh)
getCall(hh)

## ----get-set-eg3--------------------------------------------------------------
hh <- how(within = Within(type = "series"),
          plots = Plots(type = "series", strata = gl(10, 5)),
          blocks = gl(5, 10))

## ----get-set-eg4--------------------------------------------------------------
pl <- getPlots(hh)
setType(pl) <- "free"
setPlots(hh) <- pl

## ----get-set-eg5--------------------------------------------------------------
getType(hh, which = "plots")

## ----get-set-eg6--------------------------------------------------------------
getCall(getPlots(hh))

## ----get-set-eg7--------------------------------------------------------------
hh <- update(hh, plots = update(getPlots(hh), type = "series"))
getType(hh, which = "plots")

## ----seesionInfo, echo=FALSE--------------------------------------------------
sessioninfo::session_info()