# S-PLUS code to run exercises in
#        "Bootstrap Methods and Permutation Tests",
# a document which can be used independently or as Chapter 18 of
#        "The Practice of Business Statistics",
# published by W. H. Freeman.
#
#
#
# PREFACE
#
# This chapter uses computationally demanding resampling procedures, for
# which the use of a computer is essential.  We used S-PLUS while
# writing this chapter, and give commands below for performing the
# analyses in S-PLUS.  Our goal is to make it as easy as possible for
# you to focus on the statistical aspects of this chapter.  You may use
# S-PLUS for this chapter even if you use other software for other
# chapters.
#
# To use S-PLUS you need:
# * a copy of S-PLUS:  there is a free version for students (and
#   instructors?), see
#         http://elms03.e-academy.com/splus
# * the S+Resample library, free at
#         http://www.insightful.com/downloads/libraries
#   See also
#         http://www.insightful.com/Hesterberg/bootstrap
# * the data: the easiest way to get this is with the S+PBSdata library, and
# * optionally, other materials, such as the documentation in the PBSdata
#   library.  See the README file in the library to see what is available.
#
#
#
#
# EXAMPLE CODE
# To run any of the following code in the Windows version of S-PLUS,
# select the corresponding lines using the mouse, then hit F10 or
# click the triangle in the S-PLUS menu bar.
#
# There are exercises that ask you do do certain procedures if the
# procedures are reasonable to do.  Here we include the code,
# and leave it to you to decide whether to use it.

# 18.1 This is normally done by hand, but the following commands would work:
boot1 = bootstrap(Exercise18.1$Time, mean, B=20)
stem(boot1$replicates)
boot1  # The standard error is the printed "SE"


# 18.2
boot2 = bootstrap(Exercise18.2$Earnings, mean)
plot(boot2)
qqnorm(boot2)
boot2  # the SE is printed
limits.percentile(boot2)


# 18.3
hist(Exercise18.3$Spending)
qqnorm(Exercise18.3$Spending)
boot3 = bootstrap(Exercise18.3$Spending, mean)
plot(boot3)
qqnorm(boot3)


# 18.4
hist(Exercise18.4$Survival)
qqnorm(Exercise18.4$Survival)
boot4 = bootstrap(Exercise18.4$Survival, mean)
plot(boot4)
qqnorm(boot4)


# 18.5
hist(Exercise18.5$Spending)
qqnorm(Exercise18.5$Spending)
boot5 = bootstrap(Exercise18.5$Spending, mean)
par(mfrow=c(2,2))  # two rows and two columns of plots
plot(boot5)
qqnorm(boot5)
plot(boot3)
qqnorm(boot3)
par(mfrow=c(1,1))  # back to the usual single plot per page
boot5
boot3
stdev(Exercise18.5$Spending)
stdev(Exercise18.3$Spending)


# 18.6
hist(Exercise18.6$Survival)
qqnorm(Exercise18.6$Survival)
boot6 = bootstrap(Exercise18.6$Survival, mean)
par(mfrow=c(2,2))  # two rows and two columns of plots
plot(boot6); qqnorm(boot6)
plot(boot4); qqnorm(boot4)
par(mfrow=c(1,1))
boot6
boot4
stdev(Exercise18.6$Survival)
stdev(Exercise18.4$Survival)


# 18.7
stdev(Exercise18.3$Spending)
stdev(Exercise18.3$Spending) / sqrt(50)
boot3
stdev(boot3$replicates) # another way


# 18.8
boot3
mean(boot3$replicates) - boot3$observed  # another way


# 18.9
boot4
mean(boot4$replicates) - boot4$observed  # another way

# 18.10
limits.t(boot3)
t.test(Exercise18.3$Spending)


# 18.11
mean(Exercise18.3$Spending, trim=.25)
mean(Exercise18.3$Spending)
hist(Exercise18.3$Spending)
qqnorm(Exercise18.3$Spending)
boot11 = bootstrap(Exercise18.3, mean(Spending, trim=.25))
limits.t(boot11)


# 18.12
boot12 = bootstrap(Exercise18.3$Spending, median)
plot(boot12)
qqnorm(boot12)
qqnorm(boot12, pch=".") # this gives smaller points, easier to see duplicates


# 18.13
n1 = length(ILEC$Time)
n2 = length(CLEC$Time)
sqrt( var(ILEC$Time) / n1 + var(CLEC$Time) / n2)


