% Power for STA441, cut down from STA305s14, but also a bit enhanced.

% \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 by Statistical Power\footnote{See last slide for copyright information.}}
\subtitle{STA441 Spring 2020}
\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.
  \end{itemize}
\end{frame}

\section{Introduction}

\begin{frame}
\frametitle{How many subjects?} \pause
\begin{itemize}
     \item Mostly, we analyze data sets that people give us. The sample size is given. \pause
     \item It's better to decide in advance on some rational basis. \pause
     \item So how many experimental units do you need? \pause
     \item Need for what? \pause
     \item Maybe to make a parameter estimate accurate within some designated range. 
     \item Maybe to have a high probability of significant results when a relationship between variables is present.
\end{itemize}
\end{frame}


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

\section{Power}

\begin{frame}
\frametitle{Testing (null) hypotheses}
%\framesubtitle{} 
\begin{itemize}
    \item Goal is to make correct decisions with high probability. \pause
    \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. \pause
    \item When $H_0$ if false, we want to reject it with high probability. \pause
    \item The probability of rejecting the null hypothesis when the null hypothesis is false is called the \emph{power} of the test. \pause
    \item Power is one minus the probability of a Type II error. \pause
    \item It is a function of the true parameter values. \pause
    \item And also the design, including total sample size. 
\end{itemize}
\end{frame}

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

\begin{frame}
\frametitle{Power is an increasing function of sample size} \pause
%\framesubtitle{} 
\begin{itemize}
    \item Usually, when $H_0$ is false, larger sample size yields larger power. \pause
    \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}. \pause
    \item Most commonly used tests are consistent, including the general linear $F$-test.  \pause
    \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} \pause
%\framesubtitle{} 
\begin{itemize}
    \item Pick an effect you'd like to be able to detect. \pause An ``effect" means a way that $H_0$ is wrong. It should be just over the boundary of interesting and meaningful. \pause
    \item Pick a desired power -- a probability with which you'd like to be able to detect the effect by rejecting the null hypothesis. \pause
    \item Start with a fairly small $n$ and calculate the power. \pause Increase the sample size until the desired power is reached.
\end{itemize}
\end{frame}

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

\begin{frame}
\frametitle{Distribution theory} \pause
%\framesubtitle{} 
\begin{itemize}
    \item Power depends on the distribution of the test statistic when the null hypothesis is \emph{false}. \pause
    \item All the distributions you've seen ($Z$, $t$, $\chi^2$, $F$) were derived under the assumption that $H_0$ is \emph{true}. \pause
    \item For fixed-effects regression and analysis of variance, the distributions under the \emph{alternative} hypothesis are called \emph{non-central}. \pause
                \begin{itemize}
                    \item The non-central chi-squared.
                    \item The non-central $t$.
                    \item The non-central $F$.
                \end{itemize} \pause
    \item These distributions have the same parameters as their ordinary (central) versions, and they also have a \emph{non-centrality parameter}. \pause
    \item If the non-centrality parameter equals zero, the non-central distribution reduces to the ordinary central distribution.

\end{itemize}
\end{frame}

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

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

\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$. \pause
    \item That is, the bigger the non-centrality parameter, the greater the probability of getting $F^*$ above any point (such as a critical value). \pause
    \item $\displaystyle \lim_{\lambda \rightarrow \infty} Pr\{F^*>x|\lambda \} = 1$. 
\end{itemize}
\end{frame}

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

\begin{frame}
\frametitle{The greater the non-centrality parameter $\lambda$, the greater the power}
\framesubtitle{$\lambda=0$ means the null hypothesis is true} \pause 
\begin{center}
\includegraphics[width=3.5in]{Fpower}
\end{center}
\end{frame}

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

\begin{frame}
\frametitle{Non-centrality parameter looks like the test statistic}
\framesubtitle{But with true parameter values instead of estimates}  \pause
For Testing $H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{h}$\pause, like  \pause
\begin{displaymath}
    \left(
    \begin{array}{cccc}
    0 & 1 & 0 & 0 \\
    0 & 0 & 1 & 0 \\
    \end{array}
    \right)
    \left(
    \begin{array}{c}
    \beta_0 \\ \beta_1 \\ \beta_2 \\ \beta_3
    \end{array}
    \right) = 
    \left(
    \begin{array}{c}
    0 \\ 0
    \end{array}
    \right)
\end{displaymath} \pause
We have
\begin{eqnarray*}
    F & = & \frac{(\mathbf{C} \widehat{\boldsymbol{\beta}}-\mathbf{h})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C} \widehat{\boldsymbol{\beta}}-\mathbf{h})}
           {r \, MSE} \\
    && \\
    \lambda & = & 
            \frac{(\mathbf{C}\boldsymbol{\beta}-\mathbf{h})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{h})}
           {\sigma^2} 
\end{eqnarray*}
\end{frame}

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

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

\vspace{4mm}

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

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

\begin{frame}
\frametitle{An Important Special Case: Factorial ANOVA with no covariates} \pause
%\framesubtitle{} 
\begin{itemize}
    \item Use cell means coding. \pause
    \item Skipping a lot of STA305 material \ldots
