% \documentclass[serif]{beamer} % Get Computer Modern math font.
% \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{amsmath} % For \binom{n}{y}
\usepackage{graphicx} % To include pdf files!
\usepackage{fullpage}
\usepackage{comment}

\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{The Multinomial Model\footnote{See last slide for copyright information.}}
% \subtitle{STA 312: Fall 2012}
% \date{} % To suppress date


\begin{document}

% More material required for handout in article mode. Also eliminate vspace
\title{The Multinomial Model\footnote{See last page for copyright information.}}
\subtitle{STA 312: Fall 2022}
\date{} % To suppress date
\maketitle


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

% \section{Overview}

% \begin{frame}
% \frametitle{Overview}
% \tableofcontents
% \end{frame}

\section{Binomial Distribution}


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


\begin{frame}
\frametitle{The Bernoulli Distribution}
% \framesubtitle{Example 2.3.2} 

\begin{itemize}
    \item Simple probability model: Toss a coin with $P(\mbox{Head})=\pi$, one time. Let $Y$ equal the number of heads. 
    \item Probability (mass) function of $Y$: 
\begin{displaymath}
    P(y) = \left\{ \begin{array}{ll} % ll means left left
   \pi^y(1-\pi)^{1-y} & \mbox{for $y = 0$ or 1}  \\
   0     & \mbox{Otherwise} 
   \end{array}  \right.  % Need that crazy invisible right period!
\end{displaymath}  
    \item An \emph{indicator random variable} equals one if some event happens, and zero if it does not happen.  
        \begin{itemize}
            \item 1=Female, 0=Male 
            \item 1=Lived, 0=Died 
            \item 1=Passed, 0=Failed 
        \end{itemize}
    \item Indicators are usually assumed to have a Bernoulli distribution.
\end{itemize}
\end{frame}

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

\begin{frame}
\frametitle{The Binomial Distribution} 
% \framesubtitle{Example 2.3.3}   
\begin{itemize}
    \item Simple probability model: Toss a coin with $P(\mbox{Head})=\pi$.   Toss it $n$ times. Let $Y$ equal the number of heads.  
    \item Probability (mass) function of $X$:   
\begin{displaymath}
    P(y)  = \left\{ \begin{array}{ll} % ll means left left
   \binom{n}{y}\pi^y(1-\pi)^{n-y} & \mbox{for $y = 0, 1, \ldots, n$}  \\
   ~0     & \mbox{Otherwise} 
   \end{array}  \right.  % Need that crazy invisible right period!
\end{displaymath}   
    \item The Bernoulli is a special case of the Binomial, with $n=1$.
\end{itemize}
\end{frame}

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

\begin{frame} % Maybe omit or by-pass
\frametitle{Why does $P(y) = \binom{n}{y}\pi^y(1-\pi)^{n-y}$} 
\framesubtitle{For the Binomial Distribution?}    
Toss a coin $n$ times with $P(\mbox{Head})=\pi$, and let $Y$ equal the number of heads. Why does 
$P(Y=y) = \binom{n}{y}\pi^y(1-\pi)^{n-y}$? 
\begin{itemize}
    \item The sample space is the set of all strings of $n$ letters composed of H and T. 
    \item By the Multiplication Principle, there are $2^n$ elements.  % Think of a tree. 
    \item If two different strings have $y$ heads (and $n-y$ tails), they have the same probability. 
    \item For example, $P\{HHTH\} = P\{THHH\} = \pi^3(1-\pi)$ by independence. 
    \item Count the number of ways that $y$ positions out of $n$ can be chosen to have the symbol H. 
    \item $n$ choose $y$ is $\binom{n}{y} = \frac{n!}{y!(n-y)!}$. 
    \item So $P(Y=y) = \binom{n}{y}\pi^y(1-\pi)^{n-y}$ ~  $\blacksquare$
\end{itemize}
\end{frame}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Multinomial Distribution} %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Multinomial coefficient}
\framesubtitle{For $c$ categories}
From $n$ objects, number of ways to choose
  \begin{itemize}
    \item $n_1$ of type 1
    \item $n_2$ of type 2 \\
         ~~~~~~~$\vdots$
    \item $n_c$ of type $c$
  \end{itemize}

\begin{displaymath}
    \binom{n}{n_1~\cdots~n_c} = \frac{n!}{n_1!~\cdots~n_c!}
\end{displaymath}
\end{frame}

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

\begin{frame}
\frametitle{Example of a multinomial coefficient}
\framesubtitle{A counting problem} % Optional sub-title

Of 30 graduating students, how many ways are there for 15 to be employed in a job related to their field of study, 10 to be employed in a job unrelated to their field of study, and 5 unemployed?

\begin{displaymath}
    \binom{30}{15~10~5} = 465,817,912,560
\end{displaymath}
\end{frame}
% Sage: factorial(30)/(factorial(5)*factorial(25))

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


\begin{frame}
\frametitle{Multinomial Distribution}
\framesubtitle{Denote by $M(n,\boldsymbol{\pi})$, where 
$\boldsymbol{\pi} = (\pi_1, \ldots, \pi_c)$}

  \begin{itemize}
    \item Statistical experiment with $c$ outcomes
    \item Repeated independently $n$ times
    \item $Pr(\mbox{Outcome }j) = \pi_j$, $j = 1, \ldots, c$
    \item Number of times outcome $j$ occurs is $n_j$, $j = 1, \ldots, c$
    \item An integer-valued \emph{multivariate} distribution
  \end{itemize}
\begin{displaymath}
    P(n_1, \ldots, n_c) = \binom{n}{n_1~\cdots~n_c} \pi_1^{n_1} \cdots \pi_c^{n_c},
\end{displaymath}
\vspace{5mm}

where $0 \leq n_j \leq n$, $\displaystyle{\sum_{j=1}^c n_j=n}$, $0 < \pi_j < 1$, and 
$\displaystyle{\sum_{j=1}^c \pi_j=1}$.
\end{frame}

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

\begin{frame}
\frametitle{There are actually $c-1$ variables and $c-1$ parameters}
\framesubtitle{In the multinomial with $c$ categories}
\begin{eqnarray*}
    P(n_1, \ldots, n_{c-1}) &=& \frac{n!}{n_1!~\cdots~n_{c-1}!(n-\sum_{j=1}^{c-1}n_j)!} \\
& & ~~~~\times~ \pi_1^{n_1} \cdots \pi_{c-1}^{n_{c-1}}(1-\sum_{j=1}^{c-1}\pi_j)^{n-\sum_{j=1}^{c-1}n_j}
\end{eqnarray*}
\end{frame}

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

