\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{Interpreting the Inference}

\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

In each of 100 counties, data are collected measuring the IQ (measured by standard testing) of a small number of infants in the county, and at the same time, whether that infant lives in a household with a low or high level of an environmental factor.   

\medskip

The following plot displays boxplots of the IQ samples across the 100 counties.

<<chunk1, echo=FALSE,fig.height=4.5,fig.width=8.5>>=
N<-100
set.seed(2453)
n1<-rpois(N,15)
n2<-rpois(N,15); #print(c(max(n1),max(n2),max(n1+n2)))
xmat<-ymat<-matrix(NA,nrow=N, ncol=max(c(n1,n2)))
for(i in 1:N){
    xmat[i,1:n1[i]]<-rnorm(n1[i],100,12)
    ymat[i,1:n2[i]]<-rnorm(n2[i],100,12)
}
par(mar=c(4,4,3,2))
boxplot(t(xmat),pch=19,cex=0.6,xlab='County',xaxt='n',ylim=range(60,140),main='Low level')
axis(1,at=seq(5,100,by=5))
boxplot(t(ymat),pch=19,cex=0.6,xlab='County',xaxt='n',ylim=range(60,140),main='High level')
axis(1,at=seq(5,100,by=5))
@

We can also plot the differences between the IQ means of the Low and High level households:

<<chunk2, echo=FALSE,fig.height=4.8,fig.width=8.75>>=
xmeans<-apply(xmat,1,mean,na.rm=T)
ymeans<-apply(ymat,1,mean,na.rm=T)
xsds<-apply(xmat,1,sd,na.rm=T)
ysds<-apply(ymat,1,sd,na.rm=T)
par(mar=c(4,4,1,2))
plot(xmeans-ymeans,pch=19,ylim=range(-30,30),cex=0.8,xlab='County',xaxt='n',ylab='Difference in Means')
axis(1,at=seq(5,100,by=5))
abline(h=0,lty=3,col='red')
alpha<-0.05
ci<-matrix(0,nrow=N,ncol=2)
tobs<-rej<-rej1<-pval<-rep(0,N)
for(i in 1:N){
    t.al<-qt(1-alpha/2,n1[i]+n2[i]-2)
    sp<-sqrt(((n1[i]-1)*xsds[i]^2+(n2[i]-1)*ysds[i]^2)/(n1[i]+n2[i]-2))
    ci[i,]<-(xmeans[i]-ymeans[i])+c(-1,1)*t.al*sp*sqrt(1/n1[i]+1/n2[i])  
    tobs[i]<-(xmeans[i]-ymeans[i])/(sp*sqrt(1/n1[i]+1/n2[i]))
}
rej[abs(tobs)>qt(1-alpha/2,n1+n2-2)]<-1
rej1[tobs>qt(1-alpha,n1+n2-2)]<-1
pval<-1-pt(tobs,n1+n2-2)
@

For a statistical assessment, we can construct 95\% Confidence Intervals for the differences in means in each county, using the formula (in this small sample setting, under the assumption of equal population variances)
\[
(\xbar - \ybar) \pm t_{\alpha/2,m+n-2} s_p \sqrt{\dfrac{1}{m} + \dfrac{1}{n}}
\]
where $m$ and $n$ are the sample sizes for the Low and High level households.

<<chunk3, echo=FALSE,fig.height=4.8,fig.width=8.75>>=
colvec<-rep('gray',N)
colvec[ci[,1]>0 | ci[,2] < 0]<-'red'
par(mar=c(4,4,1,2))
plot(xmeans-ymeans,type='n',pch=19,cex=0.8,xlab='County',ylim=range(-30,30),xaxt='n',ylab='CI')
axis(1,at=seq(5,100,by=5))
for(i in 1:N){
    rect(i,ci[i,1],i+0.5,ci[i,2],col=colvec[i])
}
abline(h=0,lty=3,col='red')
@

