%                               Power 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}
\usepackage[makeroom]{cancel}
% \definecolor{links}{HTML}{2A1B81}
% \definecolor{links}{red}

\setbeamertemplate{footline}[frame number] 

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



\title{Selection of Sample Size\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 \emph{Data analysis with SAS}, Chapter 8.
    \item For technical details about the non-central distributions, Rencher and Schaalje's \emph{Linear models in statistics},  Section 5.4
  \end{itemize}
\end{frame}

\section{Introduction}

\begin{frame}
\frametitle{How many subjects?}
\begin{itemize}
     \item In planning an experiment, one of the biggest decisions is how many experimental units you need.
     \item[]
     \item Need for what?
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{We never know the truth}
%\framesubtitle{} 
\begin{itemize}
    \item We wish we knew the exact values of all the parameters, but we never do.
    \item All we can hope for is to make good decisions (guesses) with high probability.  
    \item We might guess the exact value of a parameter: Estimation.
    \item We might guess whether some statement about the parameter is right: Testing.
\end{itemize}
\end{frame}

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

\section{Estimation}

\begin{frame}
\frametitle{Estimate the value of a linear combination}
%\framesubtitle{}

\begin{eqnarray*}
    \ell       & = & a_1\mu_1 + a_2\mu_2 + \cdots + a_p\mu_p \\
&&\\
    \hat{\ell} & = & a_1\overline{Y}_1 + a_2\overline{Y}_2 + \cdots + a_p\overline{Y}_p
\end{eqnarray*}

\begin{itemize}
    \item The estimate will be wrong, with probability one.
    \item What is the probability that the estimate will be wrong by less than some margin of error $m$?
    \item Choose sample size to make this probability acceptably high. 
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Expected value and variance}
\framesubtitle{\begin{eqnarray*}
    \ell       & = & a_1\mu_1 + a_2\mu_2 + \cdots + a_p\mu_p \\
    \hat{\ell} & = & a_1\overline{Y}_1 + a_2\overline{Y}_2 + \cdots + a_p\overline{Y}_p
\end{eqnarray*}     } 

\begin{eqnarray*}
    E(\hat{\ell})    & = & \ell\\
    Var(\hat{\ell})  & = & \frac{\sigma^2}{n}\sum_{j=1}^p\frac{a_j^2}{f_j} 
\end{eqnarray*}
\begin{itemize}
    \item $f_j = \frac{n_j}{n}$ are relative sample sizes that add to one.
    \item $f_j$ can be chosen to minimize the variance for any $\sigma^2$ and $n$.
    \item Assume they have been chosen in a smart way, and treat them as fixed.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Probability that $\hat{\ell}$ is close to $\ell$}
%\framesubtitle{}
First, what is the distribution of $\hat{\ell}$?


\begin{eqnarray*}
    Pr\{|\hat{\ell}- \ell| < m \} 
    & = & Pr\left\{\left|\frac{\hat{\ell}- \ell}
             {\sqrt{\frac{\sigma^2}{n}\sum_{j=1}^p\frac{a_j^2}{f_j}} } \right|
          \leq \frac{m}{\sqrt{\frac{\sigma^2}{n}\sum_{j=1}^p\frac{a_j^2}{f_j}}}\right\}   \\
    & = &  Pr\left\{ |Z| < \frac{\sqrt{n}m}
                           {\sigma \sqrt{ \sum_{j=1}^p\frac{a_j^2}{f_j}} } \right\}
\end{eqnarray*}
\end{frame}

\begin{frame}
\frametitle{Make the probability as large as you like}
%\framesubtitle{} 
\begin{displaymath}
    Pr\left\{ |Z| < \frac{\sqrt{n}m}
                         {\sigma \sqrt{\sum_{j=1}^p\frac{a_j^2}{f_j}} } \right\}
\end{displaymath}
\begin{center}
\includegraphics[width=4in]{NormalCurve}
\end{center}
The probability increases to one as $n \rightarrow \infty$.
\end{frame}

\begin{frame}
\frametitle{The probability $Pr\left\{ |Z| < \frac{\sqrt{n}m}
                         {\sigma \sqrt{\sum_{j=1}^p\frac{a_j^2}{f_j}} } \right\}$}
%\framesubtitle{} 
\begin{center}
\includegraphics[width=4in]{NormalCurve}
\end{center}

\begin{itemize}
    \item The probability increases to one as $n \rightarrow \infty$.
    \item Set right hand side to 1.96 to and solve for $n$ get a probability of 0.95, etc.
    \item The $a_j$ and $f_j$ are known.
    \item The value of $m$ depends on what you mean by ``close." 
    \item But $\sigma$ is tougher.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Do you really have to know $\sigma$?}