\begin{frame}
\frametitle{Marginals of the multinomial are multinomial too}
\framesubtitle{Add over $n_{c-1}$, which goes from zero to whatever is left over from the other counts.} 
% It would be good to shift this over to the left somehow.
{\Tiny
\begin{eqnarray*}
 & & \sum_{n_{c-1}=0}^{n-\sum_{j=1}^{c-2}n_j} 
\frac{n!}{n_1!~\ldots~n_{c-1}!(n-\sum_{j=1}^{c-1}n_j)!}
\pi_1^{n_1} \cdots \pi_{c-1}^{n_{c-1}}(1-\sum_{j=1}^{c-1}\pi_j)^{n-\sum_{j=1}^{c-1}n_j} 
\times \frac{(n-\sum_{j=1}^{c-2}n_j)!}{(n-\sum_{j=1}^{c-2}n_j)!} 
\\ \\ \\
 &=& \frac{n!}{n_1!~\ldots~n_{c-2}!(n-\sum_{j=1}^{c-2}n_j)!} 
\pi_1^{n_1} \cdots \pi_{c-2}^{n_{c-2}} \\
 & & \times
\sum_{n_{c-1}=0}^{n-\sum_{j=1}^{c-2}n_j} 
\frac{(n-\sum_{j=1}^{c-2}n_j)!}{n_{c-1}!(n-\sum_{j=1}^{c-2}n_j-n_{c-1})!}
\pi_{c-1}^{n_{c-1}} (1-\sum_{j=1}^{c-2}\pi_j-\pi_{c-1})^{n-\sum_{j=1}^{c-2}n_j-n_{c-1}} 
\\ \\ \\
 &=& \frac{n!}{n_1!~\ldots~n_{c-2}!(n-\sum_{j=1}^{c-2}n_j)!} 
\pi_1^{n_1} \cdots \pi_{c-2}^{n_{c-2}} (1-\sum_{j=1}^{c-2}\pi_j)^{n-\sum_{j=1}^{c-2}n_j},
\end{eqnarray*}
}

where the last equality follows from the Binomial Theorem. It's multinomial with $c-1$ categories.
%This is too messy -- students are not responsible for it.
\end{frame}

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

\begin{frame}
\frametitle{Observe}
\framesubtitle{You are responsible for these \emph{implications} of the last slide.} 

  \begin{itemize}
    \item Adding over $n_{c-1}$ throws it into the last (``leftover") category.
    \item Labels $1, \ldots, c$ are arbitrary, so this means you can combine any 2 categories and the result is still multinomial.
    \item $c$ is arbitrary, so you can keep doing it and combine any number of categories.
    \item When only two categories are left, the result is binomial
    \item $E(n_j)=n\pi_j = \mu_j$, $Var(n_j) = n\pi_j(1-\pi_j)$
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Sample problem}
\framesubtitle{Recent university graduates} 
  \begin{itemize}
    \item Probability of job related to field of study = 0.60
    \item Probability of job unrelated to field of study = 0.30
    \item Probability of no job = 0.10
  \end{itemize}

Of 30 randomly chosen students, what is probability that 15 are employed in a job related to their field of study, 10 are employed in a job unrelated to their field of study, and 5 are unemployed?
\begin{displaymath}
    \binom{30}{15~10~5} \, 0.60^{15} 0.30^{10} 0.10^5
= \frac{4933527642332542053801}{381469726562500000000000} \approx 0.0129
\end{displaymath}
% Sage: factorial(30)/(factorial(15)*factorial(10)*factorial(5)) * ((6/10)^(15)*(3/10)^(10)*(1/10)^5)

What is the probability that exactly 5 are unemployed?
\begin{displaymath}
    \binom{30}{5~25} \, 0.10^5 0.90^{25} = 
\frac{51152385317007572507646551997}{500000000000000000000000000000} \approx 0.1023
\end{displaymath}
% Sage factorial(30)/(factorial(5)*factorial(25)) * ((9/10)^(25)*(1/10)^5)
\end{frame}



\begin{frame}
\frametitle{Conditional probabilities are multinomial too}
%\framesubtitle{} 
  \begin{itemize}
    \item Given that a student finds a job, what is the probability that the job will be in the student's field of study?
\begin{displaymath}
    P(\mbox{Field}|\mbox{Job}) = \frac{P(\mbox{Field, Job})}{P(\mbox{Job})}
    = \frac{0.60}{0.90} = \frac{2}{3}
\end{displaymath}

    \item Suppose we choose 50 students at random from those who found jobs. What is the probability that exactly $y$ of them will be employed in their field of study, for $y = 0, \ldots, 50$?
\begin{displaymath}
    P(y|\mbox{Job}) = \binom{50}{y} \left(\frac{2}{3}\right)^y \left(1-\frac{2}{3}\right)^{50-y}
\end{displaymath}
  \end{itemize}
\end{frame}




\begin{frame}[fragile]
\frametitle{Calculating multinomial probabilities with R}
%\framesubtitle{} 
Of 30 randomly chosen students, what is probability that 15 are employed in a job related to their field of study, 10 are employed in a job unrelated to their field of study, and 5 are unemployed?
\begin{displaymath}
    \binom{30}{15~10~5} \, 0.60^{15} 0.30^{10} 0.10^5
= \frac{4933527642332542053801}{381469726562500000000000} \approx 0.0129
\end{displaymath}

\vspace{10mm}

\begin{verbatim}
> dmultinom(c(15,10,5), prob=c(.6, .3, .1))
[1] 0.01293295
\end{verbatim}
\end{frame}

\section{Estimation} %%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{Hypothetical data file}
\framesubtitle{Let $Y_{i,j}$ be indicators for category membership, $i=1, \ldots, n$ and $j=1, \ldots, c$}
\begin{center}
\begin{tabular}{|c|c||c|c|c|}  \hline
Case     &   Job    &  $Y_1$   &  $Y_2$    &  $Y_3$    \\ \hline \hline
  1      &    1     &    1     &    0      &    0      \\ \hline
  2      &    3     &    0     &    0      &    1      \\ \hline
  3      &    2     &    0     &    1      &    0      \\ \hline
  4      &    1     &    1     &    0      &    0      \\ \hline
$\vdots$ & $\vdots$ & $\vdots$ & $\vdots$  & $\vdots$  \\ \hline
 $n$     &    2     &    0     &    1      &    0      \\ \hline \hline
Total    &          & $\sum_{i=1}^ny_{i,1}$ & $\sum_{i=1}^ny_{i,2}$  & $\sum_{i=1}^ny_{i,3}$  \\ \hline
\end{tabular}
\end{center}
Note that
  \begin{itemize}
    \item A real data file will almost never have the redundant variables $Y_1$, $Y_2$ and $Y_3$.
    %\item $\sum_{j=1}^cy_{i,j}=1$
    \item $\sum_{i=1}^ny_{i,j}=n_j$
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Lessons from the data file}
%\framesubtitle{} 
  \begin{itemize}
    \item Cases ($n$ of them) are independent $M(1,\boldsymbol{\pi})$, so $E(Y_{i,j}) = \pi_j$.
    \item Column totals $n_j$ count the number of times each category occurs:  Joint distribution is $M(n,\boldsymbol{\pi})$
    \item If you make a frequency table (frequency distribution)
        \begin{itemize}
            \item The $n_j$ counts are the cell frequencies!
            \item They are random variables, and now we know their joint distribution.
            \item Each individual (marginal) table frequency is $B(n,\pi_j)$.
            \item Expected value of cell frequency j is $E(n_j) = n\pi_j = \mu_j$
        \end{itemize}
    \item Tables of 2 and or more dimensions present no problems; form combination variables.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Example of a frequency table}
