%                       Power for Applied Stat I

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

\hypersetup{colorlinks,linkcolor=,urlcolor=red}

\usefonttheme{serif} % Looks like Computer Modern for non-math text -- nice!
\setbeamertemplate{navigation symbols}{} % Suppress navigation symbols
% \usetheme{Berlin} % Displays sections on top
\usetheme{Frankfurt}  % Displays section titles on top: Fairly thin but still swallows some material at bottom of crowded slides
%\usetheme{Berkeley}

\usepackage[english]{babel}
\usepackage{amsmath} % for binom
% \usepackage{graphicx} % To include pdf files!
% \definecolor{links}{HTML}{2A1B81}
% \definecolor{links}{red}
\setbeamertemplate{footline}[frame number] 

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

% Comment this out for handout
\title{Power and Sample Size for Linear Regresson\footnote{See last slide for copyright information.}}
\subtitle{STA442/2101 Fall 2016}
\date{} % To suppress date

\begin{document}

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

\begin{frame}
\frametitle{Statistical Power}
% \framesubtitle{}
\begin{itemize}
    \item Power is the probability of rejecting $H_0$ when $H_0$ is false. \pause
    \item[] 
    \item Choose the sample size so as to have decent power. \pause
    \item[] 
    \item Study the distribution of the test statistic when the null hypothesis is false.
\end{itemize}
\end{frame}

% ATS uses lambda

\begin{frame}
\frametitle{Non-central Chisquare} \pause
%\framesubtitle{} 
\begin{eqnarray*}
    Z \sim N(0,1)   & \Rightarrow & Z^2 \sim \chi^2(1) \\ \pause
&& \\
    Z \sim N(\mu,1)  \pause & \Rightarrow & Z^2 \sim  \chi^2_{nc}(1,\lambda = \mu^2)  \\ \pause
&& \\
    Z_i \stackrel{i.i.d}{\sim} N(0,1) & \Rightarrow & \sum_{i=1}^n Z_i^2 \sim \chi^2(n) \\ \pause 
&& \\
    Z_i \stackrel{ind}{\sim} N(\mu_i,1)  \pause & \Rightarrow & \sum_{i=1}^n Z_i^2 \sim \chi^2_{nc}(n,\lambda=\sum_{i=1}^n \mu_i^2) 
\end{eqnarray*} 

\end{frame}

\begin{frame}
\frametitle{$Z_i \stackrel{ind}{\sim} N(\mu_i,1)  \Rightarrow  \sum_{i=1}^n Z_i^2 \sim \chi^2_{nc}(n,\lambda=\sum_{i=1}^n \mu_i^2)$}
%\framesubtitle{} 
 
\pause 
{\Large
Let $Y_i \stackrel{ind}{\sim} N(\mu_i,\sigma^2)$. Then \pause 

\vspace{4mm}

\begin{itemize}
     \item[]$Z_i = \frac{Y_i}{\sigma} \sim N(\frac{\mu_i}{\sigma},1)$  \pause
     \item[]
     \item[] $\sum_{i=1}^n Z_i^2 \sim \chi^2_{nc}(n,\lambda=\frac{\sum_{i=1}^n \mu_i^2}{\sigma^2})$
\end{itemize}
} % End size
\end{frame}


\begin{frame}
\frametitle{Non-central $F$ distribution}  \pause
\begin{displaymath}
    \begin{array}{l}
    W_1 \sim \chi^2_{nc}(\nu_1,\lambda) \\ 
    W_2 \sim \chi^2(\nu_2) \\  
    W_1 \mbox{ and } W_2 \mbox{ independent} \\ \\  \pause
    F^* = \frac{W_1/\nu_1}{W_2/\nu_2} \sim F(\nu_1,\nu_2,\lambda) \\ \\  \pause
    E(F^*) = \frac{\nu_2(\nu_1+\lambda)}{\nu_1(\nu_2-2)} \mbox{ for } \nu_2 > 2 \\  \pause
    Var(F^*) = \left( \frac{2\nu_2^2}{\nu_1^3(\nu_2-2)^2(\nu_2-4)} \right)
         \left( \lambda^2 + (2\lambda+\nu_1)(\nu_1+\nu_2-2) \right) \\
         \mbox{~~~~~~~ for } \nu_2 > 4