# 18.14
Exercise18.14  # has columns score and group
boot14 = bootstrap2(Exercise18.14$score, mean,
		     treatment = Exercise18.14$group)
boot14
plot(boot14)
qqnorm(boot14)
limits.t(boot14)


# 18.15
Exercise18.15  # has columns status and ratio
boot15 = bootstrap2(Exercise18.15$ratio, mean,
		     treatment = Exercise18.15$status)
plot(boot15)
qqnorm(boot15)
boot15
limits.t(boot15)


# 18.16


# 18.17
boot17 = bootstrap(Exercise18.17$Price, mean)
plot(boot17)
qqnorm(boot17)
boot17


# 18.18
boot18 = bootstrap(Exercise18.18$Price, median)
plot(boot18)
qqnorm(boot18)
boot18


# 18.19
hist(Exercise18.19$x)
qqnorm(Exercise18.19$x)
boot19 = bootstrap(Exercise18.19$x, mean)
boot19
plot(boot19)
qqnorm(boot19)
limits.t(boot19)


# 18.20
hist(Exercise18.20$Salary)
qqnorm(Exercise18.20$Salary)
summary(Exercise18.20$Salary)
boot20 = bootstrap(Exercise18.20, mean(Salary, trim=.25))
plot(boot20)
qqnorm(boot20)
limits.t(boot20)
boot20b = bootstrap(Exercise18.20$Salary, median)
plot(boot20b)
qqnorm(boot20b)
limits.t(boot20b)


# 18.21
stdev(Exercise18.21$Weight)
boot21 = bootstrap(Exercise18.21$Weight, stdev)
boot21
plot(boot21)
qqnorm(boot21)


# 18.22
summary(Exercise18.21$Weight)
quantile(Exercise18.21$Weight)
# simplify things by defining our own function
interQuartile = function(x) { a = quantile(x); a[4] - a[2] }
interQuartile(Exercise18.21$Weight)
boot22 = bootstrap(Exercise18.21$Weight, interQuartile)
boot22
plot(boot22)
qqnorm(boot22)
limits.t(boot22)


# 18.23
Exercise18.23  # columns Minority and White
par(mfrow=c(2,2))  # two rows and two columns of plots
hist(Exercise18.23$Minority)
qqnorm(Exercise18.23$Minority)
hist(Exercise18.23$White)
qqnorm(Exercise18.23$White)
par(mfrow=c(1,1))
t.test(Exercise18.23$Minority, Exercise18.23$White, paired=T)
difference = Exercise18.23$Minority - Exercise18.23$White
t.test(difference)
boot23 = bootstrap(difference, mean)
# another way:
# boot23 = bootstrap(Exercise18.23, mean(Minority) - mean(White))
plot(boot23)
qqnorm(boot23)
limits.t(boot23)


# 18.24
sort(Exercise18.24$NetWorth)
hist(Exercise18.24$NetWorth)
qqnorm(Exercise18.24$NetWorth)
boot24 = bootstrap(Exercise18.24, mean(NetWorth, trim=.25))
plot(boot24)
qqnorm(boot24)
boot24
limits.t(boot24)


# 18.25
boot25 = bootstrap(CLEC$Time, mean)
plot(boot25)
qqnorm(boot25)
boot25
boot25b = bootstrap(ILEC$Time, mean)
boot25b
plot(boot25b)
qqnorm(boot25b)


# 18.26
x10 = rnorm(10, mean = 8.4, sd = 14.7)
boot26a = bootstrap(x10, mean)
plot(boot26a)
qqnorm(boot26a)
boot26a
x40 = rnorm(40, mean = 8.4, sd = 14.7)
boot26b = bootstrap(x40, mean)
plot(boot26b)
qqnorm(boot26b)
boot26b
x160 = rnorm(160, mean = 8.4, sd = 14.7)
boot26c = bootstrap(x160, mean)
plot(boot26c)
qqnorm(boot26c)
boot26c


# 18.27
y10 = sample(ILEC$Time, size=10)
boot27a = bootstrap(y10, mean)
plot(boot27a)
qqnorm(boot27a)
boot27a
y40 = sample(ILEC$Time, size=40)
boot27b = bootstrap(y40, mean)
plot(boot27b)
qqnorm(boot27b)
boot27b
y160 = sample(ILEC$Time, size=160)
boot27c = bootstrap(y160, mean)
plot(boot27c)
qqnorm(boot27c)
boot27c