\framesubtitle{For the Jobs data} 
\begin{center}
\begin{tabular}{lrr}  \hline
Job Category & Frequency & Percent      \\ \hline
Employed in field       & 106 & 53      \\
Employed outside field  &  74 & 37      \\
Unemployed              &  20 & 10      \\ \hline
Total                   & 200 & 100.0   \\ \hline
\end{tabular}
\end{center}
\end{frame}

\begin{frame}
\frametitle{Likelihood function for the multinomial}
\begin{eqnarray*}
    \ell(\boldsymbol{\pi}) 
        & = & \prod_{i=1}^n Pr\{Y_{i,1}=y_{i,1}, Y_{i,2}=y_{i,2}, \ldots, Y_{i,c}=y_{i,c}  |\boldsymbol{\pi}\} \\
        & = & \prod_{i=1}^n \pi_1^{y_{i,1}} \pi_2^{y_{i,2}} \cdots \pi_c^{y_{i,c}} \\
        & = & \pi_1^{\sum_{i=1}^n y_{i,1}} \pi_2^{\sum_{i=1}^n y_{i,2}} \cdots \pi_c^{\sum_{i=1}^n y_{i,c}} \\
        & = & \pi_1^{n_1} \pi_2^{n_2} \cdots \pi_c^{n_c}
\end{eqnarray*}

  \begin{itemize}
    \item Product of $n$ probability mass functions, each  $M(1,\boldsymbol{\pi})$
    \item Depends upon the sample data only through the vector of $c$ frequency counts: $(n_1, \ldots, n_c)$
  \end{itemize}
%    \item 
%    \item  
%  So it would be okay to think of the experiment as yielding a single $M(n,\boldsymbol{\pi})$ observation.
\end{frame}

\begin{frame}
\frametitle{All you need is the frequency table}
\begin{displaymath}
    \ell(\boldsymbol{\pi}) = \pi_1^{n_1} \pi_2^{n_2} \cdots \pi_c^{n_c}
\end{displaymath}
  \begin{itemize}
    \item Likelihood function depends upon the sample data only through the frequency counts.
    \item By the factorization theorem, $(n_1, \ldots, n_c)$ is a sufficient statistic.
    \item \emph{All} the information about the parameter in the sample data is contained in the  sufficient statistic.
    \item So everything the sample data could tell you about $(\pi_1, \ldots, \pi_c)$ is given by in $(n_1, \ldots, n_c)$. 
    \item You don't need the raw data.
  \end{itemize}
% Sufficiency is a theoretical concept with extremely practical implications: Data reduction
% But you might need the raw data if there were explanatory variables. 
\end{frame}

\begin{frame}
\frametitle{Log likelihood: $c-1$ parameters}
\begin{eqnarray*}
    \ell(\boldsymbol{\pi}) & = & \pi_1^{n_1} \cdots \pi_c^{n_c} \\
                           & = & \pi_1^{n_1} \cdots \pi_{c-1}^{n_{c-1}} 
                \left(1-\sum_{j=1}^{c-1} \pi_j\right)^{n-\sum_{j=1}^{c-1} n_j} \\
    \log\ell(\boldsymbol{\pi})  & = & \sum_{j=1}^{c-1} n_j \log\pi_j 
 + \left(n-\sum_{j=1}^{c-1} n_j\right) \log\left(1-\sum_{j=1}^{c-1} \pi_j\right)
\end{eqnarray*}

\vspace{10mm}

\begin{displaymath}
    \frac{\partial\log\ell}{\partial \pi_j} = 
\frac{n_j}{\pi_j} - \frac{n-\sum_{k=1}^{c-1} n_k}{1-\sum_{k=1}^{c-1} \pi_k},
\mbox{ for } j=1,\ldots,c-1
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{Set all the partial derivatives to zero and solve}
\framesubtitle{For $\pi_j$, j = 1 \ldots, c-1}

{\huge
\begin{displaymath}
    \widehat{\pi}_j = \frac{n_j}{n} = p_j = \frac{\sum_{i=1}^n y_{i,j}}{n} = \overline{y}_j
\end{displaymath}
}

\vspace{10mm}

So the MLE is the sample proportion, which is also a sample mean.
% HW: What about pi-hat_c?
\end{frame}

\begin{frame}
\frametitle{In matrix terms: $\widehat{\boldsymbol{\pi}} = \textbf{p} = \overline{\textbf{Y}}_n$}

\begin{displaymath}
\left( \begin{array}{c}
\widehat{\pi}_1 \\ 
\vdots \\ 
\widehat{\pi}_{c-1}
        \end{array} \right) = 
\left( \begin{array}{c}
p_1 \\ 
\vdots \\ 
p_{c-1}
        \end{array} \right) =
\left( \begin{array}{c}
\overline{Y}_{1~~} \\ 
\vdots \\ 
\overline{Y}_{c-1}
        \end{array} \right)
\end{displaymath}

%\vspace{5mm}

Remarks:

  \begin{itemize}
    \item Multivariate Law of Large Numbers says $\textbf{p} \stackrel{p}{\rightarrow} \boldsymbol{\pi}$
    \item Multivariate Central Limit Theorem says that $\overline{\textbf{Y}}_n$ is approximately multivariate normal for large $n$.
    \item Because $n_j \sim B(n,\pi_j)$, $\frac{n_j}{n} = \overline{Y}_j = p_j$ is approximately $N\left(\pi_j,\frac{\pi_j(1-\pi_j)}{n}\right)$.
    \item Approximate $\pi_j$ with $p_j$ in the variance if necessary.
    \item Can be used in confidence intervals and tests about a single parameter.
    \item We have been using $c-1$ categories only for technical convenience. 
  \end{itemize}
\end{frame}


\begin{frame}
\frametitle{Confidence interval for a single parameter $\pi_j$}
\framesubtitle{95\% confidence interval for true proportion unemployed} 
\begin{center}
{\tiny
\begin{tabular}{lrr}  \hline
Job Category & Frequency & Percent      \\ \hline
Employed in field       & 106 & 53      \\
Employed outside field  &  74 & 37      \\
Unemployed              &  20 & 10      \\ \hline
Total                   & 200 & 100.0   \\ \hline
\end{tabular}
}
\end{center}

\begin{displaymath}
    p_3 \stackrel{approx}{\sim} N\left(\pi_3,\frac{\pi_3(1-\pi_3)}{n}\right)