%\framesubtitle{} 
\begin{itemize}
    \item Maybe you have a good idea of $\sigma$ from other, similar studies. This is most likely if you are planning a follow-up study. Use $\sqrt{MSE}$. 
    \item Or maybe you can give a high guess and a low guess of $\sigma$, and get a \emph{range} of sample sizes you might need.
    \item This is a lot of guessing.
    \item There is another possibility.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Another way out}
%\framesubtitle{} 
Try to express the desired margin of error in \emph{units of the common standard deviation} $\sigma$. 
\begin{itemize}
    \item ``Say something like ``I want my estimate to be within one-tenth of a standard deviation of the correct answer, with probability $0.90$.
    \item To get an idea of what this means, the mean heights of Canadian men and women are about two standard deviations apart.
    \item A poll estimates that 40\% intend to vote for a candidate, and says ``These results are expected to be accurate within four percentage points, 19 times out of 20." They are estimating a margin of error of around 0.08 of a standard deviation, with probability 0.95. % 0.04/sqrt(.4*(1-.4)) 
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Express the desired margin of error in units of $\sigma$}
\framesubtitle{Replace $m$ with $m\sigma$} 
\begin{displaymath}
    Pr\left\{ |Z| < \frac{\sqrt{n}m\bcancel{\sigma}}
                         {\bcancel{\sigma} \sqrt{\sum_{j=1}^p\frac{a_j^2}{f_j}} } \right\}
= Pr\left\{ |Z| < \frac{\sqrt{n}m}
                         { \sqrt{\sum_{j=1}^p\frac{a_j^2}{f_j}} } \right\}
\end{displaymath}
% bcancel for back cancel. See also cancel and xcancel. Need \usepackage[makeroom]{cancel}

\vspace{3mm}
To get probability $1-\alpha$, set 
\begin{displaymath}
    z_{\alpha/2} = \frac{\sqrt{n}m}
                         { \sqrt{\sum_{j=1}^p\frac{a_j^2}{f_j}} } 
\end{displaymath}
\begin{center}
\includegraphics[width=4in]{NormalCurve}
\end{center}

\end{frame}


\begin{frame}
\frametitle{Solve for $n$}
%\framesubtitle{} 
\begin{displaymath}
z_{\alpha/2} = \frac{\sqrt{n}m}
                         { \sqrt{\sum_{j=1}^p\frac{a_j^2}{f_j}} }
~~~\Leftrightarrow~~~
n =  \frac{z_{\alpha/2}^2  \sum_{j=1}^p\frac{a_j^2}{f_j} }{m^2}
\end{displaymath}

\vspace{5mm}

And take the next higher integer. If you ``know" $\sigma$, let
\begin{displaymath}
n =  \frac{\sigma^2 \, z_{\alpha/2}^2 \sum_{j=1}^p\frac{a_j^2}{f_j} }{m^2}
\end{displaymath}

\end{frame}



\begin{frame}
\frametitle{Estimate $\mu_1-\mu_2$}
%\framesubtitle{} 
We want the estimate to be accurate to within $\frac{\sigma}{10}$, with probability 0.95.
\begin{itemize}
    \item $m=\frac{1}{10}$ standard deviations.
    \item $a_1=1$, $a_2=-1$
    \item Equal sample sizes, so $f_1=f_2=\frac{1}{2}$
    \item $z_{\alpha/2}=1.96$
\end{itemize}

\begin{eqnarray*}
    n & = & \frac{z_{\alpha/2}^2  \sum_{j=1}^p\frac{a_j^2}{f_j} }{m^2}\\
      & = &  \frac{1.96^2 (\frac{1^2}{1/2}+\frac{(-1)^2}{1/2})}{(1/10)^2}\\
      & = & 1536.64
\end{eqnarray*} 
So need $n=1537$, or $n_1=n_2=1538/2 = 769$ per group. 
\end{frame}

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

\section{Testing (Power)}

\begin{frame}
\frametitle{Testing (null) hypotheses}
%\framesubtitle{} 
\begin{itemize}
    \item Goal is to make correct decisions with high probability.
    \item When $H_0$ is true, probability of a correct decision (don't reject) is $1-\alpha$. That's guaranteed if the model is correct.
    \item When $H_0$ if false, we want to reject it with high probability.
    \item The probability of rejecting the null hypothesis when the null hypothesis is false is called the \emph{power} of the test.
    \item Power is one minus the probability of a Type II error.
    \item It is a function of the true parameter values.
    \item And also the design, including total sample size. 
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Power is an increasing function of sample size}
%\framesubtitle{} 
\begin{itemize}
    \item Usually, when $H_0$ is false, larger sample size yields larger power.
    \item If power goes to one as a limit when $H_0$ is false (regardless of the exact parameter values) the test is called \emph{consistent}.
    \item Most commonly used tests are consistent, including the general linear $F$-test. 
    \item This means that if $H_0$ is false, you can make the power as high as you wish by making the sample size bigger.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Strategy}