# 18.28


# 18.29


# 18.30
qqnorm(Exercise18.30$IQ)
stdev(Exercise18.30$IQ)/sqrt(length(Exercise18.30$IQ))
t.test(Exercise18.30$IQ)
boot30 = bootstrap(Exercise18.30$IQ, mean)
plot(boot30)
qqnorm(boot30)
boot30
limits.t(boot30)
limits.percentile(boot30)


# 18.31


# 18.32
Exercise18.32  # columns Wages LOS Size
boot32 = bootstrap(Exercise18.32, cor(Wages, LOS))
boot32
plot(boot32)
qqnorm(boot32)
limits.t(boot32)
limits.percentile(boot32)


# 18.33
# To use the following code, replace MIDDLE, STDEV, LEFT, and RIGHT
# with the appropriate numbers.
MIDDLE = what?
STDEV = what?
LEFT = what?
RIGHT = what?
MIDDLE + qt(c(.025, .975), 49) * STDEV
boot33 = bootstrap(Example18.11$Time, mean, B=1000)
plot(boot33)
segments(x1 = LEFT, x2 = RIGHT, y1 = 0.008)
text(470, 0.008, "bootstrap t")
segments(x1 = LEFT, x2 = RIGHT, y1 = 0.0075)
text(470, 0.0075, "percentile")
segments(x1 = LEFT, x2 = RIGHT, y1 = 0.0070)
text(470, 0.0070, "traditional")
segments(x1 = LEFT, x2 = RIGHT, y1 = 0.0065)
text(470, 0.0065, "BCa")
segments(x1 = LEFT, x2 = RIGHT, y1 = 0.0060)
text(470, 0.0060, "tilting")



# 18.34
# To use the following code, replace LEFT and RIGHT
# with the appropriate numbers.
LEFT = what?
RIGHT = what?
boot34 = bootstrap(Seattle2002, mean(Price, trim = .25))
limits.t(boot34)
limits.percentile(boot34)
limits.bca(boot34)
limits.tilt(boot34)
plot(boot34)
segments(x1 = LEFT, x2 = RIGHT, y1 = .020); text(300, .020, "bootstrap t")
segments(x1 = LEFT, x2 = RIGHT, y1 = .018); text(300, .018, "percentile")
segments(x1 = LEFT, x2 = RIGHT, y1 = .016); text(300, .016, "BCa")
segments(x1 = LEFT, x2 = RIGHT, y1 = .014); text(300, .014, "tilting")



# 18.35
boot35 = boot32  # see above
boot35
plot(boot35)
qqnorm(boot35)
limits.t(boot35)
limits.percentile(boot35)
limits.bca(boot35)
limits.tilt(boot35)


# 18.36
hist(CLEC$Time)
qqnorm(CLEC$Time)
mean(CLEC$Time)
boot36 = bootstrap(CLEC$Time, mean)
plot(boot36)
qqnorm(boot36)
boot36
limits.t(boot36)
limits.percentile(boot36)
limits.bca(boot36)
limits.tilt(boot36)


# 18.37
# boot37 = bootstrap2(Verizon, mean(Time), treatment = Group)
# Previous works, but next version makes BCa calculations faster
boot37 = bootstrap2(Verizon$Time, mean, treatment = Verizon$Group)
limits.bca(boot37)
limits.tilt(boot37)

# 18.38
limits.percentile(boot19)  # see problem19, where boot19 was created
limits.t(boot19)


# 18.39
t.test(Exercise18.21$Weight)
hist(Exercise18.21$Weight)
qqnorm(Exercise18.21$Weight)
boot39 = bootstrap(Exercise18.21$Weight, mean)
limits.percentile(boot39)
plot(boot39)
qqnorm(boot39)


# 18.40
hist(Exercise18.40$Earnings)
qqnorm(Exercise18.40$Earnings)
boot40 = bootstrap(Exercise18.40, mean(Earnings, trim=.25))
boot40
boot40b = bootstrap(Exercise18.40$Earnings, mean)
boot40b
plot(boot40)
qqnorm(boot40)
limits.bca(boot40); limits.tilt(boot40)
limits.t(boot40); limits.percentile(boot40)
limits.bca(boot40b); limits.tilt(boot40b)
limits.t(boot40b); limits.percentile(boot40b)


