############# R code for variation on problem 6.38, Ott & Longnecker (5th ed.)
### Comparing two populations via two sample & paired t-tests
###############################################################################
# load package "car" (needed by qq.plot)
library(car)
# Data: lung capacity of 12 mice (x1=no ozone) vs. 12 mice (x2=ozone)
x1<-c(8.7,7.9,8.3,8.4,9.2,9.1,8.2,8.1,8.9,8.2,8.9,7.5)
x2<-c(9.4,9.8,9.9,10.3,8.9,8.8,9.8,8.2,9.4,9.9,12.2,9.3)
# side-by-side boxplots
boxplot(x1,x2,names=c("no-ozone","ozone"),main="Lung Capacity")

# first test for equal variances (parametric & nonparametric)
var.test(x1, x2, ratio = 1, alternative ="two.sided") # parametric
fligner.test(x1, x2)                                  # non-par
# two-sample t-test for Ha: mu_1<mu_2, assuming unequal variances
t.test(x1, x2, alternative = "less", mu = 0, paired = FALSE, var.equal = FALSE)
# corresponding Wilcoxon rank sum test:
wilcox.test(x1, x2, alternative = "less", mu = 0, paired = FALSE)

# paired t-test for Ha: mu_1<mu_2:
t.test(x1, x2, alternative = "less", mu = 0, paired = TRUE)
# corresponding Wilcoxon signed rank test:
wilcox.test(x1, x2, alternative = "less", mu = 0, paired = TRUE)
# check normality via qq-plot
qqnorm(x1-x2); qqline(x1-x2)
# Battery of tests for normality of x1 (beyond shapiro need package "nortest")
shapiro.test(x1)
library(nortest)
cvm.test(x1)
ad.test(x1)
lillie.test(x1)
sf.test(x1)
# power
power.t.test(n=12, delta=1, sd=1, sig.level=0.05, power = NULL, type="paired", alternative="one.sided")
             
### Some quantile plots to check distribution assumptions
# simulate 100 values from N(0,1) & do qq-plot, pts should fall close 
# to a straight line through origin of slope 1:
y <- rnorm(100, mean=0, sd=1)
qqnorm(y); abline(0,1,col=2)

# do qq-plot for x1:
par(mfrow=c(3,4))
qqnorm(x1); qqline(x1)

# simulate 11 data sets of n=12 from N(mean(x1),var(x1)^.5), and
# do qq-plot to see how much variation can be expected when n=12:
for (i in 1:11){
 y <- rnorm(12, mean=mean(x1), sd=var(x1)^.5)
 qqnorm(y); qqline(y)
}

# Beware, other distributions such as a t, can also look normal:
par(mfrow=c(2,3))
for (i in 1:6){
 y <- rt(12, df=3)
 qqnorm(y); qqline(y)
}

# In general, qq.plot allows comparison with any distribution, and in addition
# plots a 95% "confidence envelope", e.g. compare t df=3 with chisquare df=3:
for (i in 1:6){
 y <- rt(100, df=3)
 qq.plot(y, dist="chisq", df=3, col="blue")
}
# e.g. compare t df=3 with t df=9:
for (i in 1:6){
 y <- rt(100, df=3)
 qq.plot(y, dist="t", df=9, col="blue")
}
# finally, use qq.plot to check normality for the rat data
par(mfrow=c(1,1))
qq.plot(x1-x2, col="blue")