%\framesubtitle{} 
\begin{itemize}
    \item Pick an effect you'd like to be able to detect. An ``effect" means a way that $H_0$ is wrong. It should be just over the boundary of interesting and meaningful.
    \item Pick a desired power -- a probability with which you'd like to be able to detect the effect by rejecting the null hypothesis.
    \item Start with a fairly small $n$ and calculate the power. Increase the sample size until the desired power is reached.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Distribution theory}
%\framesubtitle{} 
\begin{itemize}
    \item We need to study the distribution of the test statistic when the null hypothesis is \emph{false}.
    \item All the distributions you've seen ($Z$, $t$, $\chi^2$, $F$) were derived under the assumption that $H_0$ is \emph{true}.
    \item Here we go.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Intermediate Goal}
%\framesubtitle{} 
\begin{itemize}
    \item For the regression model $\mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$
    \item Testing $H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{t}$
    \item We need to know the distribution of the $F$ statistic when $H_0$ is false.
    \item It will be called the ``non-central" $F$ distribution.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Non-central chi-square with $df=1$}
%\framesubtitle{} 
Let $Z \sim N(\mu,1)$. Then $W = (Z-\mu)^2 \sim \chi^2(1)$, and $M_W(t) = (1-2t)^{-\frac{1}{2}}$. 
If $Z$ is not ``centered" by subtracting off $\mu$,
\begin{eqnarray*}
    M_{Z^2}(t) & = & E(e^{Z^2t}) \\
               & = & \frac{1}{\sqrt{2\pi}} \int_{-\infty}^\infty e^{z^2t} \, e^{-\frac{1}{2}(z-\mu)^2} \, dz \\
               & = & (1-2t)^{-\frac{1}{2}} e^{\mu^2t/(1-2t)}
\end{eqnarray*} 
for $t<1/2$.
\end{frame}

\begin{frame}
\frametitle{Definition of the non-central chi-squared distribution}
\framesubtitle{Generalizing $M_W(t) = (1-2t)^{-\frac{1}{2}} e^{\mu^2t/(1-2t)}$} 
The positive random variable $W$ is said to have a non-central chi-squared distribution with degrees of freedom $\nu>0$ and non-centrality parameter $\lambda \geq 0$ if 

\begin{displaymath}
    M_W(t) = (1-2t)^{-\frac{\nu}{2}} e^{\frac{\lambda t}{1-2t}}
\end{displaymath}
for $t<1/2$.

\begin{itemize}
    \item If $\lambda=0$, this reduces to the ordinary central chi-squared.
    \item We have seen that if $Z \sim N(\mu,1)$, then $W = Z^2 \sim \chi^2(\nu=1,\lambda=\mu^2)$.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Sum of independent chi-squares}
\framesubtitle{Use $M_W(t) = (1-2t)^{-\frac{\nu}{2}} e^{\frac{\lambda t}{1-2t}}$} 
Let $Z_1, \ldots, Z_p \stackrel{ind}{\sim} N(\mu_j,1)$, and $W = \sum_{j=1}^p Z^2_j$. Then 
\begin{eqnarray*}
    M_W(t) & = & \prod_{j=1}^p M_{Z^2_j}(t) \\
           & = & \prod_{j=1}^p (1-2t)^{-\frac{1}{2}} e^{\mu_j^2t/(1-2t)} \\
           & = & (1-2t)^{-\frac{p}{2}} e^{\frac{\left(\sum_{j=1}^p \mu^2_j\right) t}{1-2t}}
\end{eqnarray*}
So $W = \mathbf{Z}^\prime\mathbf{Z} \sim \chi^2\left(p,\lambda=\sum_{j=1}^p \mu^2_j\right)$.
\end{frame}

\begin{frame}
\frametitle{Non-central $F$}
%\framesubtitle{} 

Let $W_1 \sim \chi^2(\nu_1,\lambda)$ and $W_2 \sim \chi^2(\nu_2)$ be independent. Then

\begin{displaymath}
    F^* = \frac{Y_1/\nu_1}{Y_2/\nu_2} \sim F(\nu_1,\nu_2,\lambda) 
\end{displaymath} 
is said to have a \emph{non-central} $F$ distribution with degrees of freedom $\nu_1$ and $\nu_2$, and non-centrality parameter $\lambda$. Write $F^* \sim F(\nu_1,\nu_2,\lambda)$.