The intervals coloured in red are those that exclude zero (that is, the interval lies entirely above or entirely below the zero line).  There are seven such intervals out of 100.  We can repeat also report the test statistics for tests of the hypotheses
\begin{align*}
  H_0 & : \mu_1 = \mu_2 \\
  H_a & : \mu_1 \neq \mu_2
\end{align*}
assuming a common population variance $\sigma$ for the two populations using the test statistic
\[
T_{obs} = \dfrac{\overline{x} - \overline{y}}{s_p \sqrt{\dfrac{1}{m} + \dfrac{1}{n}}}
\]
plotting $T_{obs}$ in red for the statistics indicating rejection of $H_0$ in favour of $H_a$.  The Rejection Region in these two-sided tests take the form
\[
\{|T| \geq t_{\alpha/2,m+n-2}\} = \{T \leq -t_{\alpha/2,m+n-2}\} \cup \{T \geq t_{\alpha/2,m+n-2}\}
\]
(the union of the left-tail region and the right-tail region) and we may test using $\alpha=0.05$.

<<chunk4, echo=FALSE,fig.height=5.0,fig.width=8.75>>=
rvec<-rep('gray',N)
rvec[rej==1]<-'red'
par(mar=c(4,4,2,2))
plot(tobs,pch=19,cex=0.8,col=rvec,xlab='County',xaxt='n',ylab=expression(T[obs]))
axis(1,at=seq(5,100,by=5))
abline(h=0,lty=3,col='red')
@

We see the same seven counties identified as leading to a rejection of the null hypothesis.

\medskip

It might be considered more relevant to consider the one-sided alternative hypothesis
\begin{align*}
  H_0 & : \mu_1 = \mu_2 \\
  H_a & : \mu_1 > \mu_2
\end{align*}
for which the the Rejection Region take the form
\[
\{T \geq t_{\alpha,m+n-2}\}
\]
(the right-tail region only).

<<chunk5, echo=FALSE,fig.height=5.0,fig.width=8.75>>=
rvec<-rep('gray',N)
rvec[rej1==1]<-'blue'
par(mar=c(4,4,2,2))
plot(tobs,pch=19,cex=0.8,col=rvec,xlab='County',xaxt='n',ylab=expression(T[obs]))
axis(1,at=seq(5,100,by=5))
abline(h=0,lty=3,col='red')
@

The statistics plotted in blue indicated rejection of the null hypothesis in favour of the alternative hypothesis, again using $\alpha=0.05$.  We see that five tests lead to rejection.
\bigskip

\textbf{Note:} In fact, \textbf{these data are simulated under the null hypothesis with equal variances}:
\begin{align*}
  \textrm{Low} \quad X_1,\ldots,X_m & \sim Normal(100,12^2) \\
  \textrm{High} \quad Y_1,\ldots,Y_n & \sim Normal(100,12^2)
\end{align*}
The observed variation in the test statistic, and the rejections of $H_0$ that we see, are merely the result of random variation, and there are in fact no systematic differences between the two populations.

\bigskip