# 18.41
hist(Exercise18.41$x)
qqnorm(Exercise18.41$x)
t.test(Exercise18.41$x)
boot41 = bootstrap(Exercise18.41$x, mean)
plot(boot41)
qqnorm(boot41)
limits.percentile(boot41)


# 18.42
limits.bca(boot41)
limits.bca(boot41) - boot41$observed
limits.tilt(boot41)
limits.tilt(boot41) - boot41$observed
limits.percentile(boot41) - boot41$observed


# 18.43
Exercise18.43  # variables Price, SquareFootage, Age
hist(Exercise18.43$Price)
qqnorm(Exercise18.43$Price)
boot43 = bootstrap(Exercise18.43, mean(Price, trim=.25))
boot43
plot(boot43)
qqnorm(boot43)
limits.percentile(boot43)
limits.bca(boot43)
limits.tilt(boot43)
limits.t(boot43)


# 18.44
boot44 = bootstrap(Exercise18.44, cor(Price, SquareFootage))
plot(boot44)
qqnorm(boot44)
limits.t(boot44)
limits.percentile(boot44)
limits.bca(boot44)
limits.tilt(boot44)
limits.bca(boot44) - boot44$observed

# Regression problems follow.  I'll give command-line code, but it would
# be more natural to do this using menus, under:
# 	Statistics:Resample:Linear Regression

# 18.45
Exercise18.45 # variables Weight, Mileage
plot(Exercise18.45$Weight, Exercise18.45$Mileage)
cor(Exercise18.45)
boot45a = bootstrap(Exercise18.45, cor(Weight, Mileage))
limits.bca(boot45a)
limits.tilt(boot45a)
lm45 = lm(Mileage ~ Weight, data=Exercise18.45)
lm45
summary(lm45)
-.0082 + qt( c(.025, .975), df=58) * .0007
boot45b = bootstrap(lm45, coef)
limits.percentile(boot45b)
# The book doesn't ask you to do this, but it is interesting to put
# regression lines on the plot.  Not too many or it is a mess.
bootstrap(lm45, abline(lm45, col=2), B=20)
abline(lm45, lwd=2, col=5)  # Add the original line with a wider line, color


# 18.46
Exercise18.46 # variables Salary, Average
plot(Exercise18.46$Salary/10^6, Exercise18.46$Average)
lm46 = lm(Average ~ I(Salary/10^6), data=Exercise18.46)
lm46
abline(lm46)
boot46 = bootstrap(lm46, coef)
plot(boot46)
qqnorm(boot46)
limits.t(boot46)
limits.percentile(boot46)
limits.bca(boot46)
limits.tilt(boot46)


# 18.47
hist(Exercise18.47$Earnings)
qqnorm(Exercise18.47$Earnings)
stem(Exercise18.47$Earnings)
which( Exercise18.47$Earnings < 14000 ) # observation 10
Exercise18.47[ -10,  ] # omit observation 10
boot47a = bootstrap(Exercise18.47$Earnings, mean)
boot47b = bootstrap(Exercise18.47$Earnings[ -10 ], mean)
par(mfrow=c(2,2))
plot(boot47a)
plot(boot47b)
qqnorm(boot47a)
qqnorm(boot47b)
par(mfrow=c(1,1))
boot47a
boot47b
limits.bca(boot47a)
limits.bca(boot47b)


# 18.48  You should do this by hand, but can use the following as a check:
mean( c(24, 61)) - mean( c(42, 33, 46, 37))
perm48 = permutationTestMeans(c(24,61), data2=c(42,33,46,37),
			       alternative="greater")
perm48


# 18.49
Exercise18.49 # variables Price and Year
t.test(Exercise18.49$Price[Exercise18.49$Year == 2002],
       y=Exercise18.49$Price[Exercise18.49$Year == 2001], var.equal=F)
perm49 = permutationTestMeans(Exercise18.49, treatment = Year)
perm49
plot(perm49)
# boot49 = bootstrap2(Exercise18.49, mean(Price), treatment=Year)
boot49 = bootstrap2(Exercise18.49$Price, mean, treatment=Exercise18.49$Year)
# Either of previous two works, but second is quicker.
plot(boot49)
limits.bca(boot49)
limits.tilt(boot49)


# 18.50
Exercise18.50  # variables status and ratio
t.test(Exercise18.50$ratio[Exercise18.50$status == "Healthy"],
       Exercise18.50$ratio[Exercise18.50$status == "Failed"],
       alternative="greater", var.equal=F)