\end{array} 
\end{displaymath}  \pause
Reduces to the usual central $F$ when $\lambda=0$.
\end{frame}

\begin{frame}
\frametitle{General Linear Model when $H_0$ is false}  \pause
\begin{displaymath}
    \begin{array}{l}
    \mathbf{Y}~=~\mathbf{X} \boldsymbol{\beta}~+~\boldsymbol{\epsilon} \\\\  \pause
    H_{0}:  \mathbf{L} \boldsymbol{\beta} =  \mathbf{h} \\\\   \pause
    F^* = \frac{(\mathbf{L}\widehat{\boldsymbol{\beta}}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\widehat{\boldsymbol{\beta}}-\mathbf{h})}
           {r \, MSE} \\\\  \pause
    F^* \sim F(r,n-p,\lambda) \mbox{ where} \\\\  \pause
    \lambda = \frac{(\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})}
           {\sigma^2}
\end{array} 
\end{displaymath}  
\end{frame}

\begin{frame}
\frametitle{It's like the proof of $(\mathbf{y}-\boldsymbol{\mu})^\top
                 \boldsymbol{\Sigma}^{-1}(\mathbf{y}-\boldsymbol{\mu}) \sim \chi^2 (p)$  \pause}
\framesubtitle{Recalling $Z_i \stackrel{ind}{\sim} N(\mu_i,1)  \Rightarrow  \sum_{i=1}^n Z_i^2 \sim \chi^2_{nc}(n,\lambda=\sum_{i=1}^n \mu_i^2)$}  \pause
  \pause
\begin{enumerate}
    \item $\widehat{\boldsymbol{\beta}}  \sim N_p\left(\boldsymbol{\beta}, 
                   \sigma^2 (\mathbf{X}^\top \mathbf{X})^{-1}\right)$ \pause
    \item $\mathbf{L}\widehat{\boldsymbol{\beta}}  \sim N_r\left(\mathbf{L}\boldsymbol{\beta}, 
                   \sigma^2 \mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top\right)$ \pause
    \item $\mathbf{L}\widehat{\boldsymbol{\beta}} - \mathbf{h}  \sim 
          N_r\left(\mathbf{L}\boldsymbol{\beta} - \mathbf{h}, 
          \sigma^2 \mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top\right)$ \pause
    \item Letting $\mathbf{V} = \sigma^2\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top$, \pause
    \item $\mathbf{L}\widehat{\boldsymbol{\beta}} - \mathbf{h}  \sim 
          N_r\left(\mathbf{L}\boldsymbol{\beta} - \mathbf{h}, \mathbf{V}\right)$ \pause
    \item $\mathbf{z} = \mathbf{V}^{-1/2}(\mathbf{L}\widehat{\boldsymbol{\beta}} - \mathbf{h})    \pause
          \sim N_r\left(\mathbf{V}^{-1/2}(\mathbf{L}\boldsymbol{\beta} - \mathbf{h}), 
          \mathbf{V}^{-1/2}\mathbf{V}\mathbf{V}^{-1/2}\right)$ \pause
    \item  $\mathbf{z} \sim N_r\left(\boldsymbol{\mu},\mathbf{I}_r\right)$, where 
           $\boldsymbol{\mu} = \mathbf{V}^{-1/2}(\mathbf{L}\boldsymbol{\beta} - \mathbf{h})$ \pause
    \item $\mathbf{z}^\top\mathbf{z} =  \pause
            \frac{1}{\sigma^2}(\mathbf{L}\widehat{\boldsymbol{\beta}}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\widehat{\boldsymbol{\beta}}-\mathbf{h})$ \pause
    \item $\mathbf{z}^\top\mathbf{z} =  \sum_{j=1}^r Z_j^2 \sim \chi^2_{nc}(n,\lambda=\sum_{j=1}^r \mu_j^2)$ \pause
    \item $\lambda = \boldsymbol{\mu}^\top\boldsymbol{\mu} =  \pause
    (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top \mathbf{V}^{-1} (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})$ \pause
    \item $\lambda =   \frac{(\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})}
           {\sigma^2}$