These results remind us of the logic being appealed to in the hypothesis testing framework: after having formed the hypotheses and chosen a suitable test statistic $T$, we 
\begin{itemize}
  \item predict how $T$ \textbf{should} behave if the null hypothesis is \textbf{true} by forming the null distribution;
  \item compute the observed test statistic $T_{obs}$;
  \item check whether the observed result is \textbf{surprising}, compared with our prediction: does $T_{obs}$ lie in a low probability region of the null distribution?
  \item If the result is sufficiently surprising we \textbf{reject} the null hypothesis in favour of the alternative: the degree of `surprise' we are willing to tolerate is measured by $\alpha$, as 
  \[
  \alpha = \Pr(\textrm{$T$ lies in R. R.}|\textrm{$H_0$ is true}).
  \]
  
  \item even if $H_0$ is true, we \textbf{expect} to see $T_{obs}$ in the Rejection Region 100$\alpha$\% of the time.
  \end{itemize}

The computation of p-values does not change this interpretation.  For example, for the one-sided alternative hypothesis, we compute the p-value as
\[
\textrm{p-value} = P_0(T \geq T_{obs})
\]
and we can plot the 100 p-values in a histogram.

<<chunk6, echo=FALSE,fig.height=5.0,fig.width=8.75>>=
par(mar=c(4,4,3,2))
hist(pval,br=seq(0,1,by=0.1),xlab='p-value',freq=FALSE,
    main='Histogram of p-values from 100 counties');box()
@

The p-values are \textbf{uniformly} distributed across the interval $(0,1)$.  The uniformity is easier to see if we have $N=10000$ counties:

<<chunk7, echo=FALSE,fig.height=5,fig.width=8.5>>=
N<-10000
set.seed(2453)
n1<-rpois(N,15)
n2<-rpois(N,15)
xmat<-ymat<-matrix(NA,nrow=N, ncol=max(c(n1,n2)))
for(i in 1:N){
    xmat[i,1:n1[i]]<-rnorm(n1[i],100,12)
    ymat[i,1:n2[i]]<-rnorm(n2[i],100,12)
}
xmeans<-apply(xmat,1,mean,na.rm=T)
ymeans<-apply(ymat,1,mean,na.rm=T)
xsds<-apply(xmat,1,sd,na.rm=T)
ysds<-apply(ymat,1,sd,na.rm=T)
tobs<-rep(0,N)
for(i in 1:N){
    sp<-sqrt(((n1[i]-1)*xsds[i]^2+(n2[i]-1)*ysds[i]^2)/(n1[i]+n2[i]-2))
    tobs[i]<-(xmeans[i]-ymeans[i])/(sp*sqrt(1/n1[i]+1/n2[i]))
}
pval<-1-pt(tobs,n1+n2-2)
par(mar=c(4,4,3,2))
hist(pval,br=seq(0,1,by=0.01),freq=FALSE,xlab='p-value',
    main='Histogram of p-values from 10000 counties');box()
abline(h=1,lty=2)
@

Contrast this with the case where there is in fact a systematic difference between Low and High level households: if we simulate
\begin{align*}
  \text{Low} \quad X_1,\ldots,X_m & \sim Normal(100,12^2) \\
  \text{High} \quad Y_1,\ldots,Y_n & \sim Normal(95,12^2)
\end{align*}
the p-value histogram appears as follows:

<<chunk8, echo=FALSE,fig.height=5,fig.width=8.5>>=
N<-10000
set.seed(2453)
n1<-rpois(N,15)
n2<-rpois(N,15)
xmat<-ymat<-matrix(NA,nrow=N, ncol=max(c(n1,n2)))
for(i in 1:N){
    xmat[i,1:n1[i]]<-rnorm(n1[i],100,12)
    ymat[i,1:n2[i]]<-rnorm(n2[i],95,12)
}
xmeans<-apply(xmat,1,mean,na.rm=T)
ymeans<-apply(ymat,1,mean,na.rm=T)
xsds<-apply(xmat,1,sd,na.rm=T)
ysds<-apply(ymat,1,sd,na.rm=T)
tobs<-rep(0,N)
for(i in 1:N){
    sp<-sqrt(((n1[i]-1)*xsds[i]^2+(n2[i]-1)*ysds[i]^2)/(n1[i]+n2[i]-2))
    tobs[i]<-(xmeans[i]-ymeans[i])/(sp*sqrt(1/n1[i]+1/n2[i]))
}
pval<-1-pt(tobs,n1+n2-2)
par(mar=c(4,4,3,2))
hist(pval,br=seq(0,1,by=0.01),freq=FALSE,xlab='p',main='Histogram of p-values');box()
abline(h=1,lty=2)
@

If we use a p-value threshold of $\alpha=0.05$, it turns out that \Sexpr{sum(pval<alpha)} of the tests lead to p-values below $\alpha$ and a (correct) rejection of the null hypothesis.

\end{document}