-
Notifications
You must be signed in to change notification settings - Fork 0
/
RinA CH12 Code.txt
executable file
·140 lines (106 loc) · 3.65 KB
/
RinA CH12 Code.txt
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
#------------------------------------------------------------#
# R in Action: Chapter 12 #
# requires that the coin, multcomp, vcd, MASS, #
# lmPerm, and boot packages have been installed #
# install.packages(c('coin','multcomp', 'vcd', 'MASS', #
# 'lmPerm', 'boot')) #
#------------------------------------------------------------#
# Listing 12.1 - 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
library(MASS)
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)))
set.seed(1234)
chisq_test(Treatment ~ Improved, data = Arthritis,
distribution = approximate(B = 9999))
# independence between numeric variables
states <- as.data.frame(state.x77)
set.seed(1234)
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")
# Listing 12.2 - Permutation tests for simple linear regression
library(lmPerm)
set.seed(1234)
fit <- lmp(weight ~ height, data = women, perm = "Prob")
summary(fit)
# Listing 12.3 - Permutation tests for polynomial regression
library(lmPerm)
set.seed(1234)
fit <- lmp(weight ~ height + I(height^2), data = women, perm = "Prob")
summary(fit)
# Listing 12.4 - Permutation tests for multiple regression
library(lmPerm)
set.seed(1234)
states <- as.data.frame(state.x77)
fit <- lmp(Murder ~ Population + Illiteracy +
Income + Frost, data = states, perm = "Prob")
summary(fit)
# Listing 12.5 - Permutation test for One-Way ANOVA
library(lmPerm)
library(multcomp)
set.seed(1234)
fit <- aovp(response ~ trt, data = cholesterol,
perm = "Prob")
summary(fit)
# Listing 12.6 - Permutation test for One-Way ANCOVA
library(lmPerm)
set.seed(1234)
fit <- aovp(weight ~ gesttime + dose, data = litter,
perm = "Prob")
summary(fit)
# Listing 12.7 - Permutation test for Two-way ANOVA
library(lmPerm)
set.seed(1234)
fit <- aovp(len ~ supp * dose, data = ToothGrowth,
perm = "Prob")
summary(fit)
# bootstrapping a single statistic (R2)
# pause on each graph
par(ask = TRUE)
rsq <- function(formula, data, indices) {
d <- data[indices, ]
fit <- lm(formula, data = d)
return(summary(fit)$r.square)
}
library(boot)
set.seed(1234)
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)
bs <- function(formula, data, indices) {
d <- data[indices, ]
fit <- lm(formula, data = d)
return(coef(fit))
}
library(boot)
set.seed(1234)
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)