%                               Intro for STA305 (Experimental Design)
% Notes and comments at the end

% \documentclass[serif]{beamer} % Serif for Computer Modern math font.
\documentclass[serif, handout]{beamer} % Handout mode to ignore pause statements
\hypersetup{colorlinks,linkcolor=,urlcolor=red}


%  To create handout using article mode: Comment above and uncomment below (2 places)
%\documentclass[12pt]{article}
%\usepackage{beamerarticle}
%\usepackage[colorlinks=true, pdfstartview=FitV, linkcolor=blue, citecolor=blue, urlcolor=red]{hyperref}   
     % For live Web links with href in article mode
%\usepackage{fullpage}
 

\usefonttheme{serif} % Looks like Computer Modern for non-math text -- nice!
\setbeamertemplate{navigation symbols}{} % Supress navigation symbols
\usetheme{Berlin} % Displays sections on top
\usepackage[english]{babel}
% \definecolor{links}{HTML}{2A1B81}
% \definecolor{links}{red}

\setbeamertemplate{footline}[frame number] 

\mode<presentation>
% \mode<handout>{\setbeamercolor{background canvas}{bg=black!5}}



\title{Introduction to Experimental Design\footnote{See last slide for copyright information.}}
\subtitle{STA305 Winter 2014}
\date{} % To suppress date


\begin{document}


\begin{frame}
  \titlepage
\end{frame}

\begin{frame}
\frametitle{Background Reading}
\framesubtitle{Optional} 
  \begin{itemize}
    \item Pages 21-26 in Chapter 1 of \emph{Data analysis with SAS}: The correlation-causation issue.
    \item Chapter 11 of \emph{Data analysis with SAS}: It uses R and is at a lower technical level, but gives the general idea of permutation tests.
    \item \emph{Wikipedia} under \emph{Resampling statistics} for more detail and references on permutation tests.
  \end{itemize}
\end{frame}

\section{Goals}

\begin{frame}
\frametitle{Goal of the course}

We want to know the effect of some \emph{treatment} (or combination of treatments) on some \emph{response}. The method should be
\begin{itemize}
     \item Objective (scientific)
     \item As efficient as possible.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Examples}
\begin{center}
\renewcommand{\arraystretch}{1.5}
\begin{tabular}{ll} \hline
\textbf{Treatment} & \textbf{Response} \\ \hline
Advertising expenditures & Sales \\
Drug & Health \\
Industrial quality control process & Product quality \\
Fertilizer type & Crop yield \\
Sterilization method & Bacterial/viral load
\end{tabular}
\renewcommand{\arraystretch}{1.0}
\end{center}
\end{frame}

\begin{frame}
\frametitle{We want to be objective}
%\framesubtitle{} 
  \begin{itemize}
    \item Outcome will be measured numerically (including categories).
    \item Admit that the data will be somewhat noisy.
    \item Try it more than once (replication).
    \item Use statistical methods to decide if there was any effect, and if so how much.
    \item Combination of statistics and research design.
    \item The data collection should be planned with the statistical analysis in mind!
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Observational versus experimental studies}
\framesubtitle{Definitions from Cox and Reid's (2000) \emph{Theory of the design of experiments}.} 
  \begin{itemize}
    \item \textbf{Observational study}: Allocation of individuals to treatments is not under the control of the investigator.
    \item \textbf{Experimental study}: The system under study (including allocation of individuals to treatments) is mostly under the control of the investigator.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Experimental units}
%\framesubtitle{} 
\begin{itemize}
    \item The smallest subset of experimental material to which a separate treatment might be applied.
    \item Usually people, rats, stores, test tubes, plots of land, grocery stores, etc.
    \item But if one kind of contact lens is put in the left eye and another kind in the right eye, the experimental unit would be the eye.
\end{itemize}
\end{frame}

\section{Omitted variables}

\begin{frame}
\frametitle{Omitted variables}
%\framesubtitle{} 
  \begin{itemize}
    \item In any study, some things that affect the response will not be part of the data set.
    \item At best, they contribute background noise that makes the effect of the treatment harder to see.
    \item At worst, they are also systematically related to the treatments, and can produce misleading results.
    \item Observational studies are particularly subject to this bias.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Example}
%\framesubtitle{} 