\end{enumerate}
\end{frame}

\begin{frame}
\frametitle{Power increases to one with sample size when $H_0$ is false} \pause
%\framesubtitle{} 
{\small
\begin{itemize}
    \item The non-central $F$ distribution is \emph{stochastically increasing} in $\lambda$. \pause
    \item Stochastically increasing means $Pr\{ F > x \}$ is an increasing function of $\lambda$ for each fixed  $x >0$.  \pause
    \item The sample size is concealed in the non-centrality parameter.  \pause
\begin{displaymath}
    \lambda = n \,\frac{(\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\frac{1}{n}\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})}
           {\sigma^2}
\end{displaymath}  

\pause

    \item So the non-centrality parameter increases with $n$ unless $\mathbf{X}$ is crazy. \pause
    \item This means $Pr\{ F > x \}$ is an increasing function of $n$. \pause
    \item The critical value $x$ is a function of $n-p$, but it's a \emph{decreasing} function. \pause
    \item Finally, $\lim_{\lambda \rightarrow \infty} Pr\{ F > x \} = 1$. \pause
%    \item[]
    \item All this means that power goes to one as $n$ goes to infinity, as long as the null hypothesis is false.  \pause That is, the $F$ test is consistent.
\end{itemize}
} % End size
\end{frame}

\begin{frame}
\frametitle{Comparing two means} \pause
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$. \pause

\vspace{5mm}
\noindent
We'll do it with dummy variable regression, letting $x_i=1$ if observation $i$ is from population one, and $x_i=0$ if observation $i$ is from population two. Use cell means coding: \pause
\begin{displaymath}
    Y_i = \mu_1 x_i + \mu_2 (1-x_i) + \epsilon_i 
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{$H_0: \mathbf{L} \boldsymbol{\beta} =  \mathbf{h}$}
\framesubtitle{$Y_i = \mu_1 x_i + \mu_2 (1-x_i) + \epsilon_i $}  \pause
\begin{displaymath}
\begin{array}{cccc}
\mathbf{L} &  \boldsymbol{\beta} & = & \mathbf{h} \\
\left( \begin{array}{c c} 1 & -1
    \end{array} \right) &
\left( \begin{array}{c}
                 \mu_1 \\
                 \mu_2
    \end{array} \right) & = & 0
\end{array}
\end{displaymath} \pause

\begin{displaymath}
    \mathbf{X}^\top \mathbf{X} 
    = \left( \begin{array}{c c}
                 n_1 & 0 \\
                 0   & n_2
    \end{array} \right)  
~~~~~~
(\mathbf{X}^\top \mathbf{X})^{-1} = 
     \left( \begin{array}{c c}
                 1/n_1 & 0 \\
                 0   & 1/n_2
    \end{array} \right) 
\end{displaymath} \pause

\begin{eqnarray*}
    \lambda & = & \frac{1}{\sigma^2} (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h}) \\ \pause
     & = & \frac{1}{\sigma^2}(\mu_1-\mu_2)\, 
           \left(\frac{1}{n_1} + \frac{1}{n_2}\right)^{-1} \, 
           (\mu_1-\mu_2) 
\end{eqnarray*} 
\end{frame}

\begin{frame}
\frametitle{Non-centrality parameter}
%\framesubtitle{} 
\begin{eqnarray*}
    \lambda & = &  \frac{1}{\sigma^2}(\mu_1-\mu_2)\, 
           \left(\frac{1}{n_1} + \frac{1}{n_2}\right)^{-1} \, 
           (\mu_1-\mu_2) \\ \pause
     & = & \frac{(\mu_1-\mu_2)^2}{\sigma^2}  \left( \frac{n_1 n_2}{n_1+n_2} \right)\\ \pause
          & = & n \left( \frac{\mu_1-\mu_2}{\sigma} \right)^2
             \frac{n_1}{n}  \frac{n_2}{n} \\  \pause
          & = & n f (1-f) \left( \frac{\mu_1-\mu_2}{\sigma} \right)^2 \\ \pause
          & = & n f (1-f) d^2, 