\end{displaymath}

So a confidence interval for $\pi_3$ is

\begin{eqnarray*}
    p_3 \pm 1.96\sqrt{\frac{p_3(1-p_3)}{n}} & = & 
            0.10 \pm 1.96\sqrt{\frac{0.10(1-0.10)}{200}} \\
            & = &  0.10 \pm 0.042 \\
            & = &  (0.058,0.142)
\end{eqnarray*}
% me = 1.96*sqrt(.1*.9/200); me
% .1-me; .1+me
\end{frame}


\section{Hypothesis tests} %%%%%%%%%%%%%%%%%%%%%%%%

\begin{frame}
\frametitle{For general tests on multinomial data}
We will use mostly
  \begin{itemize}
    \item Pearson chi-squared tests
    \item Large-sample likelihood ratio tests
  \end{itemize}

\vspace{10mm}

There are other possibilities, including
  \begin{itemize}
    \item Wald tests
    \item Score tests
  \end{itemize}

\vspace{10mm}

All these are large-sample chi-squared tests, justified as $n \rightarrow \infty$
\end{frame}

\begin{frame}
\frametitle{Likelihood ratio tests}
\framesubtitle{In general}


Setup
\begin{displaymath}
    \begin{array}{l}
     Y_1, \ldots, Y_n  \stackrel{i.i.d.}{\sim} F_\beta, \, \beta \in \mathcal{B},  \\
     H_0: \beta \in \mathcal{B}_0 \mbox{ v.s. } 
     H_1: \beta \in  \mathcal{B}_1 = \mathcal{B} \cap \mathcal{B}_0^c \\
    \end{array}
\end{displaymath}

\vspace{5mm}

Test Statistic:
\begin{displaymath}
    G^2 = -2 \log \left(   
           \frac{\max_{\beta \in \mathcal{B}_0} \ell(\beta)}
                {\max_{\beta \in \mathcal{B}} \ell(\beta)}
           \right)
        = -2\log\left(\frac{\ell(\widehat{\beta}_0)}{\ell(\widehat{\beta})} \right)
\end{displaymath}
\end{frame}

\begin{frame}
\frametitle{What to do}
\framesubtitle{And how to think about it}
\begin{displaymath}
    G^2 = -2 \log \left(   
           \frac{\max_{\beta \in \mathcal{B}_0} \ell(\beta)}
                {\max_{\beta \in \mathcal{B}} \ell(\beta)}
           \right)
        = -2\log\left(\frac{\ell(\widehat{\beta}_0)}{\ell(\widehat{\beta})} \right)
\end{displaymath}

  \begin{itemize}
    \item Maximize the likelihood over the whole parameter space. You already did this to calculate the MLE, denoted by $\widehat{\beta}$. Evaluate the likelihood there. That's the denominator.
    \item Maximize the likelihood over just the parameter values where $H_0$ is true -- that is, over $\mathcal{B}_0$. This yields a restricted MLE, denoted by $\widehat{\beta}_0$. Evaluate the likelihood there. That's the numerator.
    \item The numerator cannot be larger, because $\mathcal{B}_0 \subset \mathcal{B}$.
    \item If the numerator is a \emph{lot} less than the denominator, the null hypothesis is unbelievable, and
    \begin{itemize}
        \item The ratio is close to zero
        \item The log of the ratio is a big negative number
        \item $-2$ times the log is a big positive number
        \item Reject $H_0$ when $G^2$ is large enough.
    \end{itemize}
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Distribution of $G^2$ under $H_0$}
Given some technical conditions,
  \begin{itemize}
    \item $G^2$ has an approximate chi-squared distribution under $H_0$ for large $n$.   
    \item Degrees of freedom equal number of (non-redundant) equalities specified by $H_0$.
    \item Reject $H_0$ when $G^2$ is larger than the chi-squared critical value.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Counting degrees of freedom}
%\framesubtitle{For the likelihood ratio test} 
  \begin{itemize}
    \item Express $H_0$ as a set of linear combinations of the parameters, set equal to constants (usually zeros for regression problems). 
    \item Degrees of freedom = number of non-redundant (linearly independent) linear combinations.
  \end{itemize}

\vspace{5mm}

Suppose $\boldsymbol{\beta} = (\beta_1, \ldots \beta_7)$, with 

$H_0: \beta_1=\beta_2$, $\beta_6=\beta_7, \frac{1}{3}\left(\beta_1+\beta_2+\beta_3\right) = 
\frac{1}{3}\left(\beta_4+\beta_5+\beta_6\right) $

Then $df=3$: Count the equals signs.

\vspace{5mm}

But if $\beta_1=\beta_2=\beta_3=\beta_4=\beta_5=\beta_6=\beta_7=\frac{1}{7}$, then $df=6$

\end{frame}

\begin{frame}
\frametitle{Example}
%\framesubtitle{} 
University administrators recognize that the percentage of students who are unemployed after graduation will vary depending upon economic conditions, but they claim that still, about twice as many students will be employed in a job related to their field of study, compared to those who get an unrelated job. To test this hypothesis, they select a random sample of 200 students from the most recent class, and observe 106 employed in a job related to their field of study, 74 employed in a job unrelated to their field of study, and 20 unemployed. Test the hypothesis using a large-sample likelihood ratio test and the usual 0.05 significance level State your conclusions in symbols and words.
\end{frame}


\begin{frame}
\frametitle{Some detailed questions}
\framesubtitle{To guide us through the problem} 
  \begin{itemize}
    \item What is the model? 
\begin{displaymath}
    Y_1, \ldots, Y_n  \stackrel{i.i.d.}{\sim} M(1,(\pi_1,\pi_2,\pi_3))
\end{displaymath}
    \item What is the null hypothesis, in symbols?
\begin{displaymath}
    H_0: \pi_1=2\pi_2
\end{displaymath}
    \item What are the degrees of freedom for this test?
\begin{displaymath}
    1
\end{displaymath}
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{What is the parameter space $\mathcal{B}$?}
\framesubtitle{What is the restricted parameter space $\mathcal{B}_0$?} 

\begin{center}
\includegraphics[width=2.5in]{ParameterSpace}
\end{center}
% R code for plot is after the end of the document.
\begin{eqnarray*}
    \mathcal{B}   & = & \{(\pi_1,\pi_2): 0<\pi_1<1,  0<\pi_2<1, \pi_1+\pi_2<1 \} \\
    \mathcal{B}_0 & = & \{(\pi_1,\pi_2): 0<\pi_1<1,  0<\pi_2<1, \pi_1 = 2\pi_2 \} 
\end{eqnarray*}
\end{frame}

\begin{frame}
\frametitle{What is the unrestricted MLE?}

Give the answer in both symbolic and numerical form. Just write it down. There is no need to show any work. 

\vspace{10mm}