\begin{center}
{\Large Hair length at age 25 and length of life}
\end{center}
\end{frame}

\begin{frame}
\frametitle{Regression}
  \begin{itemize}
    \item The usual conditional regression model assumes that any omitted variables are independent of the explanatory variables in the model.
    \item What happens when this is violated?
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Example of omitted variables in regression}

Independently for $i= 1, \ldots, n$,
\begin{displaymath}
    Y_i = \beta_0 + \beta_1 X_{i,1} + \beta_2 X_{i,2} + \epsilon_i,
\end{displaymath}
where $\epsilon_i \sim N(0,\sigma^2)$. 
The mean and covariance matrix of the independent variables are given by
\begin{displaymath}
    E\left( \begin{array}{c} X_{i,1} \\ X_{i,2} \end{array} \right) =
     \left( \begin{array}{c} \mu_1 \\ \mu_2 \end{array} \right)
     \mbox{~~ and ~~}
    cov\left( \begin{array}{c} X_{i,1} \\ X_{i,2} \end{array} \right) = 
 \left( \begin{array}{rr}
\phi_{11} & \phi_{12} \\ 
\phi_{12} & \phi_{22}
\end{array} \right) 
\end{displaymath}

\vspace{10mm}
But $X_{i,2}$ is not in the data set, so we just use $X_{i,1}$.
\end{frame}

\begin{frame}
\frametitle{What happens to $\widehat{\beta}_1$ as $n \rightarrow \infty$}
\framesubtitle{When $X{i,2}$ is ignored} 
\begin{eqnarray*}
    \widehat{\beta}_1 
    &=& \frac{\sum_{i=1}^n(X_{i,1}-\overline{X}_1)(Y_i-\overline{Y})}
                           {\sum_{i=1}^n(X_{i,1}-\overline{X}_1)^2} \\
&&\\
    &=& \frac{\frac{1}{n}\sum_{i=1}^n(X_{i,1}-\overline{X}_1)(Y_i-\overline{Y})}
             {\frac{1}{n}\sum_{i=1}^n(X_{i,1}-\overline{X}_1)^2} \\
&&\\
    & \rightarrow & \frac{Cov(X_{i,1},Y_i)}{Var(X_{i,1})} \\
&&\\
    &=& \beta_1 + \frac{\beta_2\phi_{12}}{\phi_{11}}
\end{eqnarray*}
% Do bias in a fixed x model as HW
\end{frame}

\begin{frame}
\frametitle{Effects of omitted variables}
\framesubtitle{Illustrated by $\widehat{\beta}_1 \rightarrow \beta_1 + \frac{\beta_2\phi_{12}}{\phi_{11}}$} 
An omitted variable that is associated with \emph{both} the explanatory and the response variable is sometimes called a \emph{confounding variable}. It can
  \begin{itemize}
    \item Produce an association between explanatory and response variable even when, conditionally on the omitted variable, they are independent
    \item Mask a real relationship between explanatory and response variable
    \item Reverse a relationship between explanatory and response variable.
  \end{itemize}
\end{frame}


\begin{frame}
\frametitle{Include everything?}
One possible solution is to include \emph{all} relevant variables, and ``control" for them somehow, maybe by regression or subdivision. But,
  \begin{itemize}
    \item You can't know what they all are.
    \item Data set will be huge and expensive to collect.
    \item Most variables will be measured with error anyway.
  \end{itemize}
\end{frame}

\section{Experiments}

\begin{frame}
\frametitle{The solution is an experiment}
\framesubtitle{Excellent for certain problems, impossible for others} 
Use control over the setting to break up the association between the treatment and potential confounding variables.
  \begin{itemize}
    \item Random assignment: Thank you Mr. Fisher.
    \item There are other ways, like alternate assigment to experimental and control treatments.
    \item The details of how it's done are important.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Random assignment of experimental units to treatments}
\framesubtitle{Say, to an experimental group versus control group} 