perm50 = permutationTestMeans(Exercise18.50, treatment=status,
			       alternative="greater")
perm50
plot(perm50)
# boot50 = bootstrap2(Exercise18.50, mean(ratio), treatment=status)
boot50 = bootstrap2(Exercise18.50$ratio, mean, # faster than previous
		     treatment=Exercise18.50$status)
plot(boot50)
limits.bca(boot50)
limits.tilt(boot50)


# 18.51 To use the following code, replace P and N by appropriate values
P = what?
N = what?
sqrt(P * (1-P) / N)


# 18.52


# 18.53
Exercise18.53 # contingency table (use this from the Proportions/Resample menu)
p1 = 108/142
p2 = 15/28
p = (108+15)/(142+28)
(p1-p2)/sqrt(p1*(1-p1)/142 + p2*(1-p2)/28)
(p1-p2)/sqrt(p*(1-p)/142 + p*(1-p)/28)
1 - pnorm( (p1-p2)/sqrt(p*(1-p)/142 + p*(1-p)/28) ) # for one-sided test
# For a permutation test we need one row for each company.  Use this:
Exercise18.53expanded # variables Exclusive and Success
perm53 = permutationTestMeans(Exercise18.53expanded, treatment=Exclusive,
			       alternative="greater")
perm53
permutationTestMeans(Exercise18.53expanded, treatment=Exclusive,
		     alternative="greater", B=99999) # more accurate
plot(perm53)
# Note - z-test fails to take ties into account.
#boot53 = bootstrap2(Exercise18.53, mean(Success), treatment=Exclusive)
boot53 = bootstrap2(Exercise18.53expanded$Success, mean,
		     treatment=Exercise18.53expanded$Exclusive)
plot(boot53)
limits.bca(boot53)
limits.tilt(boot53)


# 18.54
Exercise18.54 # variables Right and Left
perm54 = permutationTestMeans(Exercise18.54$Right, data2=Exercise18.54$Left,
			       paired=T, alternative="less")
perm54
permutationTestMeans(Exercise18.54$Right, data2=Exercise18.54$Left, B=99999,
			       paired=T, alternative="less")
# more samples for higher accuracy
plot(perm54)


# 18.55
Exercise18.55 # variables Salary and Average
# must permute only one of the variables
perm55 = permutationTest(Exercise18.55, resampleColumns = "Average",
               cor(Salary, Average), alternative="greater")
perm55
plot(perm55)


# 18.56
Exercise18.56 # variables Height and Position
stem(Exercise18.56$Height[ Exercise18.56$Position == "Forward" ],
     scale = -1, nl=2)
stem(Exercise18.56$Height[ Exercise18.56$Position == "Center" ],
     scale = -1, nl=2)
# Put those together back-to-back by hand.
perm56 = permutationTestMeans(Exercise18.56, treatment = Position,
			       alternative="less")
perm56
plot(perm56)


# 18.57
Exercise18.57 # variables Time, Distance, Hand; latter "right" or "left"
perm57a = permutationTest2(Exercise18.57, median(Time),
			     treatment=Hand)
plot(perm57a)
perm57a
perm57b = permutationTest2(Exercise18.57, mean(Time, trim=.25),
			     treatment=Hand)
plot(perm57b)
perm57b


# 18.58
Exercise18.58 # variables Weight and Age, either "<60" or ">=60"
stem(Exercise18.58$Weight[ Exercise18.58$Age == "<60" ])
stem(Exercise18.58$Weight[ Exercise18.58$Age == ">=60" ])
# do back-to-back by hand
perm58 = permutationTestMeans(Exercise18.58, treatment=Age,
             alternative="greater")
plot(perm58)
qqnorm(perm58)
perm58


# 18.59
Exercise18.59
names(Exercise18.59expanded)  # variables Sex and LabelUser
perm59 = permutationTestMeans(Exercise18.59expanded, treatment=Sex)
perm59
plot(perm59)
qqnorm(perm59)


# 18.60
names(Exercise18.60)  # Price and Year
perm60 = permutationTest2(Exercise18.60, median(Price),
			    treatment = Year)
plot(perm60)
perm60


# 18.61
names(Exercise18.61) # variables Price, SquareFootage, Age
perm61 = permutationTest(Exercise18.61, resampleColumns = "SquareFootage",
			 cor(SquareFootage, Age))
plot(perm61)
perm61