\begin{eqnarray*}
    \mathbf{p} & = & (\frac{n_1}{n},\frac{n_2}{n}, \frac{n_3}{n}) \\
      &   & \\
& = & (\frac{106}{200},\frac{74}{200}, \frac{20}{200}) \\
      &   & \\
      & = & (0.53,0.37,0.10) 
\end{eqnarray*}
\end{frame}


\begin{frame}
\frametitle{Derive the restricted MLE}

Your answer is a symbolic expression. It's a vector. Show your work.

\begin{eqnarray*}
& & \frac{\partial}{\partial \pi} \left(n_1 \log (2\pi) + n_2 \log \pi + n_3 \log(1-3\pi)\right) \\
 & = & \frac{n_1}{\pi} + \frac{n_2}{\pi} + \frac{n_3}{1-3\pi}(-3)
 \stackrel{\mbox{set}}{=} 0\\
 & \Rightarrow & \frac{n_1+n_2}{\pi} = \frac{3n_3}{1-3\pi} \\
 & \Rightarrow & (n_1+n_2)(1-3\pi) = 3\pi n_3 \\
 & \Rightarrow & n_1+n_2 = 3\pi(n_1+n_2+n_3) = 3\pi n \\
 & \Rightarrow & \pi = \frac{n_1+n_2}{3n}
\end{eqnarray*}
So $\widehat{\boldsymbol{\pi}} = \left( \frac{2(n_1+n_2)}{3n}, \frac{n_1+n_2}{3n}, \frac{n_3}{n}  \right)$. From now on, $\widehat{\boldsymbol{\pi}}$ means the \emph{restricted} MLE.
% I'm using pi-hat for the restricted MLE. I think this is going to work because unrestricted MLE is always p.
\end{frame}

\begin{frame}
\frametitle{Give the restricted MLE in numeric form}
\framesubtitle{The answer is a vector of 3 numbers} 

\begin{eqnarray*}
    \widehat{\boldsymbol{\pi}} & = & 
    \left( \frac{2(n_1+n_2)}{3n}, \frac{n_1+n_2}{3n}, \frac{n_3}{n}  \right) \\
    & = & \left(\frac{2(106+74)}{600},\frac{106+74}{600}, \frac{20}{200} \right) \\
    & = & (0.6, 0.3, 0.1)
\end{eqnarray*}
\end{frame}

\begin{frame}
\frametitle{Show the restricted and unrestricted MLEs}
\framesubtitle{Restricted is black diamond, unrestricted is red circle} 
\begin{center}
\includegraphics[width=3in]{ParameterSpaceWithMLEs}
\end{center}

\end{frame}




\begin{frame}
\frametitle{Calculate $G^2$. Show your work.}
\framesubtitle{The answer is a number.}
\begin{eqnarray*}
G^2 & = & -2\log \frac{\widehat{\pi}_1^{n_1} \widehat{\pi}_2^{n_2} p_3^{n_3}}
          {p_1^{n_1} p_2^{n_2} p_3^{n_3}} \\
    & = & -2\left(\log\left[\frac{\widehat{\pi}_1}{p_1}\right]^{n_1} + 
                  \log\left[\frac{\widehat{\pi}_2}{p_2}\right]^{n_2}  \right) \\
    & = & -2\left(n_1\log\frac{\widehat{\pi}_1}{p_1} + 
                  n_2\log\frac{\widehat{\pi}_2}{p_2}  \right) \\
    & = & -2\left(106\log\frac{0.60}{0.53} + 74\log\frac{0.30}{0.37}  \right) \\
    & = & 4.739
\end{eqnarray*}
\end{frame}


\begin{frame}[fragile]
\frametitle{Do the calculation with R.}
\framesubtitle{Display the critical value and $p$-value as well.} 
\begin{verbatim}
> G2 <- -2*(106*log(60/53)+74*log(30/37)); G2
[1] 4.739477
> qchisq(0.95,df=1) # Critical value
[1] 3.841459
> pval <- 1-pchisq(G2,df=1); pval
[1] 0.02947803
\end{verbatim}
\end{frame}


\begin{frame}
\frametitle{State your conclusions}
  \begin{itemize}
    \item \emph{In symbols}: Reject $H_0: \pi_1=2\pi_2$ at $\alpha = 0.05$.
    \item \emph{In words}: More graduates appear to be employed in jobs unrelated to their fields of study than predicted.
  \end{itemize}

The statement in words can be justified by comparing observed frequencies to those expected under $H_0$.

\vspace{2mm}

\begin{center}
\begin{tabular}{rrrr}
           & Related & Unrelated & Unemployed \\ \hline
Observed   & 106 &  74 &  20 \\
Expected   & 120 &  60 &  20 \\
Residual & -14 &  14 &   0 \\ \hline
\end{tabular}
\end{center}

\vspace{5mm}

Expected frequency is $E(n_j) = \mu_j = n\pi_j$. Estimated expected frequency is 
$\widehat{\mu}_j = n\widehat{\pi}_j$
\end{frame}


\begin{frame}
\frametitle{Write $G^2$ in terms of observed and expected frequencies}
\framesubtitle{For a general hypothesis about a multinomial} 

\begin{eqnarray*}
G^2 & = & -2 \log \left(   \frac{\ell_0}{\ell_1}  \right) \\
    & = & -2 \log \left( 
            \frac{\prod_{j=1}^c \widehat{\pi}_j^{n_j}}
                 {\prod_{j=1}^c p_j^{n_j}}
            \right) \\
    & = & -2 \log \prod_{j=1}^k \left(\frac{\widehat{\pi}_j}{p_j}\right)^{n_j} 
      =   2 \sum_{j=1}^c -\log \left(\frac{\widehat{\pi}_j}{p_j}\right)^{x_j} \\
    & = & 2 \sum_{j=1}^c n_j\log \left(\frac{\widehat{\pi}_j}{p_j}\right)^{-1}             
      =   2 \sum_{j=1}^c n_j\log \left(\frac{p_j}{\widehat{\pi}_j}\right)  \\           
    & = & 2 \sum_{j=1}^c n_j\log \left(\frac{n_j}{n\widehat{\pi}_j}\right)    
      =   2 \sum_{j=1}^c n_j\log \left(\frac{n_j}{\widehat{\mu}_j}\right)
\end{eqnarray*}
\end{frame}



\begin{frame}[fragile]
\frametitle{Likelihood ratio test for the multinomial}
\framesubtitle{Jobs data}
\begin{displaymath}
    G^2 = 2 \sum_{j=1}^c n_j\log \left(\frac{n_j}{n\widehat{\pi}_j}\right) 
        = 2 \sum_{j=1}^c n_j\log \left(\frac{n_j}{\widehat{\mu}_j}\right)
\end{displaymath}

\begin{verbatim}



> freq = c(106,74,20); n = sum(freq)
> pihat = c(0.6,0.3,0.1); muhat = n*pihat
> G2 = 2 * sum(freq*log(freq/muhat)); G2
[1] 4.739477
\end{verbatim}
\end{frame}