Use pseudo-random numbers from a computer, a table of random numbers or a physical game of chance to assign experimental units to treatments.

  \begin{itemize}
    \item Make sure there are \emph{no} systematic differences in what happens to the different groups, other than the treatment of interest.
    \item If there is a systematic relationship between treatment and any omitted variable, it's purely due to chance.
    \item If the treatment has no effect, any difference in the response must also be due to chance. 
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Completely randomized design}
%\framesubtitle{} 
  \begin{itemize}
    \item There are $p$ treatments.
    \item Randomly assign experimental units to treatments .
    \item Assign $n_1$ to Treatment $1$, \ldots, $n_p$ to treatment $p$.
    \item So that all assignments are equally likely.
    \item Could be done with a jar of marbles.
    \item Or a random permutation.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Random permutation}
\framesubtitle{An easy way to do it} 
  \begin{itemize}
    \item Number the experimental units $1, \ldots, n$, where $n = \sum_{j=1}^p n_j$.
    \item Using software, generate $n$ (pseudo) random variables from a uniform distribution on $(0,1)$.
    \item Sort the integers $1, \ldots, n$ by the random numbers.
    \item The randomly scrambled numbers $1, \ldots, n$ are a \emph{random permutation}.
    \item Assign the first $n_1$ units to treatment one, the next $n_2$ units to treatment two, and so on.
    \item All such assignments are equally likely. 
  \end{itemize}
\end{frame}

\section{Permutation tests}

\begin{frame}
\frametitle{Basis of the permutation test}
\framesubtitle{For example, compare an experimental group to a control group.} 
  \begin{itemize}
    \item Calculate a test statistic that expresses \emph{difference} in the response variable values for the different treatment groups. For example, $\overline{Y}_1-\overline{Y}_2$.
    \item Under $H_0$, the treatment does nothing.
    \item So if the test statistic is big, it's just by the luck of random assignment.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{How the permutation test goes}
%\framesubtitle{In theory} 
  \begin{itemize}
    \item Consider the response variable values to be fixed constants.
    \item Under $H_0$, all allocations of these values to treatment groups are equally likely.
    \item For each allocation, there is a value of the test statistic.
    \item This induces a probability distribution of the test statistic under $H_0$. 
    \item The $p$-value is the probability of obtaining a test statistic greater than or equal to the one from the actual experiment (perhaps in absolute value). 
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{A small example}
\framesubtitle{To illustrate the idea} 
  \begin{itemize}
    \item Have $n=6$ automobiles to crash test.
    \item Randomly assign 3 to the experimental condition (new air bag system), and 3 to the existing system.
    \item Measure damage to the crash test dummy.
    \item Damage measurements are 1.3 6.0 3.0 for the new system, and 5.6 6.5 7.1 for the existing system. 
    \item $\overline{Y}_1=3.43$, $\overline{Y}_2=6.4$, and  $\overline{Y}_2-\overline{Y}_1 = 2.97$.
    \item Is the new system better? 
    \item We need a one-sided test here.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Get the permutation distribution of the test statistic}
\framesubtitle{} 
  \begin{itemize}
    \item There are $\binom{6}{3}=20$ ways to divide the observations into 2 groups.
    \item All equally likely under $H_0$
    \item List them, and calculate $\overline{Y}_2-\overline{Y}_1$ for each one.
  \end{itemize}
\end{frame}