# 18.62
Exercise18.62
Exercise18.62expanded # variables Treatment (calcium or placebo) and Success
perm62 = permutationTestMeans(Exercise18.62expanded, treatment=Treatment,
                 alternative = "greater")
plot(perm62)
perm62


# 18.63
boot63 = bootstrap2(Exercise18.62$Success, mean,
		      treatment = Exercise18.62$Treatment)
plot(boot63)  # Unusual plot -- high-low pattern.  
# Histograms are sensitive to just where the bin edges are.
plot(boot63, nclass=10)  # more reasonable
plot(boot63, nclass=40)  # shows gaps
qqnorm(boot63)# note the bunching & gaps, explains the high-low pattern
boot63
limits.t(boot63, z=T)


# 18.64
Exercise18.64 # variables Time and Group (ILEC or CLEC)
histogram(~Time | Group, data=Exercise18.64)
qqmath(~Time | Group, data=Exercise18.64)
t.test(Exercise18.64$Time[Exercise18.64$Group == "ILEC"],
       Exercise18.64$Time[Exercise18.64$Group == "CLEC"],
       alternative="less", var.equal=F)
# It is easy to do perm test for difference in means:
perm64a = permutationTestMeans(Exercise18.64, treatment=Group,
				alternative="less")
perm64a
permutationTestMeans(Exercise18.64, treatment=Group, B=99999, # higher accuracy
				alternative="less")
# Harder to do a perm test where the test statistics is a pooled-variance t;
# need to ensure that only one variable is permuted.
# As a side note, if you use the same random numbers, then a permutation
# test using the difference in means gives exactly the same P-value as one
# based on the pooled-variance t-statistics.  So we rarely bother with the
# latter.
perm64b = permutationTest(Exercise18.64, resampleColumns = "Time",
			    alternative="less",
	    t.test(Time[Group == "ILEC"],
		   Time[Group == "CLEC"],
		   alternative="less", var.equal=T)$statistic)
perm64b
plot(perm64a)
plot(perm64b)
# permutation distributions are very skewed


# 18.65
perm65 = permutationTest2(Exercise18.64, stdev(Time), treatment=Group,
                          ratio = T)
perm65
plot(perm65)


# 18.66
Exercise18.66 # variables status (Healthy or Failed) and ratio
qqmath(~ratio | status, data=Exercise18.66)
perm66 = permutationTest2(Exercise18.66, var(ratio), treatment=status,
                          ratio = T)
perm66


# 18.67
Exercise18.67 # variables Pretest and Posttest
perm67 = permutationTestMeans(Exercise18.67$Posttest,
			       data2 = Exercise18.67$Pretest,
			       paired = T, alternative = "greater")
perm67
plot(perm67)


# 18.68
hist(Exercise18.68$Change)
qqnorm(Exercise18.68$Change)
mean(Exercise18.68$Change)
stdev(Exercise18.68$Change) / sqrt(length(Exercise18.68$Change))
boot68 = bootstrap(Exercise18.68$Change, mean)
boot68


# 18.69
x69 = runif(50)
median(x69)
boot69 = bootstrap(x69, median)
plot(boot69)
qqnorm(boot69)
boot69
limits.t(boot69)
limits.bca(boot69)
limits.tilt(boot69)


# 18.70
Exercise18.70 # variables Age and Sex (male or female)
stem(Exercise18.70$Age[ Exercise18.70$Sex == "male" ],
     scale = -1, nl=2)
stem(Exercise18.70$Age[ Exercise18.70$Sex == "female" ],
     scale = -1, nl=2)
# do back-to-back stemplot by hand
perm70 = permutationTestMeans(Exercise18.70, treatment=Sex)
perm70
permutationTestMeans(Exercise18.70, treatment=Sex, B=99999)
plot(perm70)


# 18.71
Exercise18.71 # variables Year, Overseas, U.S.
boot71 = bootstrap(Exercise18.71, cor(Overseas, U.S.))
plot(boot71)
boot71
limits.bca(boot71)
limits.tilt(boot71)


# 18.72
Exercise18.72 # variables Day, Close, Change
x72 = 100 * Exercise18.72$Change / Exercise18.72$Close  # percent change
# There we divided by the close, because it was easy.  Really should divide
# by the previous day's close:
x72b = 100 * Exercise18.72$Change / (Exercise18.72$Close - Exercise18.72$Change)
stdev(x72) # 2.21
stdev(x72b)# also 2.21
boot72 = bootstrap(x72, stdev)
boot72
plot(boot72)
limits.t(boot72)
limits.bca(boot72)
.Last.value - boot72$observed
limits.tilt(boot72)
.Last.value - boot72$observed