\end{eqnarray*} 
where $f=\frac{n_1}{n}$, the fraction of observations in group one.
\end{frame}

\begin{frame}
\frametitle{Effect size}  \pause
% \framesubtitle{}
\begin{center}
    $\lambda = n f (1-f) d^2$, where $f = \frac{n_1}{n}$ and $d =  \frac{|\mu_1-\mu_2|}{\sigma}$
\end{center} \pause
\begin{itemize}
    \item $d$ is 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. \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. Equal sample sizes yield the highest power in the 2-sample case.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{How to proceed}  \pause
% \framesubtitle{}
\begin{itemize}
    \item Pick an effect size you'd like to be able to detect.  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.   Increase the sample size until the desired power is reached.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{For the 2-sample comparison}
%\framesubtitle{} 
\begin{itemize}
    \item Suppose we want to be able to detect a half standard deviation difference between means with power = 0.80 at the $\alpha = 0.05$ significance level. \pause
    \item Definitely use equal sample sizes. \pause
    \item $\lambda = n f (1-f) d^2 = n \left(\frac{1}{2}\right) \left(1-\frac{1}{2}\right)\left(\frac{1}{2}\right)^2 = n/16$
\end{itemize}
\end{frame}

 \begin{frame}[fragile]
\frametitle{Two sample test with R}
%\framesubtitle{} 
{\footnotesize % or scriptsize
\begin{verbatim}
> n <- seq(from=120,to=140,by=2) ; lambda <- n/16 ; ddf <- n-2
> fcrit = qf(.95,1,ddf) # Critical value for each n
> cbind(n,1-pf(fcrit,1,ddf,lambda))
        n
 [1,] 120 0.7752659
 [2,] 122 0.7820745
 [3,] 124 0.7887077
 [4,] 126 0.7951683
 [5,] 128 0.8014596
 [6,] 130 0.8075844
 [7,] 132 0.8135460
 [8,] 134 0.8193475
 [9,] 136 0.8249920
[10,] 138 0.8304825
[11,] 140 0.8358223

\end{verbatim} 
} % End size
 \pause  
$n_1=n_2=64$
\end{frame}
    
\begin{frame}
\frametitle{One Factor ANOVA ($p$ means)}
%\framesubtitle{} 
\begin{displaymath}
    \lambda = \frac{\sum_{k=1}^p n_k (\mu_k-\mu_.)^2}{\sigma^2}
                 ~~~ \mbox{ where } \mu_. = \sum_{k=1}^p \frac{n_k}{n} \mu_k
\end{displaymath}   \pause

\begin{eqnarray*}
     \lambda &=& \frac{\sum_{k=1}^p n_k (\mu_k-\mu_.)^2}{\sigma^2} \\ \pause
          &=&  n \, \frac{\sum_{k=1}^p \frac{n_k}{n} (\mu_k-\mu_.)^2}
                {\sigma^2} \\ \pause
          &=&  n \, \frac{\sum_{k=1}^p f_k (\mu_k-\mu_.)^2}
                {\sigma^2} \\ \pause
          &=&  n \, \sum_{k=1}^p f_k \left(\frac{\mu_k-\mu_.}{\sigma}\right)^2
\end{eqnarray*}  \pause
{\footnotesize
Notice how division and multiplication by $n$ helps. % It allows us to write the noncentrality parameter as the product of the sample size and a quantity that is either a constant or settles down rapidly to a constant as the sample size increases.
} % End size
\end{frame}
    