\begin{frame}
\frametitle{Pearson's chi-squared test}
\framesubtitle{Comparing observed and expected frequencies} 

{\huge
\begin{displaymath}
    X^2 = \sum_{j=1}^c \frac{(n_j-\widehat{\mu}_j)^2}{\widehat{\mu}_j}
\end{displaymath}
}
where $\widehat{\mu}_j = n \widehat{\pi}_j$

  \begin{itemize}
    \item A large value means the observed frequencies are far from what is expected given $H_0$.
    \item A large value makes $H_0$ less believable.
    \item Distributed approximately as chi-squared for large $n$ if $H_0$ is true.
  \end{itemize}


\end{frame}

\begin{frame}
\frametitle{Pearson Chi-squared on the jobs data}

\begin{center}
\begin{tabular}{rrrr}
Observed   & 106 &  74 &  20 \\
Expected   & 120 &  60 &  20 \\
\end{tabular}
\end{center}

\begin{eqnarray*}
X^2 & = & \sum_{j=1}^c \frac{(n_j-\widehat{\mu}_j)^2}{\widehat{\mu}_j} \\ \\
    & = & \frac{(106-120)^2}{120} + \frac{(74-60)^2}{60} + 0 \\ \\
    & = & 4.9 \mbox{~~~~(Compare } G^2 = 4.74)
\end{eqnarray*}
\end{frame}


\begin{frame}
\frametitle{Two chi-squared test statistics}
\framesubtitle{There are plenty more.} 

\begin{displaymath}
    G^2 = 2 \sum_{j=1}^c n_j\log \left(\frac{n_j}{\widehat{\mu}_j}\right)
    ~~~~~~~~~~
    X^2 = \sum_{j=1}^c \frac{(n_j-\widehat{\mu}_j)^2}{\widehat{\mu}_j}
\end{displaymath}


  \begin{itemize}
    \item Both compare observed to expected frequencies.
    \item By expected we mean \emph{estimated} expected: $\widehat{\mu}_j = n \widehat{\pi}_j$.
    \item Both equal zero when all observed frequencies equal the corresponding expected frequencies.
    \item Both have approximate chi-squared distributions with the same $df$ when $H_0$ is true, for large $n$. 
    \item Values are close for large $n$ when $H_0$ is true.
    \item Both go to infinity when $H_0$ is false.
    \item $X^2$ works better for smaller samples.
    \item $X^2$ is specific to multinomial data; $G^2$ is more general.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Rules of thumb}
  \begin{itemize}
    \item Small expected frequencies can create trouble by inflating the test statistic.
    \item $G^2$ is okay if all (estimated) expected frequencies are at least 5.
    \item $X^2$ is okay if all (estimated) expected frequencies are at least 1.
  \end{itemize}
\end{frame}


\begin{frame}
\frametitle{One more example: Is a die fair?}

Roll the die 300 times and observe these frequencies:
\begin{center}
\begin{tabular}{cccccc}
 1  &  2  &  3  &  4  &  5  &  6 \\ \hline
72  & 39  & 54  & 44  & 44  & 47 \\ \hline
\end{tabular}
\end{center}

{\small
  \begin{itemize}
    \item State a reasonable model for these data.
    \item Without any derivation, estimate the probability of rolling a $1$. Your answer is a number. 
    \item Give an approximate 95\% confidence interval for the probability of rolling a $1$. Your answer is a set of two numbers. 
    \item What is the null hypothesis corresponding to the \emph{main question}, in symbols?
    \item What is the parameter space $\mathcal{B}$?
    \item What is the restricted parameter space $\mathcal{B}_0$?
    \item What are the degrees of freedom? The answer is a number.
    \item What is the critical value of the test statistic at $\alpha=0.05$?  The answer is a number.
  \end{itemize}
} % End small
\end{frame}
% # Testing fairness of a die using R
% set.seed(999)
% n = 300
% mu = n * (1+numeric(6))/6; mu
% critval = qchisq(0.95,5); critval # 11.07050
% x = as.numeric( rmultinom(1,n,c(3,2,2,2,2,2)) ); x 
% # Got x = 72 39 54 44 44 47
% p = x/n; p 
% # Got p = 0.2400000 0.1300000 0.1800000 0.1466667 0.1466667 0.1566667
% # CI for pi1
% me = 1.96 * sqrt(p[1]*(1-p[1])/n); me
% p[1]-me  # 0.191671
% p[1]+me  # 0.288329
% # 1/6 = 0.1666667 is not in there
% X2 = sum((x-mu)^2/mu);   X2 # 14.04
% 1-pchisq(X2,5) # 0.01535732
% G2 = 2*sum(x*log(x/mu)); G2 # 13.12545
% 1-pchisq(G2,5) # 0.02223107


\begin{frame}
\frametitle{Questions continued}
%\framesubtitle{} 
  \begin{itemize}
    \item What are the expected frequencies under $H_0$? Give 6 numbers.
    \item Carry out the likelihood ratio test.
        \begin{itemize}
            \item What is the value of the test statistic? Your answer is a number. Show some work.
            \item Do you reject $H_0$ at $\alpha=0.05$? Answer Yes or No.
            \item Using R, calculate the $p$-value.
            \item Do the data provide convincing evidence against the null hypothesis?
        \end{itemize}
    \item Carry out Pearson test. Answer the same questions you did for the likelihood ratio test.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{More questions}
\framesubtitle{To help with the plain language conclusion} 
  \begin{itemize}
    \item Does the confidence interval for $\pi_1$ allow you to reject $H_0: \pi_1=\frac{1}{6}$ at $\alpha = 0.05$? Answer Yes or No.
    \item In plain language, what do you conclude from the test corresponding to the confidence interval? (You need not actually carry out the test.)
    \item Is there evidence that the chances of getting 2 through 6 are unequal? This question requires its own slide.
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Is there evidence that the chances of getting 2 through 6 are unequal?}
  \begin{itemize}
    \item What is the null hypothesis?
    \item What is the restricted parameter space $\mathcal{B}_0$? It's convenient to make the first category the residual category.
    \item Write the likelihood function for the restricted model. How many free parameters are there in this model?
    \item Obtain the restricted MLE $\widehat{\boldsymbol{\pi}}$. Your final answer is a set of 6 numbers. % 0.240 0.152 0.152 0.152 0.152 0.152
    \item Give the estimated expected frequencies $(\widehat{\mu}_1, \ldots, \widehat{\mu}_6)$. % mu = c(72,45.6,45.6,45.6,45.6,45.6)
    \item Calculate the likelihood ratio test statistic. Your answer is a number. % 2.797361
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{Questions continued}
  \begin{itemize}
    \item What are the degrees of freedom of the test? The answer is a number.
    \item What is the critical value of the test statistic at $\alpha=0.05$?  The answer is a number. % qchisq(0.95,4) = 9.487729
    \item Do you reject $H_0$ at $\alpha = 0.05$? Answer Yes or No.
    \item In plain language, what (if anything) do you conclude from the test.