\begin{itemize}
    \item Reduces to the ordinary central $F$ when $\lambda=0$.
    \item Good numerical numerical methods are available for calculating the probabilities.
\end{itemize}

\end{frame}

\begin{frame}
\frametitle{Theorem}
\framesubtitle{Proof omitted} 
If $F^* \sim F(\nu_1,\nu_2,\lambda)$, then

\begin{itemize}
    \item $F^*$ is stochastically increasing in $\lambda$, meaning that for every $x>0$, $Pr\{F^*>x|\lambda \}$ is an increasing function of $\lambda$.
    \item That is, the bigger the non-centrality parameter, the greater the probability of getting $F^*$ above any point (such as a critical value).
    \item $lim_{\lambda \rightarrow \infty} Pr\{F^*>x|\lambda \} = 1$. 
\end{itemize}
\end{frame}


\begin{frame}
\frametitle{Heading for the general linear test of $H_0: \mathbf{C}\boldsymbol{\beta}=\mathbf{t}$}
\framesubtitle{With $F = \frac{(\mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t})}
           {q \, MSE}$ } 
\begin{itemize}
    \item Let $\mathbf{Y} \sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})$.
    \item Recall $(\mathbf{Y}-\boldsymbol{\mu})^\prime
           \boldsymbol{\Sigma}^{-1}(\mathbf{Y}-\boldsymbol{\mu}) \sim \chi^2(p)$
    \item When $\mathbf{Y}$ is not centered, 
    \end{itemize}
    
\vspace{4mm}    
\begin{displaymath}
    \mathbf{Y}^\prime
           \boldsymbol{\Sigma}^{-1}\mathbf{Y} \sim \chi^2(p, 
           \boldsymbol{\mu}^\prime \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}).
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{Let $\mathbf{Y} \sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})$. Then
$\mathbf{Y}^\prime\boldsymbol{\Sigma}^{-1}\mathbf{Y} \sim \chi^2(p, 
           \boldsymbol{\mu}^\prime \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu})$.}
\framesubtitle{Proof}
{\small
\begin{itemize}
    \item Because $\boldsymbol{\Sigma}$ is symmetric and positive definite, the symmetric matrix $\boldsymbol{\Sigma}^{-1/2}$ exists.
    \item Let $\mathbf{Z} = \boldsymbol{\Sigma}^{-1/2}\mathbf{Y}$. Note $\mathbf{Z}^\prime\mathbf{Z} = \mathbf{Y}^\prime\boldsymbol{\Sigma}^{-1}\mathbf{Y}$.
    \item $\mathbf{Z}$ is multivariate normal with  $E(\mathbf{Z}) = \boldsymbol{\Sigma}^{-1/2}\boldsymbol{\mu}$ and
\begin{eqnarray*}
    cov(\mathbf{Z}) 
& = & \boldsymbol{\Sigma}^{-1/2}cov(\mathbf{Y})(\boldsymbol{\Sigma}^{-1/2})^\prime \\
& = & \boldsymbol{\Sigma}^{-1/2}\boldsymbol{\Sigma\Sigma}^{-1/2} \\
& = & \boldsymbol{\Sigma}^{-1/2} \boldsymbol{\Sigma}^{1/2} 
      \boldsymbol{\Sigma}^{1/2}\boldsymbol{\Sigma}^{-1/2} \\
& = & \mathbf{I\cdot I} =  \mathbf{I}
\end{eqnarray*} 

    \item Thus the elements of $\mathbf{Z}$ are independent normal with variance one, and $\sum_{j=1}^p Z^2_j = \mathbf{Z}^\prime \mathbf{Z}$ is non-central chi-squared with $df=p$ and non-centrality parameter $\lambda = \sum_{j=1}^p E(Z_j)^2 =  E(\mathbf{Z})^\prime E(\mathbf{Z}) = \left(\boldsymbol{\Sigma}^{-1/2}\boldsymbol{\mu}\right)^\prime \boldsymbol{\Sigma}^{-1/2}\boldsymbol{\mu} = \boldsymbol{\mu}^\prime \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}$.  $\blacksquare$
\end{itemize}
} % End size
\end{frame}

\begin{frame}
\frametitle{General linear test}
% \framesubtitle{When $H_0$ is false} 
\begin{itemize}
    \item Assume the linear model $\mathbf{Y} = \mathbf{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}$
    \item $H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{t}$ may be false.
\end{itemize}

\begin{displaymath}
    F^* = \frac{(\mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t})/q}
           {MSE} \sim F\left(q,n-p,\lambda \right)
\end{displaymath} 
where
\begin{displaymath}
    \lambda = \frac{(\mathbf{C}\boldsymbol{\beta}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{t})}
           {\sigma^2} ~~~~~~~~~~~~~~~~~~~~~
