# Here is R code for some simple analysis of the 
# mouse data given in Ch. 2 of E&T.

Treatment <- c(94, 197, 16, 38, 99, 141, 23)
Control <- c(52, 104, 146, 10, 51, 30, 40, 27, 46)

# First I'll create a simple dot plot of the data.
group <- c( rep(1, length(Control)), rep(2, length(Treatment)) )
data <- c( Control, Treatment )
stripchart( data ~ group )

# Next I'll create an empirical Q-Q plot.
qqplot(Treatment, Control)
abline(0, 1)

# It seems safe to assume that one distribution is
# stochastically larger than the other one if they 
# differ.  So I can use the Wilcoxon rank sum test
# to do a test about the means.  (R gives an exact
# p-value in this case, because the sample sizes 
# are small and there are no ties.)
wilcox.test(Treatment, Control, alternative="greater")

# Alternatively, one might choose to rely on the 
# robustness of Welch's test.  (R's t.test function,
# when applied to two samples, performs Welch's test
# unless one specifies that the variances should be 
# assumed to be equal.)
t.test(Treatment, Control, alternative="greater")

# Note that even though the dot plot may have suggested   
# that the treament is effective, with such a small amount
# of data we aren't able to claim that there is 
# statistically significant evidence that this is the case.

# Since we're not able to claim that the treatment is
# effective using either test it isn't as important
# to check on the approximate normality of the underlying
# distributions --- but I'll do it anyway.  One way is to
# use R's "canned" probit plot function.
qqnorm(Control)

# Unfortunately (for me) R's canned function doesn't make
# the plot the way I want it, so I'll write my own function
# that will make a plot with the axes reversed.  
cdsqqnorm <- function(x) {
  ox <- sort(x)
  n <- length(x)
  i <- c(1:n)
  q <- (i - 3/8)/(n + 1/4)
  aesnos <- qnorm(q)
  return( plot( ox, aesnos, ylab=c("normal scores"), xlab=c("observed data") ) )
 }

# Now I will use my function to make the desired plots --- one
# for each sample.
cdsqqnorm(Control)

cdsqqnorm(Treatment)

# *** Note: It was pointed out to me during class that there was 
#     no need for me to write my own function to create a probit 
#     plot the way I like it --- I can get what I want using R's
#     qqnorm function by overriding the default value of the 
#     argument datax.  I'll demonstrate below.
qqnorm(Treatment, datax=TRUE)