# R code for chapter 53 resampling: permutations, jacknives and bootstraps
# Book details:
# Health Science Statistics using R and R Commander
# by Robin Beaumont
# Paperback - ISBN 9781907904318
# http://www.amazon.co.uk/Health-Science-Statistics-using-Commander/dp/190790431X
# or direct from publisher at: 
# http://www.scionpublishing.com/
# book website: http://www.robin-beaumont.co.uk/rbook
###########################################################
###
### producing all possible combinations:
install.packages("gtools", dependencies=TRUE) 
update.packages(ask=FALSE) 
library (gtools)
results<- c('a','b','c','d','e','f','g','h','i','j')
results

combins <- combinations (10, 5 , v= results, set=FALSE, repeats=FALSE)
head(combins, 5)
tail(combins, 5)
length(combins)
dim(combins)
#####
# calculate summary count for each and produce a histogram:
combinN <- combins
combinN[combinN=='a']<-2 
combinN[combinN=='b']<- 3
combinN[combinN=='c']<- 1
combinN[combinN=='d']<- 4
combinN[combinN=='e']<- 3
combinN[combinN=='f']<- 5
combinN[combinN=='g']<- 6
combinN[combinN=='h']<- 4
combinN[combinN=='i']<- 7
combinN[combinN=='j']<- 6
# mode(combinN) <- 'numeric'  

combinN<- apply(combinN, 2, as.numeric)  
head(combinN)
resultsnum<- sum(c(2,3,1,4,3,5,6,4,7,6))
resultsnum
combinsum <- apply(combinN, 1, sum) 
combinsum 
statistic <-  combinsum - (resultsnum -  combinsum)
statistic
table(statistic)
hist(statistic)
##################
# example of obtaining a p value from above results:
4 x 0.003968254
#
# 53.1.1	Performing permutation tests in R - small samples
install.packages("coin", dependencies=TRUE) 
library(coin)
results<- c(2,3,1,4,3,5,6,4,7,6)
group <- factor(c(1,1,1,1,1,2,2,2,2,2))
rst<- wilcox_test(results ~ group, 
distribution="exact", conf.int=TRUE)
rst
##  Wilcoxon matched pairs
library(coin) 
post<-c(7, 6,8,7,6,9)
pre<-c(1,4,3,5,4,5) 
# you you can request a ci here
rst<- wilcoxsign_test(pre ~ post, distribution="exact")
rst
# to get the result:
statistic(rst)
# quartiles for ci:
qperm(rst, 0.025) 
qperm(rst, 0.975) 

# giving the 95% CI:
statistic(rst) + qperm(rst , 0.025) 
statistic(rst) + qperm(rst, 0.975) 
# gives:
# statistic(rst) + qperm(rst , 0.025) 
#     neg 
# -4.027992 
# > statistic(rst) + qperm(rst, 0.975) 
#       neg 
# -0.4239992 
#######
# oneway tests
mydataframe <- read.delim ("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/ttest_2independent_long1.txt", header=TRUE)
###
library(coin) 
model1<- oneway_test(values~ ind, distribution=approximate(B=10000), 
         data= mydataframe)
model1
# draw a histogram
hist(support(model1))
# asymptotic approach:
model2<- oneway_test(values ~ ind, data= mydataframe)  
model2
# quartiles for ci:
qperm(model2, 0.025) 
qperm(model2, 0.975)
# giving the 95% CI:
statistic(model2) + qperm(model2, 0.025) 
statistic(model2) + qperm(model2, 0.975)
# gives:
# > statistic(model2) + qperm(model2, 0.025) 
#  control 
# -3.953726 
# >statistic(model2) + qperm(model2, 0.975)
#  control 
# -0.03379762 
####
# 53.2	Bootstrapping 
index <- c(6, 11, 11, 1, 8, 14, 3, 7, 15, 3, 3, 7, 2, 4, 9)
data1 <- c(3,5,4,7,6,2,4,3,9,10,2,1,1,3,4)
data1[index]

install.packages("boot", dependencies=TRUE) 
library(boot)
samplemean <- function(thedata, indices) 
{
value <- mean(thedata[indices])
return (value)
}
data1 <- c(3,5,4,7,6,2,4,3,9,10,2,1,1,3,4)
bootmean <-  boot(data1, samplemean, R = 100)
bootmean

boot.ci(bootmean, type = "all")
###################
####### correlation example
library(boot)
boot_corr <- function (thedata, index)
{
thesample <- thedata[index,]
value <- cor(thesample$HBA1C , thesample$fastingBG)
return(value)
}
mydataframe <- read.delim("http://www.robin-beaumont.co.uk/virtualclassroom/book2data/regression1.dat", header=TRUE)
names(mydataframe)
#  "HBA1C"     "fastingBG"
attach(mydataframe)
cor(HBA1C, fastingBG)
bootcorr <-  boot(mydataframe, boot_corr, R = 1000)
bootcorr
plot(bootcorr)
boot.ci(bootcorr, type = "all")
###########
# 53.3	The Jackknife
jack.after.boot(bootmean, index = 1)
###########
# 53.4	Memory management 
memory.size(FALSE)
memory.size(TRUE)
memory.limit()
######### random seed - 626 is the default value
# to find out retrospectively what the last random seed number was
# use summary(bootobject) e.g:
summary(bootcorr)

save.seed <- .Random.seed  # to save your current random seed number
set.seed(number)  # to set the random seed number to number 


########### end ##############