\end{displaymath} 
\end{frame}

\begin{frame}
\frametitle{To show it}
%\framesubtitle{} 
Recall that if $\mathbf{Y} \sim N_p(\boldsymbol{\mu}, \boldsymbol{\Sigma})$, then
$\mathbf{Y}^\prime\boldsymbol{\Sigma}^{-1}\mathbf{Y} \sim \chi^2(p, 
           \boldsymbol{\mu}^\prime \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu})$.

\vspace{3mm}

\begin{itemize}
    \item $\widehat{\boldsymbol{\beta}} \sim N_p\left(\boldsymbol{\beta}, 
                   \sigma^2 (\mathbf{X}^\prime \mathbf{X})^{-1}\right)$
    \item $\mathbf{Y} = \mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t}
           \sim N_q(\mathbf{C}\boldsymbol{\beta}-\mathbf{t}, 
                    \sigma^2 \mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)$
    \item So with $\boldsymbol{\mu}=\mathbf{C}\boldsymbol{\beta}-\mathbf{t}$ and $\boldsymbol{\Sigma} = \sigma^2 \mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime$, the numerator is non-central chi-squared divided by $df$.
    
\begin{displaymath}
    F^* = \frac{(\mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t})^\prime
            (\mathbf{C}(\sigma^2\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\widehat{\boldsymbol{\beta}}-\mathbf{t})/q}
           {MSE/\sigma^2}
\end{displaymath} 
    \item $SSE/\sigma^2 \sim \chi^2(n-p)$ regardless of $H_0$.
    \item And numerator and denominator are independent because $\widehat{\boldsymbol{\beta}}$ and $SSE$ are independent. 
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{The greater the non-centrality parameter $\lambda$, the greater the power}
\framesubtitle{If $H_0$ is false} 
\begin{center}
\includegraphics[width=3.5in]{Fpower}
\end{center}
\end{frame}

 \begin{frame}
\frametitle{What makes $\lambda$ big?}
%\framesubtitle{} 
\begin{displaymath}
    \lambda = \frac{(\mathbf{C}\boldsymbol{\beta}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{t})}
           {\sigma^2}
\end{displaymath} 

\vspace{4mm}

\begin{itemize}
    \item Small $\sigma^2$.
    \item Null hypothesis very wrong.
    \item {\footnotesize Relative sample sizes.}
    \item Total sample size big.
    \item But sample size is hidden in the $\mathbf{X}^\prime \mathbf{X}$ matrix.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{With cell means coding}
%\framesubtitle{} 
\begin{itemize}
    \item Assume there are $p$ treatment combinations.
    \item The $\mathbf{X}$ matrix has exactly one 1 in each row, and all the rest zeros.
    \item There are $n_j$ ones in each column.
\end{itemize}

\vspace{3mm}

\begin{displaymath}
    \mathbf{X}^\prime \mathbf{X} = 
     \left[ \begin{array}{c c c c c}
                 n_1    & 0      &   \cdots  &  0         \\
                 0      & n_2    &   \cdots  &  0         \\
                 \vdots & \vdots &   \ddots  &  \vdots    \\
                 0      & 0      &   \cdots  &  n_p       \\
    \end{array} \right]
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{Multiplying and dividing by $n$}
\framesubtitle{$\lambda = \frac{(\mathbf{C}\boldsymbol{\beta}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{t})}
           {\sigma^2}$} 

\begin{displaymath}
    \lambda = n  \times (\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma})^\prime
            (\mathbf{C}
                 \left[ \begin{array}{c c c c c}
                 1/f_1  & 0      &   \cdots  &  0         \\
                 0      & 1/f_2  &   \cdots  &  0         \\
                 \vdots & \vdots &   \ddots  &  \vdots    \\
                 0      & 0      &   \cdots  &  1/f_p     \\
                 \end{array} \right]
            \mathbf{C}^\prime)^{-1}
            (\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma})
\end{displaymath}
{\footnotesize
\begin{itemize}
    \item $f_1, \ldots f_p$ are relative sample sizes: $f_j = n_j/n$
    \item $\mathbf{C}\boldsymbol{\beta}- \mathbf{t}$ is an \emph{effect}, a particular way in which the null hypothesis is wrong. It is naturally expressed in units of the common within-treatment standard deviation sigma, and in general there is no reasonable way to avoid it.
    \item Almost always, $\mathbf{t} = \mathbf{0}$.
    \item The non-centrality parameter is sample size times a quantity that is sometimes called ``effect size." 
    \item The idea is that effect size represents how wrong $H_0$ is.
%    \item Here, effect size is squared distance in a space that is scaled and stretched by the relative sample sizes -- an aspect of design. 
\end{itemize}
} % End size
\end{frame}



