%\VignetteEngine{knitr}
%\VignetteIndexEntry{Chapter 06, Horton et al. using mosaic}
%\VignettePackage{Sleuth2}
\documentclass[11pt]{article}

\usepackage[margin=1in,bottom=.5in,includehead,includefoot]{geometry}
\usepackage{hyperref}
\usepackage{language}
\usepackage{alltt}
\usepackage{fancyhdr}
\pagestyle{fancy}
\fancyhf{}

%% Now begin customising things. See the fancyhdr docs for more info.

\chead{}
\lhead[\sf \thepage]{\sf \leftmark}
\rhead[\sf \leftmark]{\sf \thepage}
\lfoot{}
\cfoot{Statistical Sleuth in R: Chapter 6}
\rfoot{}

\newcounter{myenumi}
\newcommand{\saveenumi}{\setcounter{myenumi}{\value{enumi}}}
\newcommand{\reuseenumi}{\setcounter{enumi}{\value{myenumi}}}

\pagestyle{fancy}

\def\R{{\sf R}}
\def\Rstudio{{\sf RStudio}}
\def\RStudio{{\sf RStudio}}
\def\term#1{\textbf{#1}}
\def\tab#1{{\sf #1}}


\usepackage{relsize}

\newlength{\tempfmlength}
\newsavebox{\fmbox}
\newenvironment{fmpage}[1]
     {
   \medskip
   \setlength{\tempfmlength}{#1}
   \begin{lrbox}{\fmbox}
	   \begin{minipage}{#1}
		 \vspace*{.02\tempfmlength}
		 \hfill
	   \begin{minipage}{.95 \tempfmlength}}
		 {\end{minipage}\hfill
		 \vspace*{.015\tempfmlength}
		 \end{minipage}\end{lrbox}\fbox{\usebox{\fmbox}}
	 \medskip
	 }


\newenvironment{boxedText}[1][.98\textwidth]%
{%
\begin{center}
\begin{fmpage}{#1}
}%
{%
\end{fmpage}
\end{center}
}

\newenvironment{boxedTable}[2][tbp]%
{%
\begin{table}[#1]
  \refstepcounter{table}
  \begin{center}
\begin{fmpage}{.98\textwidth}
  \begin{center}
	\sf \large Box~\expandafter\thetable. #2
\end{center}
\medskip
}%
{%
\end{fmpage}
\end{center}
\end{table}		% need to do something about exercises that follow boxedTable
}


\newcommand{\cran}{\href{http://www.R-project.org/}{CRAN}}

\title{The Statistical Sleuth in R: \\
Chapter 6}

\author{
Ruobing Zhang \and Kate Aloisio \and Nicholas J. Horton\thanks{Department of Mathematics, Amherst College, nhorton@amherst.edu}
} 

\date{\today}

\begin{document}


\maketitle
\tableofcontents

%\parindent=0pt


<<pvalues, echo=FALSE, message=FALSE>>=
print.pval = function(pval) {
  threshold = 0.0001
    return(ifelse(pval < threshold, paste("p<", sprintf("%.4f", threshold), sep=""),
                ifelse(pval > 0.1, paste("p=",round(pval, 2), sep=""),
                       paste("p=", round(pval, 3), sep=""))))
}
@

<<setup0, include=FALSE, cache=FALSE>>=
require(knitr)
opts_chunk$set(
  dev="pdf",
  fig.path="figures/",
        fig.height=3,
        fig.width=4,
        out.width=".47\\textwidth",
        fig.keep="high",
        fig.show="hold",
        fig.align="center",
        prompt=TRUE,  # show the prompts; but perhaps we should not do this 
        comment=NA    # turn off commenting of ouput (but perhaps we should not do this either
  )
@
<<setup,echo=FALSE,message=FALSE>>=
require(Sleuth2)
require(mosaic)
trellis.par.set(theme=col.mosaic())  # get a better color scheme 
set.seed(123)
# this allows for code formatting inline.  Use \Sexpr{'function(x,y)'}, for exmaple.
knit_hooks$set(inline = function(x) {
if (is.numeric(x)) return(knitr:::format_sci(x, 'latex'))
x = as.character(x)
h = knitr:::hilight_source(x, 'latex', list(prompt=FALSE, size='normalsize'))
h = gsub("([_#$%&])", "\\\\\\1", h)
h = gsub('(["\'])', '\\1{}', h)
gsub('^\\\\begin\\{alltt\\}\\s*|\\\\end\\{alltt\\}\\s*$', '', h)
})
showOriginal=FALSE
showNew=TRUE
@ 

\section{Introduction}

This document is intended to help describe how to undertake analyses introduced as examples in the Second Edition of the \emph{Statistical Sleuth} (2002) by Fred Ramsey and Dan Schafer.
More information about the book can be found at \url{http://www.proaxis.com/~panorama/home.htm}.
This
file as well as the associated \pkg{knitr} reproducible analysis source file can be found at
\url{http://www.amherst.edu/~nhorton/sleuth}.


This work leverages initiatives undertaken by Project MOSAIC (\url{http://www.mosaic-web.org}), an NSF-funded effort to improve the teaching of statistics, calculus, science and computing in the undergraduate curriculum. In particular, we utilize the mosaic package, which was written to simplify the use of R for introductory statistics courses. A short summary of the R needed to teach introductory statistics can be found in the mosaic package vignette (\url{http://cran.r-project.org/web/packages/mosaic/vignettes/MinimalR.pdf}). 

To use a package within R, it must be installed (one time), and loaded (each session). The package can be installed using the following command:
<<install_mosaic,eval=FALSE>>=
install.packages('mosaic')               # note the quotation marks
@
Once this is installed, it can be loaded by running the command:
<<load_mosaic,eval=FALSE>>=
require(mosaic)
@
This
needs to be done once per session.

In addition the data files for the \emph{Sleuth} case studies can be accessed by installing the \pkg{Sleuth2} package.
<<install_Sleuth2,eval=FALSE>>=
install.packages('Sleuth2')               # note the quotation marks
@
<<load_Sleuth2,eval=FALSE>>=
require(Sleuth2)
@

We also set some options to improve legibility of graphs and output.
<<eval=TRUE>>=
trellis.par.set(theme=col.mosaic())  # get a better color scheme for lattice
options(digits=3)
@
The specific goal of this document is to demonstrate how to calculate the quantities described in Chapter 6: Linear Combinations and Multiple Comparisons of Means using R.

\section{Discrimination Against the Handicapped}

Do equivalent candidates with the same qualifications but different disabilities get treated differentially?  This is the question addressed in case study 6.1 in the \emph{Sleuth}.  

\subsection{Summary statistics and graphical display}

We begin by reading the data and summarizing the variables.

<<>>=
case0601$Handicap = relevel(case0601$Handicap, ref="Amputee")
summary(case0601)
favstats(Score ~ Handicap, data=case0601)
@

A total of \Sexpr{nrow(case0601)} undergraduate students from a U.S. university were randomly assigned to view the tapes, \Sexpr{nrow(subset(case0601, Handicap=="None"))} to each tape. The five kinds of tapes are: \emph{None}, \emph{Amputee}, \emph{Crutches}, \emph{Hearing} and \emph{Wheelchair}. After reviewing the tape, each subject rated the qualifications of the application on 0-10 scale. Among the five handicap conditions, the \emph{Crutches} group gave the highest mean score, while the \emph{Hearing} group gave the lowest mean score. This is summarized on page 150 and in Display 6.1 of the \emph{Sleuth}. 

<<>>=
with(subset(case0601, Handicap=="None"), stem(Score, scale=2))
with(subset(case0601, Handicap=="Amputee"), stem(Score, scale=2))
with(subset(case0601, Handicap=="Crutches"), stem(Score, scale=1))
with(subset(case0601, Handicap=="Hearing"), stem(Score, scale=2))
with(subset(case0601, Handicap=="Wheelchair"), stem(Score, scale=2))

<<fig.height=6, fig.width=10>>=
bwplot(Handicap ~ Score, data=case0601)
@

<<fig.height=6, fig.width=10>>=
densityplot(~ Score, groups=Handicap, auto.key=TRUE, data=case0601)
@

The stem plots show the applicant qualification scores given by objectives. The boxplots and the density plots show that all the distributions are approximately normally distributed.

\subsection{One-way ANOVA}

First we fit the one way analysis of variance (ANOVA) model, using all of the groups.
This corresponds to the interpretations on page 151.
<<>>=
anova(lm(Score ~ Handicap, data=case0601))
@

The p-value of $p=0.03$
provides some evidence that subjects rate qualifications differently according to handicap status.

By default, the use of the linear model (regression) function displays the pairwise differences between the first group and each of the other groups.  Note that the overall test of the model is the same.

<<>>=
summary(lm(Score ~ Handicap, data=case0601))
@

The reference group here is \emph{Amputee}, followed by \emph{None}, \emph{Crutches}, \emph{Hearing} and \emph{Wheelchair}.

Another way of viewing these results is through a model table, which displays the differences between the grand mean and the group means.

<<>>=
model.tables(aov(Score ~ Handicap, data=case0601))
@

Or by:

<<>>=
mean(Score ~ Handicap, data=case0601)-mean(~ Score, data=case0601)
@

\subsection{Contrasts and linear combination}

The Tukey-Kramer test is a reasonable method for these data.  We can
use this to verify the
calculation on page 151.
<<>>=
tuk = TukeyHSD(aov(lm(Score ~ Handicap, data=case0601)), "Handicap", ordered=TRUE, conf.level=0.95)
tuk
@
There is a plot function for a {\tt TukeyHSD} object, which can be displayed by running:
<<eval=FALSE>>=
plot(tuk)
@

<<echo=FALSE>>=
tuk = TukeyHSD(aov(lm(Score ~ Handicap, data=case0601)), "Handicap", ordered=TRUE, conf.level=0.95)
@

Based on the Tukey-Kramer procedure, the difference is estimated to be \Sexpr{round(tuk$Handicap["Crutches-Hearing", "diff"], 2)} points higher for the \emph{Crutches} tapes, with a 95\% confidence from \Sexpr{round(tuk$Handicap["Crutches-Hearing", "lwr"], 2)} to \Sexpr{round(tuk$Handicap["Crutches-Hearing", "upr"], 2)}.

Next, we calculate the comparison of \emph{Amputee/Hearing} to \emph{Crutches/Wheelchair}.

<<>>=
require(gmodels)
fit.contrast(lm(Score ~ Handicap, data=case0601), "Handicap", c(-1, 0, 1, -1, 1), conf.int=0.95)
@

The $t$-value=\Sexpr{round(fit.contrast(lm(Score ~ Handicap, data=case0601), "Handicap", c(-1, 0, 1, -1, 1), conf.int=0.95)[, "t value"], 2)} supports a difference between the average scores given to the \emph{Wheelchair} and  \emph{Crutches} handicaps and the average scores given to the \emph{Amputee} and \emph{Hearing} handicaps.

To verify the calculations on page 155 we used the following contrast:
<<>>=
fit.contrast(lm(Score ~ Handicap, data=case0601), "Handicap", c(-0.5, 0, 0.5, -0.5, 0.5), conf.int=0.95)
@

Other multiple comparison procedures could also be implemented. The following shows the calculation on page 164.
<<>>=
require(agricolae)
LSD.test(aov(lm(Score ~ Handicap, data=case0601)), "Handicap")   # LSD
HSD.test(aov(lm(Score ~ Handicap, data=case0601)), "Handicap")   # Tukey-Kramer
LSD.test(aov(lm(Score ~ Handicap, data=case0601)), "Handicap", p.adj=c("bonferroni"))  # Bonferroni
scheffe.test(aov(lm(Score ~ Handicap, data=case0601)), "Handicap")    # Scheffe
@

The ``Significant Difference" in each test result is the ``95\% interval half-width" described in the book.

\section{Pre-existing Preference of Fish}

Was Charles Darwin right that sexual selection is driven by females?  This is the question addressed in case study 6.2 in the \emph{Sleuth}.  

\subsection{Summary statistics and graphical display}

We begin by reading the data and summarizing the variables.

<<>>=
summary(case0602)
favstats(Proportion ~ Pair, data=case0602)
@

A total of \Sexpr{nrow(case0602)} female fish were involved in this experiment, which is displayed on page 150.

<<fig.height=6, fig.width=10>>=
bwplot(Proportion ~ Pair, data=case0602)
@

Besides the distribution of pair 5, all distributions of other pairs are approximately normally distributed. 

\subsection{One-way ANOVA}

First we fit the one way analysis of variance (ANOVA) model, using all of the groups:

<<>>=
anova(lm(Proportion ~ Pair, data=case0602))
@

The p-value is 0.56, which doesn't provide much evidence that the mean percentage of time with the yellow-sword male differed from one male pair to another.

By default, the use of the linear model (regression) function displays the pairwise differences between the first group and each of the other groups. Note that the overall test of the model is the same.

<<>>=
summary(lm(Proportion ~ Pair, data=case0602))
@

The reference group here is pair 1, followed by pairs 2-6.
Another way of viewing these results is through a model table, which displays the differences between the grand mean and the group means.

<<>>=
model.tables(aov(Proportion ~ Pair, data=case0602))
@

Or by:

<<>>=
mean(Proportion ~ Pair, data=case0602)-mean(~ Proportion, data=case0602)
@

\subsection{Contrasts and linear combination}

We can calculate the values on page 152 and Display 6.5 on page 158 using contrasts.
<<>>=
require(gmodels)
lc = fit.contrast(lm(Proportion ~ Pair, data=case0602), "Pair", c(5, -3, 1, 3, -9, 3), conf.int=0.95); lc
t=round(lc[, "t value"], 2); t
pt(t, 78, lower.tail=TRUE)
@

The $t$-value is \Sexpr{round(fit.contrast(lm(Proportion ~ Pair, data=case0602), "Pair", c(5, -3, 1, 3, -9, 3), conf.int=0.95)[, "t value"], 2)} and the one-sided $p$-value is \Sexpr{round(pt(t, 78, lower.tail=TRUE), 2)}.

<<>>=
mean(mean(Proportion ~ Pair, data=case0602))
t.test(mean(Proportion ~ Pair, data=case0602))
@

The estimated mean percentage of time spent with the yellow-sword male is \Sexpr{round(mean(mean(Proportion ~ Pair, data=case0602)), 3)*100}\%. The one-sided $p$-value$<0.0001$, and the 95\% confidence interval is (\Sexpr{round(t.test(mean(Proportion ~ Pair, data=case0602))$conf.int[1:1], 3)*100}\%, \Sexpr{round(t.test(mean(Proportion ~ Pair, data=case0602))$conf.int[2:2], 3)*100}\%).


\end{document}