\begin{frame}
\frametitle{$\lambda = n \, \sum_{k=1}^p f_k \left(\frac{\mu_k-\mu_.}{\sigma}\right)^2$}  \pause
%\framesubtitle{} 
\begin{itemize}
    \item In the two-sample case, the non-centrality parameter was the product of sample size, effect size and the configuration of relative sample sizes.  \pause
    \item Once we get beyond two groups, effect and design are mixed together in a way that's impossible to separate.   \pause
    \item For a fixed sample size, $\lambda$ (and hence power) is maximized by splitting the sample equally between the two treatments whose means are farthest apart, and giving zero observations to the other treatments.   \pause
    \item Reluctantly, we will still call $\lambda/n$ ``effect size"  \pause
    \begin{itemize}
        \item Even though it does not reduce to what we called effect size before, if $p = 2$. 
          \pause It's $d^2/4$.  \pause
        \item And it is ``size" in a metric strongly influenced by the allocation of relative sample size to treatments.
    \end{itemize}
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Interpretation as signal to noise ratio} \pause
%\framesubtitle{} 
{\LARGE
\begin{eqnarray*}
    \lambda & = &  n \, \sum_{k=1}^p f_k \left(\frac{\mu_k-\mu_.}{\sigma}\right)^2 \\ \pause
            && \\ && \\
            & = & n \, \left(\frac{\sum_{k=1}^p f_k (\mu_k-\mu_.)^2}{\sigma^2} \right)
\end{eqnarray*} 
} % End size
\end{frame}


\begin{frame}[fragile]
\frametitle{Example using $\lambda = n \, \sum_{k=1}^p f_k \left(\frac{\mu_k-\mu_.}{\sigma}\right)^2$} \pause
Suppose we have four treatments, and that the four population treatment means are equally spaced, one-quarter of a standard deviation apart. \pause
We'd like to be able to detect the differences among treatment means with probability 0.80, using the conventional significance level of $\alpha=0.05$. Use equal sample sizes. \pause

Without loss of generality, we'll let the four population treatment means be 0, 0.25, 0.50 and 0.75. Using R as a calculator, and remembering that the \texttt{var} function divides by the number of observations minus one, we'll calculate the effect size as  \pause
\begin{verbatim}
> 3 * var(c(0,.25,.5,.75)) / 4
[1] 0.078125
\end{verbatim}
\end{frame}

\begin{frame}[fragile]
\frametitle{An R function}
%\framesubtitle{} 
{\scriptsize
\begin{verbatim}
fpow3 <- function(p,r,effsize,wantpow=0.80,alpha=0.05)
#############################################################################
# Power for the general multiple regression model, testing H0: L Beta = h   #
#       p       is the number of beta parameters                            #
#       r       Number rows in the L matrix = numerator df                  #
#       effsize is ncp/n, a squared distance between L Beta and h           #
#       wantpow is the desired power, default = 0.80                        #
#       alpha   is the significance level, default = 0.05                   #
#############################################################################
    {
    pow <- 0 ; nn <- p+1 ; oneminus <- 1 - alpha
    while(pow < wantpow)
        {
        nn <- nn+1
        lambda <- nn * effsize
        ddf <- nn-p
        pow <- 1 - pf(qf(oneminus,r,ddf),r,ddf,lambda)
        }#End while
    return(nn)
    }      # End of function fpow3
\end{verbatim}
 
} % End size
\end{frame}

\begin{frame}[fragile]
\frametitle{Finding sample size with \texttt{fpow3}}
%\framesubtitle{} 
{\footnotesize % or scriptsize
\begin{verbatim}
> source("http://www.utstat.utoronto.ca/~brunner/Rfunctions/fpow3.txt")
> fpow3(p=4,r=3,effsize=0.078125)
[1] 144
\end{verbatim} 
} % End size
\end{frame}

\begin{frame}
\frametitle{The Substitution Method}
%\framesubtitle{} 
\begin{itemize}
    \item Does the numerator of $\lambda = \frac{\sum_{k=1}^p n_k (\mu_k-\mu_.)^2}
                {\sigma^2}$  look familiar?  \pause
    \item It’s the standard elementary formula for the Between-Groups sum of squares in a one-way ANOVA, except with $\mu$ values substituted for sample means.  \pause
    \item This happens because the general formulas for $F$ and $\lambda$ are so similar.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{Formulas are similar}