\end{itemize}
\end{frame}

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

\begin{frame}
\frametitle{For cell means coding, $\lambda$ simplifies to}\pause
%\framesubtitle{}  

\begin{displaymath}
    \lambda = n  \times (\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{h}}{\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{h}}{\sigma})
\end{displaymath} \pause
{\footnotesize
\begin{itemize}
    \item $f_1, \ldots f_p$ are relative sample sizes: $f_j = n_j/n$. \pause
    \item $\mathbf{C}\boldsymbol{\beta}- \mathbf{h}$ is an \emph{effect}, a particular way in which the null hypothesis is wrong. \pause It is naturally expressed in units of the common within-treatment standard deviation $\sigma$, and in general there is no reasonable way to avoid this. \pause
    \item The magnitude of a difference does not matter, because it depends on the units. \pause What counts is the size of the difference relative to how spread out the data are.
\end{itemize}
} % End size
\end{frame}

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

\begin{frame}
\frametitle{More about the non-centrality parameter}
%\framesubtitle{}  

\begin{displaymath}
    \lambda = n  \times (\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{h}}{\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{h}}{\sigma})
\end{displaymath} \pause
{\footnotesize
\begin{itemize}
    \item Almost always, $\mathbf{h} = \mathbf{0}$. \pause
    \item The non-centrality parameter is sample size times a quantity that is sometimes called ``effect size."  \pause
    \item The idea is that effect size represents how wrong $H_0$ is. \pause
    \item Here, effect size is squared distance between $\mathbf{C}\boldsymbol{\beta}$ and $\mathbf{h}$, 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} \pause
%\framesubtitle{} 
{\small
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_1: \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.  \pause
} % End size

\begin{center}
\includegraphics[width=2.5in]{halfSD}
\end{center}

\end{frame}

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

\begin{frame}
\frametitle{Two-sample $t$-test, or $F$-test}
\framesubtitle{Non-central $t$ or non-central $F$}  \pause
\begin{itemize}
    \item We'll use $F$. \pause
    \item Skipping the derivation, get
\end{itemize}
{\Large
\begin{displaymath}
    \lambda  = n f (1-f) d^2
\end{displaymath} \pause
} % End size

where $f = \frac{n_1}{n}$ and $d =  \frac{|\mu_1-\mu_2|}{\sigma}$.
\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. \pause 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. \pause
    \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$. \pause
    \item Power depends on sample size, effect size and an aspect of design -- allocation of relative sample size to treatments.  \pause
    \item 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_1: \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. \pause

\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 *********************/
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}
\frametitle{Output}
%\framesubtitle{} 
\begin{center}
\includegraphics[width=2.5in]{fpow1}
\end{center}
\end{frame}

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

\begin{frame}
\frametitle{To do a power analysis for any factorial design}
\framesubtitle{$H_0: \mathbf{C}\boldsymbol{\beta} = \mathbf{h}$} \pause
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{h}$ in units of $\sigma$. \pause

\begin{displaymath}
    \lambda = n  \times \left(\frac{\mathbf{C}\boldsymbol{\beta}-\mathbf{h}}{\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{h}}{\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{h}$ \\ 
 &&&\\
 $\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}  \pause
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$. \pause Of course there are infinitely many sets of means satisfying these constraints, even if they are expressed in standard deviation units. \pause But they will all have the same effect size. \pause One such pattern is the following.   \pause

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

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}  \pause
~~~~~~ $\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{h})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{h})}
           {\sigma^2}$}  \pause
{\scriptsize 
\begin{verbatim}
/*****************  fpow2.sas *********************/
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{h})^\prime
            (\mathbf{C}(\mathbf{X}^\prime \mathbf{X})^{-1}\mathbf{C}^\prime)^{-1}
            (\mathbf{C}\boldsymbol{\beta}-\mathbf{h})}
           {\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 " Required sample size is " n " for a power of " power;
/**********************************************************\end{verbatim} 
} % End size
\end{frame}

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

\begin{frame}[fragile]
\frametitle{Output}
%\framesubtitle{} 
\begin{center}
\includegraphics[width=3.5in]{fpow2}
\end{center} \pause

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

\end{frame}

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

\begin{frame}
\frametitle{Beyond the $F$-tests} \pause
%\framesubtitle{} 
\begin{itemize}
    \item  Lots of large-sample chi-squared tests (like Wald and Likelihood Ratio) have limiting non-central chi-squared distributions under $H_1$.  
    \item Non-centrality parameters look like the test statistics. \pause
    \item For unfamiliar tests, simulation may be easier.
\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 Mathematical and Computational Statistics, University of Toronto Mississauga. 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/441s20} {\footnotesize \texttt{http://www.utstat.toronto.edu/$^\sim$brunner/oldclass/441s20}}

\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 plotting 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)

# Treatment means mu1 and mu2 half a standard deviation apart
rm(list=ls())
x = seq(from=-4,to=4,length=100)
Density = dnorm(x,mean=0.25); y = dnorm(x,mean=-0.25)
plot(x,Density,type='l',main='Expected values half a standard deviation apart')
lines(x,y,lty=1)




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

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