\begin{frame}
\frametitle{Example: Comparing two means}
%\framesubtitle{} 
Suppose we have a random sample of size $n_1$ from a normal distribution with mean $\mu_1$ and variance $\sigma^2$, and independently, a second random sample from a normal distribution with mean $\mu_2$ and variance $\sigma^2$. We wish to test $H_0: \mu_1 = \mu_2$ versus the alternative $H_a: \mu_1 \neq \mu_2$. If the true means are a half a standard deviation apart, we want to be able to detect it with probability 0.80.

\vspace{5mm}
\noindent
We'll do it with cell means coding, letting $x_{i,1}=1$ if observation $i$ is from treatment one (and zero otherwise), and $x_{i,2}=1$ if observation $i$ is from treatment two (and zero otherwise). 

\begin{itemize}
    \item The model is $Y_i = \beta_1x_{i,1} + \beta_2 x_{i,2} + \epsilon_i$
    \item $\beta_1 = \mu_1$, $\beta_2 = \mu_2$
    \item $H_0: \beta_1-\beta_2=0$
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{$H_0: \mathbf{C}\boldsymbol{\beta}=\mathbf{t}$}
%\framesubtitle{} 
\begin{eqnarray*}
    \mathbf{C}         & = & (1, -1) \\
    \boldsymbol{\beta} & = & \left( 
        \begin{array}{c} \beta_1 \\ \beta_2  \end{array} \right) \\
    t                  & = & (0) \\
    \lambda            & = &  n  \times \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma}\right)^\prime
            (\mathbf{C}
                 \left( \begin{array}{c c}
                 1/f_1  & 0        \\
                 0      & 1/(1-f_1)  
                 \end{array} \right)
            \mathbf{C}^\prime)^{-1}
          \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma}\right) \\
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Continuing the calculation}
%\framesubtitle{} 
\begin{eqnarray*}
    \lambda            & = &  n  ~ \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma}\right)^\prime
            (\mathbf{C}
                 \left( \begin{array}{c c}
                 1/f_1  & 0        \\
                 0      & 1/(1-f_1)  
                 \end{array} \right)
            \mathbf{C}^\prime)^{-1}
          \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma}\right) \\
& = &   n ~ \left( \frac{\mu_1-\mu_2}{\sigma}\right)
        \left( \frac{1}{f_1} + \frac{1}{1-f_1} \right)^{-1}
        \left( \frac{\mu_1-\mu_2}{\sigma}\right) \\
& = &   n ~ f_1(1-f_1) \left( \frac{\mu_1-\mu_2}{\sigma}\right)^2 \\
& = &   n ~ f_1(1-f_1) ~ d^2
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{$ \lambda  = n f (1-f) d^2$, where $f = \frac{n_1}{n}$ and $d =  \frac{|\mu_1-\mu_2|}{\sigma}$}
%\framesubtitle{} 

\begin{itemize}
    \item For two-sample problems, $d$ is usually called effect size. The effect size specifies how wrong the null hypothesis is, by expressing the absolute difference between means in units of the common within-cell standard deviation.
    \item The non-centrality parameter (and hence, power) depends on the three parameters $\mu_1$, $\mu_2$ and $\sigma^2$ only through the effect size $d$.
    \item Power depends on sample size, effect size and an aspect of design – allocation of relative sample size to treatments. Equal sample sizes yield the highest power in the 2-sample case.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Back to the problem}
\framesubtitle{$ \lambda  = n f (1-f) d^2$} 
We wish to test $H_0: \mu_1 = \mu_2$ versus the alternative $H_a: \mu_1 \neq \mu_2$. If the true means are a half a standard deviation apart, we want to be able to detect it with probability 0.80.

\begin{eqnarray*}
    \lambda & = & n f (1-f) \left(\frac{|\mu_1-\mu_2|}{\sigma} \right)^2 \\
            & = & n \frac{1}{2}\left(1-\frac{1}{2} \right) \left(\frac{1}{2} \right)^2 \\
            & = & \frac{n}{16}
\end{eqnarray*}
\end{frame}