\framesubtitle{Basis of the substitution method}  \pause
\begin{eqnarray*}
    F^* & = & (\mathbf{L}\widehat{\boldsymbol{\beta}}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\widehat{\boldsymbol{\beta}}-\mathbf{h}) / (r MSE)\\
    \lambda^{~} & = & (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h}) / \sigma^2
\end{eqnarray*} \pause
{\footnotesize
\begin{itemize}
    \item Any re-expression of the numerator of $F^*$ in terms of the sample cell means corresponds to a re-expression of the numerator of $\lambda$ in terms of population cell means.  \pause
    \item So, to obtain a formula for the non-centrality parameter, all you have to do is locate a convenient formula for the $F$-test of interest.   \pause In the expression for the numerator sum of squares, replace sample cell means by population cell means. Then divide by $\sigma^2$. The result is a formula for the non-centrality parameter.  \pause
    \item This applies to any $F$-test in any fixed effects factorial ANOVA.  \pause
    \item See Scheff\'e (1959), page 39 for a more general version of this rule.
\end{itemize}
} % End size
\end{frame}

\begin{frame}
\frametitle{Example: a 2-factor design}
%\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}   \pause

For equal sample sizes, I found the formula
\begin{displaymath}
    SSAB = n \sum_i \sum_j 
(\overline{Y}_{ij.}-\overline{Y}_{i..}-\overline{Y}_{.j.}+\overline{Y}_{...})^2 
\end{displaymath}  \pause
Which yields
\begin{displaymath}
    \lambda = \frac{n}{6} \frac{  \sum_{i=1}^3 \sum_{j=1}^2
(\mu_{ij}-\mu_{i.}-\mu_{.j}+\mu_{..})^2  } {\sigma^2}
\end{displaymath}  \pause

Different $n$~!
\end{frame}

\begin{frame}
\frametitle{How to use $\lambda = \frac{n}{6} \frac{  \sum_{i=1}^3 \sum_{j=1}^2
(\mu_{ij}-\mu_{i.}-\mu_{.j}+\mu_{..})^2  } {\sigma^2}$?}
\framesubtitle{}   \pause
\begin{itemize}
    \item What is a meaningful effect size? As far as I can tell, the only solution is to make up a meaningful effect, and apply the formula to it.  \pause
    \item In general, special purpose formulas may yield insight, but maybe not.  \pause
    \item Locating a special-purpose formula can be time consuming.   \pause
    \item You have to be sure of the notation, too.   \pause
    \item It can require some calculator work or a little programming.  Errors are possible.   \pause
    \item Often, a matrix approach is better, especially if you have to make up an effect and calculate its size anyway.
\end{itemize}
\end{frame}

\begin{frame}
\frametitle{An approach for factorial designs}
\framesubtitle{Use $\lambda = (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})/\sigma^2$}  \pause
\begin{itemize}
    \item Use cell means coding with one indicator for each treatment combination. \pause
    \item There are $p$ treatment combinations. \pause
    \item All the usual tests are tests of contrasts. \pause
    \item The $\mathbf{X}$ matrix has exactly one 1 in each row, and all the rest zeros. \pause
    \item There are $n_j$ ones in column $j$. \pause
\end{itemize}
\begin{displaymath}
    \mathbf{X}^\top \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) =   \pause n 
     \left( \begin{array}{c c c c c}
                 f_1    & 0      &   \cdots  &  0         \\
                 0      & f_2    &   \cdots  &  0         \\
                 \vdots & \vdots &   \ddots  &  \vdots    \\
                 0      & 0      &   \cdots  &  f_p       \\
    \end{array} \right) 
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{Multiplying and dividing by n}
\framesubtitle{ In $\lambda = (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})^\top
            (\mathbf{L}(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{L}^\top)^{-1}
            (\mathbf{L}\boldsymbol{\beta}-\mathbf{h})/\sigma^2$}  \pause
{\footnotesize
\begin{displaymath}
    \lambda  = n  \times (\frac{\mathbf{L}\boldsymbol{\beta}-\mathbf{h}}{\sigma})^\top
            (\mathbf{L}
                 \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_r     \\
                 \end{array} \right)
            \mathbf{L}^\top)^{-1}
            (\frac{\mathbf{L}\boldsymbol{\beta}-\mathbf{h}}{\sigma})
