\documentclass[12pt]{article}

\usepackage{Math203}

\usepackage{amsmath,amssymb}
\usepackage{hyperref}
\usepackage{color}
\usepackage{xcolor}
\usepackage{enumitem}
\usepackage{xcolor}
\usepackage{longtable}
\usepackage{booktabs}


\pagestyle{plain}

\renewcommand{\labelenumi}{\arabic{enumi}.}
\renewcommand{\labelenumii}{\arabic{enumi}.\arabic{enumii}}
\renewcommand{\labelenumiii}{\arabic{enumi}.\arabic{enumii}.\arabic{enumiii}}

\setlength{\topmargin}{-0.75 in} \setlength{\textheight}{9.5 in}

\parindent0in

\newenvironment{packed_enum}{
\begin{enumerate}
  \setlength{\itemsep}{12pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
  \setlength{\topsep}{0pt}
}{\end{enumerate}}

\renewcommand{\labelitemii}{$\blacktriangleright$}

\def\ctitle{Math 203: Principles of Statistics 1}
\def\ttitle{Fall 2024}
\def\stitle{Hypothesis testing with unequal variances}

\makeatletter
\def\@maketitle{%
  \newpage
  \null
  \begin{center}%
    {\textsclarge{\ctitle} \par}%
    \medskip
    {\textsclarge{\ttitle} \par}%
    \medskip
    {\textsc{\stitle} \par}
  \end{center}%
  \par
  \vskip 1.5em}
\makeatother

\begin{document}

<<funcs,echo=FALSE>>=
myrnd <- function(x,d) {formatC(x, format="f", digits=d)}
library(shape)
library(diagram)
@

\maketitle

<<chunk0, echo=FALSE,fig.height=5.0,fig.width=8.75>>=
N<-10000
set.seed(2453)
m<-9
n<-12
mu1<-mu2<-10
sig1<-sig2<-2
@

To test the null hypothesis $H_0 : \mu_1 = \mu_2$ against a selected alternative 

\begin{itemize}
  \item when $m$ or $n$ is \textbf{small}, we must assume \textbf{Normality} of the two populations, and that the two population \textbf{standard deviations} are the \textbf{same}, that is $\sigma_1 = \sigma_2 = \sigma$ say.  In this case, we use the test statistic
  \begin{equation}\label{eq:T1}
T = \dfrac{(\Xbar - \Ybar)}{S_P \sqrt{\dfrac{1}{m} + \dfrac{1}{n}}}
\end{equation}

and compare this with the t-distribution with $m+n-2$ degrees of freedom, where $m$ and $n$ are the two sample sizes.

  \item  when $m$ and $n$ are \textbf{large},  we use the test statistic
\begin{equation}\label{eq:T2}
T = \dfrac{(\Xbar - \Ybar)}{\sqrt{\dfrac{S_1^2}{m} + \dfrac{S_2^2}{n}}}
\end{equation}
and compare this with the standard Normal distribution.  We \textbf{do not} need to assume Normality of the original data, nor assume that the population variances are equal.

\end{itemize}
In the first simulation, $m=\Sexpr{m}$ and $n=\Sexpr{n}$, in $N=\Sexpr{myrnd(N,0)}$ replicated studies with $\mu_1 = \mu_2 = \Sexpr{mu1}$ and $\sigma_1 = \sigma_2 = \Sexpr{sig1}$, with the data Normally distributed.

<<chunk1, echo=FALSE,fig.height=4.8,fig.width=8.75>>=
xmat<-matrix(NA,nrow=N, ncol=m)
ymat<-matrix(NA,nrow=N, ncol=n)
for(i in 1:N){
    xmat[i,]<-rnorm(m,mu1,sig1)
    ymat[i,]<-rnorm(n,mu2,sig2)
}
xmeans<-apply(xmat,1,mean)
ymeans<-apply(ymat,1,mean)
xsds<-apply(xmat,1,sd)
ysds<-apply(ymat,1,sd)
sps<-((m-1)*xsds^2+(n-1)*ysds^2)/(m+n-2)
t.obs.eq<-(xmeans-ymeans)/sqrt(sps*(1/m+1/n))
t.obs.uneq<-(xmeans-ymeans)/sqrt(xsds^2/m+ysds^2/n)
par(mfrow=c(1,2),mar=c(4,4,2,2))
hist(t.obs.eq,nclass=100,xlim=range(-5,5),ylim=range(0,0.5),main='Statistic (1)',xlab=expression(T[obs]),freq=FALSE);box()
xv<-seq(-5,5,by=0.001)
yv1<-dt(xv,m+n-2)
lines(xv,yv1,col='red')
hist(t.obs.uneq,nclass=100,xlim=range(-5,5),ylim=range(0,0.5),main='Statistic (2)',xlab=expression(T[obs]),freq=FALSE);box()
yv2<-dnorm(xv)
lines(xv,yv2,col='red')
@

The distribution of the two statistics follow the presumed exact and approximate null distributions plotted in red: $t_{\Sexpr{m+n-2}}$ for statistic (1), and standard Normal for statistic (2).  


\medskip

Now we repeat the calculation with $\sigma_1 = 2$ but $\sigma_2 = 5$, still with the data Normally distributed.

<<chunk2, echo=FALSE,fig.height=4.8,fig.width=8.75>>=
mu1<-mu2<-10
sig1<-2; sig2<-5
xmat<-matrix(NA,nrow=N, ncol=m)
ymat<-matrix(NA,nrow=N, ncol=n)
for(i in 1:N){
    xmat[i,]<-rnorm(m,mu1,sig1)
    ymat[i,]<-rnorm(n,mu2,sig2)
}
xmeans<-apply(xmat,1,mean)
ymeans<-apply(ymat,1,mean)
xsds<-apply(xmat,1,sd)
ysds<-apply(ymat,1,sd)
sps<-((m-1)*xsds^2+(n-1)*ysds^2)/(m+n-2)
t.obs.eq<-(xmeans-ymeans)/sqrt(sps*(1/m+1/n))
t.obs.uneq<-(xmeans-ymeans)/sqrt(xsds^2/m+ysds^2/n)
par(mfrow=c(1,2),mar=c(4,4,2,2))
hist(t.obs.eq,nclass=100,xlim=range(-5,5),ylim=range(0,0.5),main='Statistic (1)',xlab=expression(T[obs]),freq=FALSE);box()
xv<-seq(-5,5,by=0.001)
yv1<-dt(xv,m+n-2)
lines(xv,yv1,col='red')
hist(t.obs.uneq,nclass=100,xlim=range(-5,5),ylim=range(0,0.5),main='Statistic (2)',xlab=expression(T[obs]),freq=FALSE);box()
yv2<-dnorm(xv)
lines(xv,yv2,col='red')
nu.val<-(sig1^2/m+sig2^2/n)^2/((sig1^2/m)^2/(m-1)+(sig2^2/n)^2/(n-1))
@

The null distribution of statistic (1) is no longer the t-distribution computed previously, as expected. The distribution of statistic (2) is still apparently well-approximated by the standard Normal.  We can check this with a QQ plot: 

<<chunk3, echo=FALSE,fig.height=5,fig.width=8>>=
par(pty='s',mar=c(4,4,2,2))
qqnorm(t.obs.uneq,pch=19,cex=0.7)
abline(0,1)
@

The QQ plot reveals that in fact the null distribution of statistic (2) is heavier tailed than the standard Normal.  This is important, because when we carry out hypothesis tests, we are most interested in the tails of the null distribution.  A standard approximation to the null distribution in this case is the t-distribution with
\begin{equation}\label{eq:nuval}
\nu = \dfrac{\left(\dfrac{s_1^2}{m} + \dfrac{s_2^2}{n}\right)^2 }{\dfrac{1}{m-1} \left(\dfrac{s_1^2}{m}\right)^2 + \dfrac{1}{n-1} \left(\dfrac{s_2^2}{n}\right)^2  }
\end{equation}
Note that this degrees of freedom depends on $m$ and $n$, but also the two observed sample variances (and therefore the on true values of $\sigma_1$ and $\sigma_2$) and is not in general a whole number. The maximum value $\nu$ can take is $m+n-2$.

\medskip
We can compute the value of $\nu$ for each of the previous replicates:  in this case, the maximum value $\Sexpr{m+n-2}$, and the value we obtain by using $\sigma_1^2$ and $\sigma_2^2$ into \eqref{eq:nuval} is \Sexpr{myrnd(nu.val,4)}.

<<chunk4, echo=FALSE,fig.height=6,fig.width=8.75>>=
nu.samp<-(xsds^2/m+ysds^2/n)^2/((xsds^2/m)^2/(m-1)+(ysds^2/n)^2/(n-1))
hist(nu.samp,br=seq(9,20,by=0.1),xlab=expression(nu),main=expression(paste('Histogram of ',nu,' (m=9,n=12)')));box()
abline(v=nu.val,lty=3,col='red')
@

As $\nu$ varies across replicates, we check the uniformity of p-values derived in a test against a two-sided alternative.  That is, for each we compute
\[
P_0(|T| \geq T_{obs})
\]
with $P_0$ being the t-distribution with $\nu$ degrees of freedom, with $\nu$ computed using \eqref{eq:nuval}.  This probability should be have a uniform distribution if the null hypothesis is true.

<<chunk5, echo=FALSE,fig.height=5,fig.width=8.75>>=
pval<-2*(1-pt(abs(t.obs.uneq),nu.samp))
par(mar=c(4,4,2,2))
hist(pval,br=seq(0,1,by=0.02),freq=FALSE,
    xlab='p-value',main=expression('Histogram of p-values'));box()
abline(h=1,lty=2)
@

This confirms that the approximation is a good one at least for these values of $m$ and $n$.  We can therefore use this approximation to carry out the test, and also to compute confidence intervals.

\medskip

If the sample sizes increase, the distribution of test statistic \eqref{eq:T2} becomes more like the Normal distribution
 
<<chunk6, echo=FALSE,fig.height=5.5,fig.width=8.75>>=
m<-45
n<-60
mu1<-mu2<-10
sig1<-2; sig2<-5
xmat<-matrix(NA,nrow=N, ncol=m)
ymat<-matrix(NA,nrow=N, ncol=n)
for(i in 1:N){
    xmat[i,]<-rnorm(m,mu1,sig1)
    ymat[i,]<-rnorm(n,mu2,sig2)
}
xmeans<-apply(xmat,1,mean)
ymeans<-apply(ymat,1,mean)
xsds<-apply(xmat,1,sd)
ysds<-apply(ymat,1,sd)
sps<-((m-1)*xsds^2+(n-1)*ysds^2)/(m+n-2)
t.obs.eq<-(xmeans-ymeans)/sqrt(sps*(1/m+1/n))
t.obs.uneq<-(xmeans-ymeans)/sqrt(xsds^2/m+ysds^2/n)
nu.samp<-(xsds^2/m+ysds^2/n)^2/((xsds^2/m)^2/(m-1)+(ysds^2/n)^2/(n-1))
nu.val<-(sig1^2/m+sig2^2/n)^2/((sig1^2/m)^2/(m-1)+(sig2^2/n)^2/(n-1))
pval<-2*(1-pt(abs(t.obs.uneq),nu.samp))
par(mar=c(4,4,2,2))
hist(t.obs.uneq,nclass=100,xlim=range(-5,5),ylim=range(0,0.5),main='Statistic (2) (m=45,n=60)',xlab=expression(T[obs]),freq=FALSE);box()
yv2<-dnorm(xv)
lines(xv,yv2,col='red')
@

We can verify the Normality with a QQ plot:

<<chunk7, echo=FALSE,fig.height=5.5,fig.width=8.75>>=
par(pty='s',mar=c(4,4,2,2))
qqnorm(t.obs.uneq,pch=19,cex=0.7)
abline(0,1)
@

However, the degrees of freedom computed using \eqref{eq:nuval} when $\sigma_1 = 2$ and $\sigma_2 = 5$ is still below $m+n-2=\Sexpr{m+n-2}$.  The value we obtain by using $\sigma_1^2$ and $\sigma_2^2$ in \eqref{eq:nuval} is \Sexpr{myrnd(nu.val,4)}.

<<chunk8, echo=FALSE,fig.height=5.5,fig.width=8.75>>=
par(pty='m',mar=c(4,4,2,2))
hist(nu.samp,br=seq(50,120,by=1),xlab=expression(nu),main=expression(paste('Histogram of ',nu,' (m=45,n=60)')));box()
abline(v=m+n-2,lty=3)
abline(v=nu.val,lty=3,col='red')
@

Finally, we confirm that the null distribution of statistic (2) changes as $\sigma_1$ and $\sigma_2$ change: if $\sigma_1 = 8$ and $\sigma_2 = 5$, we obtain the following computed values of $\nu$ in the \Sexpr{myrnd(N,0)} replicates:

<<chunk9, echo=FALSE,fig.height=5.5,fig.width=8.75>>=
m<-45
n<-60
mu1<-mu2<-10
sig1<-8; sig2<-5
xmat<-matrix(NA,nrow=N, ncol=m)
ymat<-matrix(NA,nrow=N, ncol=n)
for(i in 1:N){
    xmat[i,]<-rnorm(m,mu1,sig1)
    ymat[i,]<-rnorm(n,mu2,sig2)
}
xmeans<-apply(xmat,1,mean)
ymeans<-apply(ymat,1,mean)
xsds<-apply(xmat,1,sd)
ysds<-apply(ymat,1,sd)
sps<-((m-1)*xsds^2+(n-1)*ysds^2)/(m+n-2)
t.obs.eq<-(xmeans-ymeans)/sqrt(sps*(1/m+1/n))
t.obs.uneq<-(xmeans-ymeans)/sqrt(xsds^2/m+ysds^2/n)
nu.samp<-(xsds^2/m+ysds^2/n)^2/((xsds^2/m)^2/(m-1)+(ysds^2/n)^2/(n-1))
nu.val<-(sig1^2/m+sig2^2/n)^2/((sig1^2/m)^2/(m-1)+(sig2^2/n)^2/(n-1))
pval<-2*(1-pt(abs(t.obs.uneq),nu.samp))
par(pty='m',mar=c(4,4,2,2))
hist(nu.samp,br=seq(50,120,by=1),xlab=expression(nu),main=expression(paste('Histogram of ',nu,' (m=45,n=60)')));box()
abline(v=nu.val,lty=3,col='red')
@

The degrees of freedom computed using \eqref{eq:nuval} when $\sigma_1 = 8$ and $\sigma_2 = 5$ using $\sigma_1^2$ and $\sigma_2^2$ in \eqref{eq:nuval} is \Sexpr{myrnd(nu.val,4)}.


\end{document}