\begin{frame}[fragile]
\frametitle{SAS \texttt{proc iml}}
%\framesubtitle{} 
{\scriptsize 
\begin{verbatim}
/*****************  fpow1.sas *********************/
options linesize=79 noovp formdlim='_' nodate;
title 'Two-sample power analysis';

proc iml; /* Replace alpha, q, p, d and wantpow below    */
     alpha = 0.05;  /* Signif. level for testing H0: C Beta = t  */
     q = 1;         /* Numerator df = # rows in C matrix         */
     p = 2;         /* There are p beta parameters               */
     d = 1/2;       /* d = |mu1-mu2|/sigma */
     wantpow = .80; /* Find n to yield this power                */
     power = 0; n = p; oneminus = 1-alpha; /* Initializing ... */
     do until (power >= wantpow);
        n=n+1 ;
        ncp = n * 1/4 * d**2;
        df2 = n-p;
        power = 1-probf(finv(oneminus,q,df2),q,df2,ncp);
     end;
     print alpha p q d wantpow;
     print "Required sample size is " n;
     print "For a power of " power;
/**********************************************************

\end{verbatim} 
} % End size
\end{frame}

\begin{frame}[fragile]
\frametitle{Output}
%\framesubtitle{} 
{\footnotesize % or scriptsize
\begin{verbatim}
 
                           Two-sample power analysis                          1

                   alpha         p         q         d   wantpow

                    0.05         2         1       0.5       0.8


                                                       n

                      Required sample size is        128


                                               power

                           For a power of  0.8014596


\end{verbatim} 
} % End size
\end{frame}

\begin{frame}
\frametitle{To do a power analysis for any factorial design}
\framesubtitle{$H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{t}$} 
All you need is a vector of relative sample sizes and a vector of numbers representing the differences between $\mathbf{C}\boldsymbol{\beta}$ and  $\mathbf{t}$ in units of sigma.

\begin{displaymath}
    \lambda = n  \times \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma}\right)^\prime
            (\mathbf{C}
                 \left[ \begin{array}{c c c c c}
                 1/f_1  & 0      &   \cdots  &  0         \\
                 0      & 1/f_2  &   \cdots  &  0         \\
                 \vdots & \vdots &   \ddots  &  \vdots    \\
                 0      & 0      &   \cdots  &  1/f_p     \\
                 \end{array} \right]
            \mathbf{C}^\prime)^{-1}
            \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{t}}{\sigma}\right)
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{Example: Test for interaction}
%\framesubtitle{} 
\begin{center}
\begin{tabular}{|c|c|c||c|}  \hline
            & \multicolumn{2}{c||}{Level of B} &             \\ \hline
Level of A  &      1      &      2            & Average     \\ \hline
1           & $\mu_{11}$ & $\mu_{12}$       & $\mu_{1.}$  \\ \hline
2           & $\mu_{21}$ & $\mu_{22}$       & $\mu_{2.}$  \\ \hline
3           & $\mu_{31}$ & $\mu_{32}$       & $\mu_{3.}$  \\ \hline \hline
Average     & $\mu_{.1}$ & $\mu_{.2}$       & $\mu_{..}$  \\ \hline
\end{tabular}
\end{center} 

$H_0: \mu_{11}-\mu_{12} = \mu_{21}-\mu_{22} = \mu_{31}-\mu_{32}$

\end{frame}

\begin{frame}
\frametitle{$H_0: \mu_{11}-\mu_{12} = \mu_{21}-\mu_{22} = \mu_{31}-\mu_{32}$}
%\framesubtitle{} 
\begin{center}
\begin{tabular}{cccc}  
 $\mathbf{C}$  & $\boldsymbol{\beta}$ &  $=$ & $\mathbf{t}$ \\ 
 &&&\\
 $\left( \begin{array}{r r r r r r}
                 1  & -1 & -1 &  1 &  0 &  0 \\
                 0  & 0 &  1 & -1 & -1 &  1 
                \end{array} \right)$
& $\left( \begin{array}{c}
                  \mu_{11} \\
                  \mu_{12} \\
                  \mu_{21} \\
                  \mu_{22} \\
                  \mu_{31} \\
                  \mu_{32} 
    \end{array} \right)$
& $=$ & $\left( \begin{array}{c} 0 \\ 0 \end{array} \right)$
\end{tabular}
\end{center} 
Suppose this null hypothesis is false in a particular way that we want to be able to detect.
\end{frame}

\begin{frame}
\frametitle{Null hypothesis is wrong}
%\framesubtitle{} 
{\small
Suppose that for $A=1$ and $A=2$, the population mean of $Y$ is a quarter of a standard deviation higher when $B=2$, but if $A=3$, the population mean of $Y$ is a quarter of a standard deviation higher for $B=1$. Of course there are infinitely many sets of means satisfying these constraints, even if they are expressed in standard deviation units. But they will all have the same effect size. One such pattern is the following.  

\begin{center}
\begin{tabular}{|c|r|r|}  \hline
            & \multicolumn{2}{c|}{Level of B} \\ \hline
Level of A  &      1      &      2     \\ \hline
1           &  0.000    &  0.250         \\ \hline
2           &  0.000    &  0.250         \\ \hline
3           &  0.000    & -0.250         \\ \hline \hline
\end{tabular} 
\end{center}

Sample sizes are all equal, and we want to be able to detect an effect of this magnitude with probability at least 0.80.  
} % End size
\end{frame}

\begin{frame}
\frametitle{All we need is true $\mathbf{C}\boldsymbol{\beta}$}
\framesubtitle{ $H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{0}$ is wrong.} 
%{\small
\begin{center}
\begin{tabular}{|c|r|r|}  \hline
            & \multicolumn{2}{c|}{Level of B} \\ \hline
Level of A  &      1      &      2     \\ \hline
1           &  0.000    &  0.250         \\ \hline
2           &  0.000    &  0.250         \\ \hline
3           &  0.000    & -0.250         \\ \hline \hline
\end{tabular} 
~~~~~~ $\mathbf{C}\boldsymbol{\beta} = 
\left( \begin{array}{r} 0 \\ -0.5 \end{array} \right)$
% } % End size
\end{center}
\end{frame}

\begin{frame}[fragile]
\frametitle{Matrix calculations with \texttt{proc iml}}
\framesubtitle{$\lambda = \frac{(\mathbf{C}\boldsymbol{\beta}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{t})}
           {\sigma^2}$} 
{\scriptsize 
\begin{verbatim}
/*****************  fpow2.sas *********************/
options linesize=79 noovp formdlim='_';
title 'Sample size calculation for the interaction example';

proc iml;
/********* Edit this input: Rows of matrices are separated by commas  ********/
     alpha = 0.05; wantpow = .80;
     f = {1,1,1,1,1,1};          /* Relative sample sizes       */
     C = { 1 -1 -1  1  0  0,       /* Contrast matrix             */
           0  0  1 -1 -1  1};
     eff = {0, 0.5};    /* In standard deviation units */
/*****************************************************************************/

\end{verbatim} 
} % End size
\end{frame}