# 18.73
Exercise18.73
par(mfrow=c(2,2))
hist(Table18.5$Price)
qqnorm(Table18.5$Price)
hist(Exercise18.73$Price)
qqnorm(Exercise18.73$Price)
par(mfrow=c(1,1))
perm73 = permutationTestMeans(Table18.5, data2 = Exercise18.73)
perm73


# 18.74
t.test(Exercise18.74$Radon)
boot74 = bootstrap(Exercise18.74$Radon, mean)
limits.tilt(boot74)
.Last.value - boot74$observed
plot(boot74)
qqnorm(boot74)


# 18.75


# 18.76
Exercise18.76 # variables Before and After
perm76 = permutationTestMeans(Exercise18.76$Before,
			       data2 = Exercise18.76$After, paired=T,
			       alternative = "less")
perm76


# 18.77
boot77 = bootstrap(Exercise18.76, mean(After)-mean(Before))
plot(boot77)
qqnorm(boot77)
limits.t(boot77)


# 18.78
Exercise18.78 # variables X2001 and X2002
par(mfrow=c(2,2))
hist(Exercise18.78$X2001, main="2001")
qqnorm(Exercise18.78$X2001, main="2001")
hist(Exercise18.78$X2002, main="2002")
qqnorm(Exercise18.78$X2002, main="2002")
par(mfrow=c(1,1))
t.test(Exercise18.78$X2001, Exercise18.78$X2002, paired=T)
perm78 = permutationTestMeans(Exercise18.78$X2002,
			       data2 = Exercise18.78$X2001, paired=T)
perm78
permutationTestMeans(Exercise18.78$X2002, B = 99999,
		     data2 = Exercise18.78$X2001, paired=T)
plot(perm78)


# 18.79
Exercise18.79
40/50
30/45
Exercise18.79expanded # variables Sex (girls and boys) and Like (0 or 1)
perm79 = permutationTestMeans(Exercise18.79expanded, treatment=Sex)
perm79


# 18.80
Exercise18.80 # variables WordCount and Education (High, Medium)
histogram(~WordCount | Education, data=Exercise18.80)
qqmath(~WordCount | Education, data=Exercise18.80)
groupMeans(Exercise18.80$WordCount, group=Exercise18.80$Education)
x80 = split(Exercise18.80$WordCount, Exercise18.80$Education)
boot80a = bootstrap(x80$High, mean)
boot80b = bootstrap(x80$Medium, mean)
boot80a
boot80b
par(mfrow=c(2,2))
plot(boot80a);qqnorm(boot80a)
plot(boot80b);qqnorm(boot80b)
par(mfrow=c(1,1))
limits.percentile(boot80a)
limits.percentile(boot80b)
limits.tilt(boot80a)
limits.tilt(boot80b)
boot80c = bootstrap2(Exercise18.80$WordCount, mean,
		      treatment = Exercise18.80$Education)
limits.percentile(boot80c)


# 18.81
perm81 = permutationTestMeans(Exercise18.80, treatment = Education,
			      alternative = "greater")
perm81


# 18.82
Exercise18.82 # variables Burglaries, When (Before and After)
groupMeans(Exercise18.82$Burglaries, Exercise18.82$When)
histogram(~Burglaries | When, data=Exercise18.82)
qqmath(~Burglaries | When, data=Exercise18.82)
x82 = split(Exercise18.82$Burglaries, Exercise18.82$When)
t.test(x82$Before, x82$After, alternative="greater")
perm82a = permutationTestMeans(Exercise18.82, treatment = When,
			      alternative = "greater")
perm82a
plot(perm82a)
qqnorm(perm82a)
perm82b = permutationTestMeans(Exercise18.82, treatment = When,
			      alternative = "less")
perm82b


# 18.83
boot83 = bootstrap2(Exercise18.82$Burglaries, mean,
		     treatment=Exercise18.82$When)
plot(boot83)
qqnorm(boot83)
boot83
limits.t(boot83)
limits.t(boot83, df = (nrow(Exercise18.82)-2)) # degrees of freedom if pooling
t.test(x82$Before, x82$After)
limits.percentile(boot83)