\vspace{3mm}

    \item In plain language, what are your overall conclusion about this die? 
  \end{itemize}

For most statistical analyses, your final conclusions should be regarded as hypotheses that need to be tested on a new set of data.
\end{frame}

\begin{comment}
\section{Power}

\begin{frame}
\frametitle{Power and sample size}
\framesubtitle{Using the non-central chi-squared distribution} 

If $X \sim N(\mu,\sigma^2)$, then
  \begin{itemize}
    \item $Z = \left(\frac{X-\mu}{\sigma}\right)^2 \sim \chi^2(1)$
    \item $Y = \frac{X^2}{\sigma^2}$ is said to have a \emph{non-central chi-squared} distribution with degrees of freedom one and \emph{non-centrality parameter} $\lambda = \frac{\mu^2}{\sigma^2}$.
    \item Write $Y \sim \chi^2(1,\lambda)$
  \end{itemize}
\end{frame}


\begin{frame}
\frametitle{Facts about the non-central chi-squared distribution}
\framesubtitle{With one $df$}
{\huge
 \begin{displaymath}
    Y \sim \chi^2(1,\lambda), \mbox{ where } \lambda \geq 0
\end{displaymath}
}
  \begin{itemize}   
    \item $Pr\{Y>0\}=1$, of course.
    \item If $\lambda=0$, the non-central chi-squared reduces to the ordinary central chi-squared.
    \item The distribution is ``stochastically increasing" in $\lambda$, meaning that if $Y_1 \sim \chi^2(1,\lambda_1)$ and $Y_2 \sim \chi^2(1,\lambda_2)$ with $\lambda_1>\lambda_2$, then $Pr\{Y_1>y\}$ $>$ $Pr\{Y_2>y\}$ for any $y>0$. 
    \item $\lim_{\lambda \rightarrow \infty} Pr\{Y>y\} = 1$
    \item There are efficient algorithms for calculating non-central chi-squared probabilities. R's \texttt{pchisq} function does it. 
  \end{itemize}
\end{frame}

\begin{frame}
\frametitle{An example}
\framesubtitle{Back to the coffee taste test} 

$Y_1, \ldots, Y_n \stackrel{i.i.d.}{\sim} B(1,\pi)$ 

\vspace{3mm}

$H_0:\pi=\pi_0=\frac{1}{2}$

%\vspace{3mm}

Reject $H_0$ if $|Z_2| = \left|\frac{\sqrt{n}(p-\pi_0)}{\sqrt{p(1-p)}} \right| > z_{\alpha/2}$

\vspace{10mm}

Suppose that in the population, 60\% of consumers would prefer the new blend. If we test 100 consumers, what is the probability of obtaining results that are statistically significant?

\vspace{3mm}

That is, if $\pi=0.60$, what is the power with $n=100$?
\end{frame}

\begin{frame}
%\frametitle{Non-central chi-square reasoning}
\frametitle{Recall that if $X \sim N(\mu,\sigma^2)$, then $\frac{X^2}{\sigma^2} \sim \chi^2(1,\frac{\mu^2}{\sigma^2})$.} 

Reject $H_0$ if
\begin{displaymath}
    |Z_2| = \left|\frac{\sqrt{n}(p-\pi_0)}{\sqrt{p(1-p)}} \right| > z_{\alpha/2} \Leftrightarrow  Z_2^2 > z_{\alpha/2}^2 = \chi^2_{\alpha}(1) 
\end{displaymath}


For large $n$, $X = p-\pi_0$ is approximately normal, with $\mu=\pi-\pi_0$ and 
          $\sigma^2 = \frac{\pi(1-\pi)}{n}$. So, 

\begin{eqnarray*}
    Z_2^2 & = & \frac{(p-\pi_0)^2}{p(1-p)/n} \approx \frac{(p-\pi_0)^2}{\pi(1-\pi)/n}
                 = \frac{X^2}{\sigma^2} \\
          && \\
          & \stackrel{approx}{\sim} &  \chi^2\left(1, n\,\frac{(\pi-\pi_0)^2}{\pi(1-\pi)}\right)
\end{eqnarray*}
\end{frame}

\begin{frame}
\frametitle{We have found that}
%\framesubtitle{} 
The Wald chi-squared test statistic of $H_0:\pi=\pi_0$
\begin{displaymath}
    Z_2^2 = \frac{n(p-\pi_0)^2}{p(1-p)}
\end{displaymath}
has an asymptotic non-central chi-squared distribution with $df=1$ and non-centrality parameter
\begin{displaymath}
    \lambda = n\,\frac{(\pi-\pi_0)^2}{\pi(1-\pi)}
\end{displaymath}
Notice the similarity, and also that
  \begin{itemize}
    \item If $\pi=\pi_0$, then $\lambda=0$ and $Z_2^2$ has a central chi-squared distribution.
    \item The probability of exceeding any critical value (power) can be made as large as desired by making $\lambda$ bigger.
    \item There are 2 ways to make $\lambda$ bigger.
  \end{itemize}
\end{frame}


\begin{frame}[fragile]
\frametitle{Power calculation with R}
\framesubtitle{For $n=100$, $\pi_0=0.50$ and $\pi=0.60$} 

\begin{verbatim}
> # Power for Wald chisquare test of H0: pi=pi0
> n=100; pi0=0.50; pi=0.60
> lambda = n * (pi-pi0)^2 / (pi*(1-pi))
> critval = qchisq(0.95,1)
> power = 1-pchisq(critval,1,lambda); power
[1] 0.5324209
\end{verbatim}
\end{frame}


\begin{frame}
\frametitle{General non-central chi-squared}
%\framesubtitle{Definition}

Let $X_1, \ldots, X_n$ be independent $N(\mu_i,\sigma^2_i)$. Then 
\begin{displaymath}
    Y = \sum_{i=1}^n \frac{X_i^2}{\sigma^2_i} \sim \chi^2(n,\lambda),
    \mbox{ where } \lambda = \sum_{i=1}^n \frac{\mu_i^2}{\sigma^2_i}
\end{displaymath}

  \begin{itemize}
    \item Density is a bit messy
    \item Reduces to central chi-squared when $\lambda=0$.
    \item Generalizes to $Y \sim \chi^2(\nu,\lambda)$, where $\nu>0$ as well as $\lambda>0$
    \item Stochastically increasing in $\lambda$, meaning $Pr\{Y>y\}$ can be increased by increasing $\lambda$.
    \item $\lim_{\lambda \rightarrow \infty} Pr\{Y>y\} = 1$
    \item Probabilities are easy to calculate numerically.    
  \end{itemize}
\end{frame}