\begin{frame}[fragile]
\frametitle{\texttt{fpow2.sas} continued}
\framesubtitle{$\lambda = \frac{(\mathbf{C}\boldsymbol{\beta}-\mathbf{t})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{t})}
           {\sigma^2}$} 
{\scriptsize 
\begin{verbatim}

     p = nrow(f) ; q = nrow(eff); f = f/sum(f);
     core = inv(C*inv(diag(f))*C`);
     effsize = eff`*core*eff;
     power = 0; n = p; oneminus = 1-alpha;  /* Initializing ...*/
     do until (power >= wantpow);
        n = n+1 ;
        ncp = n * effsize;
        df2 = n-p;
        power = 1-probf(finv(oneminus,q,df2),q,df2,ncp);
     end; /*  End Loop */
     print " ";
     print " Required sample size is " n " for a power of " power;
     print " " ;
/**********************************************************\end{verbatim} 
} % End size
\end{frame}



\begin{frame}[fragile]
\frametitle{Output}
%\framesubtitle{} 
{\scriptsize
\begin{verbatim}
             Sample size calculation for the interaction example             1

                                        


                                          n                      power

         Required sample size is        697  for a power of  0.8001726



\end{verbatim} 
} % End size

$697/6 = 116.1667$ and $117*6 = 702$, so a total of $n=702$ experimental units are needed for equal sample sizes.

\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}

# Make picture of the normal
z = seq(from=-3,to=3,by=0.1)
Density = dnorm(z)
# Found these details in help(plot.default)
plot(z,Density,type="l", xlab="", ylab="", frame.plot=F, axes=F)
lines(c(-3,3),c(0,0),type='l') # Draw the bottom
# Now draw the lines showing tail areas
x = 1; ht = dnorm(x)
lines(c(-x,-x),c(0,ht),type='l')
lines(c(x,x),c(0,ht),type='l')
# Writing below the line is out unless I raise the whole thing up because it's outside the lotting area.

# F power picture
df1=5; df2=100; alpha = 0.05; lambda=15
crit = qf(1-alpha,df1,df2)
truepower = 1-pf(crit,df1,df2,ncp=lambda); truepower
x = seq(from=0.01,to=10,by=0.05)
Density = df(x,df1,df2)  
d2 = df(x,df1,df2,ncp=lambda)
plot(x,Density,type='l')
lines(x,d2,lty=2)
lines(c(crit,crit),c(-1,df(crit,df1,df2,ncp=lambda)))
tstring = expression(paste("Power of the F test with ",lambda," = 15"))
title(tstring)





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

\begin{frame}
\frametitle{}
%\framesubtitle{} 
\begin{itemize}
    \item 
    \item 
    \item 
\end{itemize}
\end{frame}