\begin{frame}[fragile]
\frametitle{$\binom{6}{3}=20$. We don't need all $6!=720$ permutations.}
 \framesubtitle{} 
{\scriptsize
\begin{verbatim}
                                diff
(1.3 6.0 7.1) (3.0 5.6 6.5)  0.2333333
(1.3 3.0 7.1) (5.6 6.0 6.5)  2.2333333
(1.3 5.6 7.1) (3.0 6.0 6.5)  0.5000000
(1.3 6.5 7.1) (3.0 5.6 6.0) -0.1000000
(3.0 6.0 7.1) (1.3 5.6 6.5) -0.9000000
(5.6 6.0 7.1) (1.3 3.0 6.5) -2.6333333
(6.0 6.5 7.1) (1.3 3.0 5.6) -3.2333333
(3.0 5.6 7.1) (1.3 6.0 6.5) -0.6333333
(3.0 6.5 7.1) (1.3 5.6 6.0) -1.2333333
(5.6 6.5 7.1) (1.3 3.0 6.0) -2.9666667
(3.0 5.6 6.5) (1.3 6.0 7.1) -0.2333333
(5.6 6.0 6.5) (1.3 3.0 7.1) -2.2333333
(3.0 6.0 6.5) (1.3 5.6 7.1) -0.5000000
(3.0 5.6 6.0) (1.3 6.5 7.1)  0.1000000
(1.3 5.6 6.5) (3.0 6.0 7.1)  0.9000000
(1.3 3.0 6.5) (5.6 6.0 7.1)  2.6333333
(1.3 3.0 5.6) (6.0 6.5 7.1)  3.2333333
(1.3 6.0 6.5) (3.0 5.6 7.1)  0.6333333
(1.3 5.6 6.0) (3.0 6.5 7.1)  1.2333333
(1.3 3.0 6.0) (5.6 6.5 7.1)  2.9666667
\end{verbatim} 
} % End size
\end{frame}

\begin{frame}
\frametitle{Permutation distribution of $\overline{Y}_2-\overline{Y}_1$}
%\framesubtitle{} 
\begin{center}
\includegraphics[width=3in]{PermutationDist}
\end{center}
\end{frame}

\begin{frame}[fragile]
\frametitle{Observed  $\overline{Y}_2-\overline{Y}_1 = 2.97$}
%\framesubtitle{} 
{\scriptsize
\begin{verbatim}
> sort(diff)
 [1] -3.2333333 -2.9666667 -2.6333333 -2.2333333 -1.2333333 -0.9000000
 [7] -0.6333333 -0.5000000 -0.2333333 -0.1000000  0.1000000  0.2333333
[13]  0.5000000  0.6333333  0.9000000  1.2333333  2.2333333  2.6333333
[19]  2.9666667  3.2333333
\end{verbatim} 
} % End size

\vspace{30mm}

So one-sided $p = 0.10$
\end{frame}

\begin{frame}
\frametitle{This is beautiful}
\framesubtitle{Thank you, Mr. Fisher} 
  \begin{itemize}
    \item The probability theory is elementary.
    \item There is no pretence of random sampling from some population.
    \item It still works if there is random sampling.
    \item All the randomness in the model comes from random assignment.
    \item Distribution-free.
    \item Small samples are okay.
    \item Any test statistic you want is okay under $H_0$; it's up to you.
    \item Some test statistics are better than others under $H_1$; it depends on \emph{how} $H_0$ is wrong.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{The only problem}
%\framesubtitle{} 
  \begin{itemize}
    \item For larger samples, listing the permutations or combinations of the data and calculating the test statistic for each one can be a huge task.
    \item It was impossible before electronic computers, in most cases.
    \item Now, it's possible to \emph{estimate} the $p$-value of a permutation test even when it can't be obtained exactly.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Monte Carlo estimation of the permutation $p$-value}
\framesubtitle{See the \emph{Wikipedia} under \emph{Resampling statistics}.} 
\begin{itemize}
    \item Place the values of the response variable in a random order.
    \item Compute the test statistic for the randomly shuffled data.
    \item We have randomly sampled a value of the test statistic from its permutation distribution. 
    \item Carry out this procedure a large number of times.
    \item By the Law of Large Numbers, the $p$-value is approximated by the
proportion of randomly generated values that exceed or equal the observed value of the test statistic.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Another kind of approximation}
%\framesubtitle{} 

 Fisher himself considered permutation tests to be entirely theoretical. In his classic \emph{Statistical Methods for Research Workers} (1936)  he wrote, after describing the procedure,

\begin{quote}
Actually, the statistician does not carry out this very tedious process but his conclusions have no justification beyond the fact they could have been arrived at by this very elementary method.
\end{quote}
\end{frame}

\begin{frame}
\frametitle{Analysis of variance is an approximation}
%\framesubtitle{} 
  \begin{itemize}
    \item Standard ANOVA methods can be justified as approximate permutation tests.
    \item Use the theory of Sample Surveys.
    \item Random assignment to a treatment group is like sampling without replacement from a relatively small population.
    \item There is a Central Limit Theorem.
    \item Convergence is fast, but it's still a large-sample argument.
    \item See Cox and Reid's (2000) \emph{Theory of the design of experiments} and the references therein.
  \end{itemize}
\end{frame}