\begin{frame} 
\frametitle{Non-centrality parameters for Pearson and likelihood Ratio chi-squared tests}
\framesubtitle{Re-writing the test statistics in terms of $p_j$ \ldots} 

\begin{center}
\begin{tabular}{ll}
$X^2 = n\displaystyle{\sum_{j=1}^c 
            \frac{\left(p_j-\widehat{\pi}_j\right)^2}{\widehat{\pi}_j}   }$ & 
$\lambda = n\displaystyle{\sum_{j=1}^c 
            \frac{\left[\pi_j-\pi_j(M)\right]^2}{\pi_j(M)}   }$ \\
          &  \\
$G^2 = 2n\displaystyle{\sum_{j=1}^c 
            p_j \log\left(\frac{p_j}{\widehat{\pi}_j} \right)   }$ & 
$\lambda = 2n\displaystyle{\sum_{j=1}^c 
            \pi_j \log\left(\frac{\pi_j}{\pi_j(M)} \right)   }$, \\
& \\
\end{tabular}
\end{center}
% Agresti CDA p. 241
{ \small
  \begin{itemize}
    \item Where $\widehat{\pi}_j \rightarrow \pi_j(M)$ as $n
      \rightarrow \infty$ under the \emph{alternative} hypothesis.

    \item That is, $\pi_j(M)$ is the large-sample target of the
      restricted MLE.

    \item The technical meaning is convergence in probability, or $
      \widehat{\pi}_j \stackrel{p}{\rightarrow} \pi_j(M)$

    \item $\pi_j(M)$ is a function of $\pi_1, \ldots,
      \pi_c$. 
% \emph{Sometimes}, $\pi_j(M)=\pi_j$

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

\begin{frame}
\frametitle{For the fair (?) die example}
Suppose we want to test whether the die is fair, but it is not.  The true $\boldsymbol{\pi} = \left(\frac{2}{13}, \frac{2}{13}, \frac{3}{13}, \frac{2}{13}, \frac{2}{13}, \frac{2}{13} \right)$. What is the power of the Pearson chi-squared test for $n=300$?

\vspace{3mm}

Because $\mathcal{B}_0 = \{(\frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6})\}$, It's easy to see $\pi_j(M)=\frac{1}{6}$.

\begin{eqnarray*}
    \lambda & = & n \sum_{j=1}^c 
            \frac{\left[\pi_j-\pi_j(M)\right]^2}{\pi_j(M)} \\
         & = & 300 \left(\frac{(\frac{3}{13}-\frac{1}{6})^2}{\frac{1}{6}} + 5\frac{(\frac{2}{13}-\frac{1}{6})^2}{\frac{1}{6}} \right) \\
         & = & 8.87574
\end{eqnarray*}
\end{frame}

\begin{frame}[fragile]
\frametitle{Calculate power with R}

\begin{verbatim}
> piM = 1/6 + numeric(6); piM
[1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
> pi = c(2,2,3,2,2,2)/13; pi
[1] 0.1538462 0.1538462 0.2307692 0.1538462 0.1538462 0.1538462
> lambda = 300 * sum( (pi-piM)^2/piM ); lambda
[1] 8.87574
> critval = qchisq(0.95,5); critval
[1] 11.0705
> power = 1-pchisq(critval,5,lambda); power
[1] 0.6165159
\end{verbatim}
\end{frame}


\begin{frame}[fragile]
\frametitle{Calculate required sample size}
\framesubtitle{To detect $\boldsymbol{\pi} = \left(\frac{2}{13}, \frac{2}{13}, \frac{3}{13}, \frac{2}{13}, \frac{2}{13}, \frac{2}{13} \right)$ as unfair with high probability}
Power of 0.62 is not too impressive. How many times would you have to roll the die to detect this degree of unfairness with probability 0.90?

{\small
\begin{verbatim}
> piM = 1/6 + numeric(6)
> pi = c(2,2,3,2,2,2)/13
> critval = qchisq(0.95,5)
> 
> n = 0; power = 0.05
> while(power < 0.90)
+     {   n = n+1
+         lambda = n * sum( (pi-piM)^2/piM )
+         power = power = 1-pchisq(critval,5,lambda)
+     }
> n; power
[1] 557
[1] 0.9001972
\end{verbatim}
}
% This is just approximate: Got 0.8702, 0.8727, 0.8706 by simulation.
\end{frame}


\begin{frame}
\frametitle{Sometimes, finding $\boldsymbol{\pi}(M)$ is more challenging}
\framesubtitle{Recall the Jobs example, with $H_0: \pi_1 = 2\pi_2$} 

\begin{eqnarray*}
    \widehat{\boldsymbol{\pi}} & = & \left( \frac{2(n_1+n_2)}{3n}, \frac{n_1+n_2}{3n}, \frac{n_3}{n}  \right) \\
&& \\
    & = & \left( \frac{2}{3}\left(\frac{n_1}{n}+\frac{n_2}{n} \right),
                 \frac{1}{3}\left(\frac{n_1}{n}+\frac{n_2}{n} \right),
                 \frac{n_3}{n}  \right) \\
&& \\
    & = &            \left( \frac{2}{3}\left(p_1+p_2 \right),
                 \frac{1}{3}\left(p_1+p_2 \right),
                 p_3  \right)  \\
&& \\
    & \stackrel{p}{\rightarrow} &
          \left( \frac{2}{3}\left(\pi_1+\pi_2 \right),
                 \frac{1}{3}\left(\pi_1+\pi_2 \right),
                 \pi_3  \right) = \boldsymbol{\pi}(M)
\end{eqnarray*}
By the Law of Large Numbers.
% In general, you can always do this kind of thing when you have an explicit formula for the restricted MLE, because it will be possible to write it in terms of the sample proportions, which converge to the true probabilities.
\end{frame}

\end{comment}

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

%\section{Copyright}

\begin{frame}
\frametitle{Copyright Information}

This document 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/312f22} {\texttt{http://www.utstat.toronto.edu/brunner/oldclass/312f22}}

\end{frame}


\end{document}

# Plotting jobs parameter space with R
pi1 = seq(from=0,to=1,by=0.05); pi2=pi1
plot(pi1,pi2,pch=' ', frame.plot=F,
     xlab=expression(pi[1]), ylab=expression(pi[2])) 
# Draw boundaries of parameter space
xloc1 = c(0,0); yloc1 = c(0,1); lines(xloc1,yloc1,lty=1)
xloc2 = c(0,1); yloc2 = c(0,0); lines(xloc2,yloc2,lty=1)
xloc3 = c(0,1); yloc3 = c(1,0); lines(xloc3,yloc3,lty=1)
# Restricted parameter space is a line segment
xloc4 = c(0,2/3); yloc4 = c(0,1/3); lines(xloc4,yloc4,lty=2)
points(0.53,0.37, pch=19, col = "red1") # Unrestricted MLE
points(0.60,0.30, pch=23, bg="black") # Restricted MLE


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

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



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