\end{displaymath}  \pause
} % End size
\begin{itemize}
    \item $f_1, \ldots f_p$ are relative sample sizes: $f_j = n_j/n$. \pause
    \item As usual, the non-centrality parameter is sample size times a quantity that we reluctantly call effect size. \pause
    \item $\mathbf{L}\boldsymbol{\beta}-\mathbf{h}$ is an \emph{effect} -- a particular way in which the null hypothesis is wrong.  \pause
    \item The effect is naturally expressed in units of the common within-treatment standard deviation $\sigma$, and in general this is hard to avoid. \pause
    \item Almost always, $\mathbf{h} = \mathbf{0}$.
\end{itemize}
\end{frame}


\begin{frame}
\frametitle{To actually do a power analysis} \pause
%\framesubtitle{} 
\begin{itemize}
    \item All you need is a vector of relative sample sizes,  \pause
    \item The contrast matrix $\mathbf{L}$,  \pause
    \item And a vector of numbers representing the differences between $\mathbf{L}\boldsymbol{\beta}$ and $\mathbf{h}$ in units of $\sigma$.
\end{itemize}  \pause

\vspace{4mm}

{\footnotesize
\begin{displaymath}
    \lambda  = n  \times (\frac{\mathbf{L}\boldsymbol{\beta}-\mathbf{h}}{\sigma})^\top
            (\mathbf{L}
                 \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_r     \\
                 \end{array} \right)
            \mathbf{L}^\top)^{-1}
            (\frac{\mathbf{L}\boldsymbol{\beta}-\mathbf{h}}{\sigma})
\end{displaymath} 
} % End size
\end{frame}
            
\begin{frame}
\frametitle{Recall the two-factor 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}

 \pause
 
\vspace{4mm}
$H_0: \mu_{11}-\mu_{12} = \mu_{21}-\mu_{22} = \mu_{31}-\mu_{32}$ 
\end{center} \pause

\begin{displaymath}
    \mathbf{L} = 
         \left( \begin{array}{r r r r r r}
                 1  & -1 & -1 &  1 &  0 &  0 \\
                 0  & 0 &  1 & -1 & -1 &  1 
                \end{array} \right) 
~~~~~
    \boldsymbol{\beta} = 
    \left( \begin{array}{c}
                  \mu_{11} \\
                  \mu_{12} \\
                  \mu_{21} \\
                  \mu_{22} \\
                  \mu_{31} \\
                  \mu_{32} 
    \end{array} \right)
\end{displaymath}

\end{frame}

\begin{frame}
\frametitle{A numerical example}
%\framesubtitle{} 
{\footnotesize
Suppose that for A=1 and A=2, the population mean of $Y$ is a quarter of a standard deviation higher for 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.  \pause

\vspace{3mm}

\begin{center}
\begin{tabular}{|r|r|r||r|}  \hline
            & \multicolumn{2}{c||}{Level of B} & \\ \hline
Level of A  &      1      &      2     & Average \\ \hline
1           &  0.000    &  0.250       &  0.125  \\ \hline
2           &  0.000    &  0.250       &  0.125  \\ \hline
3           &  0.000    & -0.250       & -0.125  \\ \hline \hline
Average     &  0.000    &  0.083       &  0.042  \\ \hline
\end{tabular}
\end{center}  \pause
\begin{itemize}
    \item $\frac{ \mathbf{L}\boldsymbol{\beta}-\mathbf{h}} {\sigma} = \left( 
        \begin{array}{c} 0 \\ 0.5 \end{array} \right)$  \pause
    \item Cell 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{itemize}
} % End size
\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/appliedf16} {\small\texttt{http://www.utstat.toronto.edu/$^\sim$brunner/oldclass/appliedf16}}

\end{frame}


\end{document}


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



{\LARGE
\begin{displaymath}
    
\end{displaymath}
 } % End Size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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