# R in Action (2nd ed): Chapter 12
# Resampling statistics and bootstrapping
# requires packages coin, multcomp, vcd, MASS, lmPerm, boot
# install.packages(c("coin","multcomp", "vcd", "MASS", "lmPerm", "boot"))
#------------------------------------------------------------------------

par(ask=TRUE)
set.seed(1234) # make results reproducible

# t-test vs. oneway permutation test for the hypothetical data in table 12.1
library(coin)
score <- c(40, 57, 45, 55, 58, 57, 64, 55, 62, 65)
treatment <- factor(c(rep("A",5), rep("B",5)))
mydata <- data.frame(treatment, score)
t.test(score~treatment, data=mydata, var.equal=TRUE)
oneway_test(score~treatment, data=mydata, distribution="exact")

# Wilcoxon Mann-Whitney U test
UScrime <- transform(UScrime, So = factor(So))
wilcox_test(Prob ~ So, data=UScrime, distribution="exact")

# k sample test
library(multcomp)
oneway_test(response~trt, data=cholesterol, 
            distribution=approximate(B=9999))

# independence in contingency tables
library(coin)
library(vcd)
Arthritis <- transform(Arthritis, 
                       Improved = as.factor(as.numeric(Improved)))
chisq_test(Treatment~Improved, data=Arthritis,
           distribution=approximate(B=9999))


# independence between numeric variables
states <- as.data.frame(state.x77)
spearman_test(Illiteracy ~ Murder, data=states, 
              distribution=approximate(B=9999))


# dependent 2-sample and k-sample tests
library(coin)
library(MASS)
wilcoxsign_test(U1 ~ U2, data=UScrime, distribution="exact")

# Permutation tests for simple linear regression
library(lmPerm)
fit <- lmp(weight ~ height, data=women, perm="Prob")
summary(fit)

# Permutation tests for polynomial regression
library(lmPerm)
fit <- lmp(weight ~ height + I(height^2), data=women, perm="Prob")
summary(fit)

# Permutation tests for multiple regression
library(lmPerm)
states <- as.data.frame(state.x77)
fit <- lmp(Murder ~ Population + Illiteracy+Income+Frost,data=states, perm="Prob")
summary(fit)

# Permutation test for One-Way ANOVA
library(lmPerm)
library(multcomp)
fit <- lmp(response ~ trt, data=cholesterol, perm="Prob")
anova(fit)

#  Permutation test for One-Way ANCOVA
library(lmPerm)
fit <- lmp(weight ~ gesttime + dose, data=litter, perm="Prob")
anova(fit)

# Permutation test for Two-way ANOVA
library(lmPerm)
fit <- lmp(len ~ supp*dose, data=ToothGrowth, perm="Prob")
anova(fit)


# bootstrapping a single statistic (R2)
library(boot)
rsq <- function(formula, data, indices) {
         d <- data[indices,]
         fit <- lm(formula, data=d)
         return(summary(fit)$r.square)
} 

results <- boot(data=mtcars, statistic=rsq, 
                R=1000, formula=mpg~wt+disp)
print(results)
plot(results)
boot.ci(results, type=c("perc", "bca"))

# bootstrapping several statistics (regression coefficients)
library(boot)
bs <- function(formula, data, indices) {                
        d <- data[indices,]                                    
        fit <- lm(formula, data=d)                                                                                                
        return(coef(fit))                                    
}
results <- boot(data=mtcars, statistic=bs,             
                R=1000, formula=mpg~wt+disp) 
print(results)
plot(results, index=2)
boot.ci(results, type="bca", index=2)