\begin{frame}
\frametitle{Vocabulary}
\framesubtitle{You need to know these terms.} 
  \begin{itemize}
    \item Observational study
    \item Experimental study
    \item Experimental unit
    \item Confounding variable
    \item Completely randomized design
    \item Random permutation
    \item Permutation distribution
    \item Permutation test
  \end{itemize}
\end{frame}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\begin{frame}
\frametitle{Copyright Information}

This slide show was prepared by  \href{http://www.utstat.toronto.edu/~brunner}{Jerry Brunner},
Department of Statistics, University of Toronto. It is licensed under a 
\href{http://creativecommons.org/licenses/by-sa/3.0/deed.en_US}
     {Creative Commons Attribution - ShareAlike 3.0 Unported License}. Use any part of it as you like and share the result freely. The \LaTeX~source code is available from the course website:
\href{http://www.utstat.toronto.edu/~brunner/oldclass/305s14} {\footnotesize \texttt{http://www.utstat.toronto.edu/$^\sim$brunner/oldclass/305s14}}

\end{frame}


\end{document}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


R for the example

set.seed(9999)
x = round(rgamma(6,5,1),1); x # Exponential, mean = 5
x = rev(x) # Less damage for the new system.
# Oops, I didn't set the seed before I did that. Enter the data.
x <- c(1.3, 3.0, 6.0, 5.6, 6.5, 7.1)

sumx = sum(x)
m1 = mean(x[1:3]); m1
m2 = mean(x[4:6]); m2
m1-m2
# Just need ybar1 to get ybar2-ybar1
sumx/3 - 2*m1

index = read.table(stdin()) # Crtl-D ends
1 2 6
1 3 6
1 4 6
1 5 6
2 3 6
2 4 6
2 5 6
3 4 6
3 5 6
4 5 6
3 4 5
2 4 5
2 3 5
2 3 4
1 4 5
1 3 5
1 3 4
1 2 5
1 2 4
1 2 3



index = as.matrix(index)

expl = control = NULL; diff = numeric(20)
for(j in 1:20)
    {
    expl = rbind(expl,sort(x[index[j,]]))
    control = rbind(control,sort(x[-index[j,]]))
    diff[j] = mean(control[j,])-mean(expl[j,])
    }

 cbind(expl,control,diff)

                                diff
(1.3 6.0 7.1) (3.0 5.6 6.5)  0.2333333
(1.3 3.0 7.1) (5.6 6.0 6.5)  2.2333333
(1.3 5.6 7.1) (3.0 6.0 6.5)  0.5000000
(1.3 6.5 7.1) (3.0 5.6 6.0) -0.1000000
(3.0 6.0 7.1) (1.3 5.6 6.5) -0.9000000
(5.6 6.0 7.1) (1.3 3.0 6.5) -2.6333333
(6.0 6.5 7.1) (1.3 3.0 5.6) -3.2333333
(3.0 5.6 7.1) (1.3 6.0 6.5) -0.6333333
(3.0 6.5 7.1) (1.3 5.6 6.0) -1.2333333
(5.6 6.5 7.1) (1.3 3.0 6.0) -2.9666667
(3.0 5.6 6.5) (1.3 6.0 7.1) -0.2333333
(5.6 6.0 6.5) (1.3 3.0 7.1) -2.2333333
(3.0 6.0 6.5) (1.3 5.6 7.1) -0.5000000
(3.0 5.6 6.0) (1.3 6.5 7.1)  0.1000000
(1.3 5.6 6.5) (3.0 6.0 7.1)  0.9000000
(1.3 3.0 6.5) (5.6 6.0 7.1)  2.6333333
(1.3 3.0 5.6) (6.0 6.5 7.1)  3.2333333
(1.3 6.0 6.5) (3.0 5.6 7.1)  0.6333333
(1.3 5.6 6.0) (3.0 6.5 7.1)  1.2333333
(1.3 3.0 6.0) (5.6 6.5 7.1)  2.9666667

hist(diff,main="Permutation distribution",xlab="Difference between sample means")

> sort(diff)
 [1] -3.2333333 -2.9666667 -2.6333333 -2.2333333 -1.2333333 -0.9000000 -0.6333333 -0.5000000
 [9] -0.2333333 -0.1000000  0.1000000  0.2333333  0.5000000  0.6333333  0.9000000  1.2333333
[17]  2.2333333  2.6333333  2.9666667  3.2333333