lecture2.Rnw

\documentclass[compress,red]{beamer}
%\documentclass[handout]{beamer}

%\usepackage{beamerthemesplit}
\usepackage{fancyvrb}
\usepackage{graphicx}
\usepackage{gb4e}
\usepackage{qtree}
\usepackage{hyperref}
\usepackage{ulem}

\usepackage{amsmath,amssymb,amsfonts}

\setbeamerfont{page number in head/foot}{size=\large}
\setbeamertemplate{footline}[frame number]

%\setbeamertemplate{footline}%
%{%
%\hfill\insertpagenumber\ of \ref{TotPages}\hspace{.5cm}\vspace{.5cm}
%\hfill\insertpagenumber\ of 28\hspace{.5cm}\vspace{.5cm}
%}%



\mode<presentation>
{
%\usetheme{Singapore}
%\usetheme{Berlin}

%\setbeamercovered{transparent}

}


%\mode<handout>
%{
%\usepackage{pgfpages}
%\pgfpagesuselayout{4 on 1}[a4paper,landscape,border shrink=5mm]
%}

\usetheme{Antibes}
%\usecolortheme{beetle}
%\usecolortheme{seagull}
\usecolortheme{lily}

\title[Lecture 2]{Statistical methods for linguistic research: Foundational Ideas}

\author{Shravan Vasishth}

\institute{Universit\"at Potsdam\\
vasishth@uni-potsdam.de\\
http://www.ling.uni-potsdam.de/$\sim$vasishth
}

\date{\today}

\begin{document}
\maketitle



<<setup,include=FALSE,cache=FALSE>>=
library(knitr)
library(coda)

# set global chunk options, put figures into folder
options(replace.assign=TRUE,show.signif.stars=FALSE)
opts_chunk$set(fig.path='figures/figure-', fig.align='center', fig.show='hold')
options(replace.assign=TRUE,width=75)
opts_chunk$set(dev='postscript')
#library(rstan)
#set.seed(9991)
# save workspace image, if you want
#the.date <- format(Sys.time(), "%b%d%Y")
#save.image(file=paste0("homework01-",the.date,".RData")
@

\section{The story so far}

\begin{frame}

\begin{enumerate}
\item We defined random variables.
\item We learnt about pdfs and cdfs, and learnt how to compute $P(X<x)$.
\item We learnt about Maximum Likelihood Estimation.
\item We learnt about the sampling distribution of the sample means.
\end{enumerate}

This prepares the way for null hypothesis significance testing (NHST).

\end{frame}

\section{Statistical inference}

\begin{frame}[fragile]\frametitle{Hypothesis testing}

Suppose we have a random sample of size $n$, and the data come from a $N(\mu,\sigma)$ distribution. 

We can estimate sample mean $\bar{x}=\hat \mu$ and $\hat\sigma$, which in turn allows us to estimate the sampling distribution of the mean under (hypothetical) repeated sampling:

\begin{equation}
N(\bar{x},\frac{\hat \sigma}{\sqrt{n}})
\end{equation}

\end{frame}

\begin{frame}[fragile]\frametitle{The one-sample hypothesis test}

Imagine taking an \textbf{independent} random sample from a random variable X that is normally distributed, with mean 12 and standard deviation 10, sample size 11. We estimate the mean and SE:

<<>>=
sample <- rnorm(11,mean=12,sd=10)
(x_bar<-mean(sample))
(SE<-sd(sample)/sqrt(11))
@ 

\end{frame}

\subsection{The one-sample t-test}

\begin{frame}[fragile]\frametitle{The one-sample test}

The NHST approach is to set up a null hypothesis that $\mu$ has some fixed value. For example:

\begin{equation}
H_0: \mu = 0
\end{equation}

This amounts to assuming that the true distribution of sample means is (approximately*) normally distributed and centered around 0, \textit{with the standard error estimated from the data}.

\bigskip
* I will make this more precise in a minute.

\end{frame}


\begin{frame}[fragile]\frametitle{Null hypothesis distribution}

<<nullhyp,echo=FALSE,fig.height=4>>=
x<-seq(-20,20,by=0.1)
plot(x,dt(x,df=10),type="l",main="",
     ylab="density")
points(x_bar/SE,0,col="red",pch=20)
text(x=x_bar/SE,y=0.05,label="sample mean",col="red")
@

\end{frame}

\begin{frame}[fragile]\frametitle{NHST}

The intuitive idea is that 

\begin{enumerate}
\item
if the sample mean $\bar{x}$ is near the hypothesized $\mu$ (here, 0), the data are (possibly) ``consistent with'' the null hypothesis distribution.
\item
if the sample mean $\bar{x}$ is far from the hypothesized $\mu$, the data are inconsistent with the null hypothesis distribution.
\end{enumerate}

We formalize ``near'' and ``far'' by determining how many standard errors the sample mean is from the hypothesized mean:

\begin{equation}
t \times SE = \bar{x} - \mu 
\end{equation}

This quantifies the distance of sample mean from $\mu$ in SE units.

\end{frame}

\begin{frame}[fragile]\frametitle{NHST}

So, given a sample and null hypothesis mean $\mu$, we can compute the quantity: 

\begin{equation}
t  = \frac{\bar{x} - \mu}{SE}
\end{equation}

Call this the \textbf{t-value}. Its relevance will just become clear.

%calculating the probability of obtaining a sample mean $\bar{x}$ (or something more extreme), under this null hypothesis distribution.

\end{frame}

\begin{frame}[fragile]\frametitle{NHST}

The quantity

\begin{equation}
T  = \frac{\bar{X} - \mu}{SE}
\end{equation}

has a t-distribution, which is defined in terms of the sample size $n$. 
We will express this as: $T \sim t(n-1)$ 

Note also that, for large $n$, $T\sim N(0,1)$. 


\end{frame}

\begin{frame}[fragile]\frametitle{NHST}

Thus, given a sample size $n$, and given our null hypothesis, we can draw t-distribution corresponding to the null hypothesis distribution.

For large $n$, we could even use N(0,1), although it is traditional in psychology and linguistics to always use the t-distribution no matter how large $n$ is.

\end{frame}


\begin{frame}[fragile]\frametitle{The t-distribution vs the normal}

\begin{enumerate}
\item The t-distribution takes as parameter the degrees of freedom $n-1$, where $n$ is the sample size (cf.\ the normal, which takes the mean and variance/standard deviation).
\item The t-distribution has fatter tails than the normal for small $n$, say $n<20$, but for large n, the t-distribution and the normal are essentially identical.
\end{enumerate}

\end{frame}

\begin{frame}[fragile]\frametitle{The t-distribution vs the normal}

<<tversusnorm,echo=FALSE,fig.height=4>>=
range <- seq(-4,4,.01)  

op<-par(mfrow=c(2,2),pty="s")

op<-par(mar=c(2,2,3,2),pty="s")

 for(i in c(2,5,15,20)){
   plot(range,dnorm(range),type="l",lty=1,
        xlab="",ylab="",
        cex.axis=1,cex.axis=0.8)
   lines(range,dt(range,df=i),lty=2,lwd=1)
   mtext(paste("df=",i),cex=1.2)
 }
@

\end{frame}

\begin{frame}[fragile]\frametitle{t-test: Rejection region}

So, the null hypothesis testing procedure is:

\begin{enumerate}
\item Define the null hypothesis: for example, $H_0: \mu = 0$.
\item Given data of size $n$, estimate $\bar{x}$, standard deviation $s$, standard error $s/\sqrt{n}$.
\item Compute the t-value:

\begin{equation}
t=\frac{\bar{x}-\mu}{s/\sqrt{n}}
\end{equation}
\item Reject null hypothesis if t-value is large (to be made more precise next).
\end{enumerate}

\end{frame}


\begin{frame}[fragile]\frametitle{t-test}

How to decide when to reject the null hypothesis? Intuitively, when the t-value from the sample is so large that we end up far in \textit{either} tail of the distribution.

\end{frame}

\begin{frame}[fragile]\frametitle{t-test}

<<nullhyprepeat,echo=FALSE,fig.height=4>>=
x<-seq(-20,20,by=0.1)
plot(x,dt(x,df=10),type="l",main="t(n-1)",
     ylab="density")
points(x_bar/SE,0,col="red",pch=20)
text(x=x_bar/SE,y=0.01,label="sample mean",col="red")

lower<-qt(0.025,df=10)
upper<-qt(0.975,df=10)
abline(v=lower)
abline(v=upper)

x1 <- seq(upper,20,abs(0.975)/5)
y1 <- dt(x1,df=10)
polygon(c(x1, rev(x1)), 
        c(rep(0, length(x1)), rev(y1)), 
        col = gray(0.3))

x1 <- seq(-20,lower,abs(0.975)/5)
y1 <- dt(x1,df=10)
polygon(c(x1, rev(x1)), 
        c(rep(0, length(x1)), rev(y1)), 
        col = gray(0.3))

@

\end{frame}

\begin{frame}[fragile]\frametitle{Rejection region}

\begin{enumerate}
\item
For a given sample size n, we can identify the ``rejection region'' by using the \texttt{qt} function (see lecture 1).
\item
Because the shape of the t-distribution depends on the degrees of freedom (n-1), the \textbf{critical t-value} beyond which we reject the null hypothesis will change depending on sample size. 
\item
For large sample sizes, say $n>50$, the rejection point is about 2.
\end{enumerate}

<<>>=
abs(qt(0.025,df=15))
abs(qt(0.025,df=50))
@


\end{frame}


\begin{frame}[fragile]\frametitle{t-test: Rejection region}

Consider the t-value from our sample in our running example:

<<>>=
## null hypothesis mean:
mu<-0
(t_value<-(x_bar-mu)/SE)
@

Recall that the t-value from the sample is simply telling you the distance of the sample mean from the null hypothesis mean $\mu$ in standard error units.

\begin{equation}
t=\frac{\bar{x}-\mu}{s/\sqrt{n}} \hbox{ or } t\frac{s}{\sqrt{n}}=\bar{x}-\mu
\end{equation}

\end{frame}

\begin{frame}[fragile]\frametitle{t-test: Rejection region}

So, for large sample sizes, if $\mid t\mid >2$ (approximately), we can reject the null hypothesis. 

For a smaller sample size $n$, you can compute the exact critical t-value:

\begin{verbatim}
qt(0.025,df=n-1)
\end{verbatim}

This is the critical t-value on the \textbf{left}-hand side of the t-distribution.
The corresponding value on the right-hand side is:

\begin{verbatim}
qt(0.975,df=n-1)
\end{verbatim}

Their absolute values are of course identical (the distribution is symmetric).

\end{frame}

\begin{frame}[fragile]\frametitle{The t-distribution vs the normal}

Given the relevant degrees of freedom,
one can compute the area under the curve as for the Normal distribution:

<<probttest>>=
pt(-2,df=10)
pt(-2,df=20)
pt(-2,df=50)
@

Notice that with large degrees of freedom, the area under the curve to the left of -2 is approximately 0.025.

\end{frame}


\begin{frame}[fragile]\frametitle{The t.test function}

The t.test function in R delivers the t-value:

<<>>=
## from t-test function:
## t-value
t.test(sample)$statistic
@

\end{frame}


\section{Type I, Type II, power}


\begin{frame}[fragile]\frametitle{Type I, Type II error, power}

 When we do a hypothesis test, the sample mean 
\begin{enumerate}
\item
will either fall in the rejection region $\rightarrow$ reject null
\item 
or it will not  $\rightarrow$ fail to reject null
\end{enumerate}

But the null hypothesis is either true or not true.
\textit{We don't know which of those two is the reality}.

\end{frame}
\begin{frame}[fragile]\frametitle{Type I, Type II error, power}

\begin{tabular}{ccc}
\hline
Reality: & $H_0$ TRUE & $H_0$ FALSE \\
\hline
Decision: `reject': & $\alpha$ & $1~-~\beta$ \\
                                     & \textbf{Type I error}                         & \textbf{Power} \\                                      
                                     & & \\
\hline
Decision: `fail to reject': & $1 - \alpha$ & $\beta$ \\                                    &                                 & \textbf{Type II error}\\
\hline
\end{tabular}

\bigskip

Consider the situation where the true $\mu=2$. Now the $H_0$ is false.


\end{frame}

\begin{frame}[fragile]\frametitle{Type I, Type II error, Power}

Type I error is conventionally held at 0.05. Power is 1-Type II error.

<<echo=FALSE,fig.height=4>>=
## function for plotting CDF:
plot.prob<-function(x,
                           x.min,
                           x.max,
                           prob,
                           mean,
                           sd,
                           gray.level,main){

        plot(x,dnorm(x,mean,sd), 
                     type = "l",xlab="",
             ylab="",main=main)
        abline(h = 0)

## shade X<x    
    x1 = seq(x.min, qnorm(prob), abs(prob)/5)
    y1 = dnorm(x1, mean, sd)

    polygon(c(x1, rev(x1)), 
            c(rep(0, length(x1)), rev(y1)), 
            col = gray.level)
  }

shadenormal<- 
function (prob=0.5,
          gray1=gray(0.3),
          x.min=-6,
          x.max=abs(x.min),
          x = seq(x.min, x.max, 0.01),
          mean=0,
          sd=1,main="P(X<0)") 
{

     plot.prob(x=x,x.min=x.min,x.max=x.max,
               prob=prob,
                      mean=mean,sd=sd,
     gray.level=gray1,main=main)     
}
shadenormal(prob=0.025,main="Type I, II error")

x1 <- seq(qnorm(0.975),6,abs(0.975)/5)
y1 <- dnorm(x1)

polygon(c(x1, rev(x1)), 
        c(rep(0, length(x1)), rev(y1)), 
        col = gray(0.3))

x<-seq(-6,6,by=0.1)
lines(x,dnorm(x,mean=2),col="red",lwd=2)
abline(v=2)
abline(v=-2)

x1 <- seq(-2,2,0.01)
y1 <- dnorm(x1,mean=2)

polygon(c(x1, rev(x1)), 
        c(rep(0, length(x1)), rev(y1)), 
        col = "orange")

@

\end{frame}

\subsection{Example}

\begin{frame}[fragile]\frametitle{The typical statistical test}
\framesubtitle{t-test}

Given data:

<<>>=
## Sampling from Normal(0,1)
(sample<-rnorm(10))
@

\end{frame}

\begin{frame}[fragile]\frametitle{The typical statistical test}
\framesubtitle{t-test}


If we do a t-test to test the null hypothesis that $\mu = 0$:

<<>>=
n<-length(sample)
x_bar<-mean(sample)
stddev<-sd(sample)
(t_value<- (x_bar - 0)/(stddev/sqrt(n)))
@

\end{frame}

\begin{frame}[fragile]\frametitle{The typical statistical test}
\framesubtitle{t-test}

We can also compute the probability of getting a t-value like the one we got from the sample, or something more extreme, given the null hypothesis.

This can be computed, as done earlier, simply by calculating the area under the curve in the rejection region for the relevant t-distribution:

<<>>=
pt(-abs(t_value),df=n-1)
@

\end{frame}

\begin{frame}[fragile]\frametitle{The typical statistical test}
\framesubtitle{t-test}


I just took the absolute t-value from the sample and took its negation in order to compute the probability on the left tail. I could have also written:

<<>>=
pt(abs(t_value),df=n-1,lower.tail=FALSE)
@


\end{frame}

\begin{frame}[fragile]\frametitle{The typical statistical test}
\framesubtitle{t-test}

The convention is to compute the probability of getting a t-value like the one we got or something more extreme --- we look at both sides of the t-distribution.

<<>>=
2*pt(-abs(t_value),df=n-1)
@

Conventionally, we reject the null if $p<0.05$. This is because we set the Type I error at 0.05.

\end{frame}


\begin{frame}[fragile]\frametitle{The typical statistical test}

<<>>=
2*pt(-abs(t_value),df=n-1)
@

Conventionally, we reject the null if this probability $<0.05$. 

This probability is called the p-value.

\end{frame}

\begin{frame}[fragile]\frametitle{The typical statistical test}
\framesubtitle{t-test}

You can use the built-in function in R to do such a t-test:

<<>>=
t.test(sample)
@

\end{frame}



\begin{frame}[fragile]\frametitle{Some cautionary notes about the p-value}

The p-value is widely misunderstood, even by veteran scientists. 

Here are some things people \textbf{incorrectly} think is true of p-values:

\medskip

\textbf{Mistake}: A lower p-value gives me more confidence in the specific alternative hypothesis I am interested in verifying.

In fact, a lower p-value only gives me stronger evidence against the null; it doesn't necessarily give me any more evidence than p=0.05 for my \textbf{specific} favored alternative. 

\end{frame}

\begin{frame}[fragile]\frametitle{Some cautionary notes about the p-value}

\textbf{Mistake}: A p-value greater than 0.05 tells me that the null hypothesis is true.

Psychology and linguistics is littered with invalid claims like these. The issue here is lack of statistical power (to be explained next).

\end{frame}

\begin{frame}[fragile]\frametitle{Some cautionary notes about the p-value}

\textbf{Mistake}: It is widely assumed that if $p<0.05$, we have found out that the alternative is true, i.e., that there is a true effect. 

One can always be wrong, typically we allow  probability 0.05 to be wrong.
The only currency we will recognize is replicability.

\end{frame}

\begin{frame}[fragile]\frametitle{Some cautionary notes about the p-value}

\textbf{Mistake}: It is widely believed that the p-value is the probability of the null hypothesis being true.

The p-value is a \textbf{conditional} probability: the probability of seeing the t-value that you got, or something more extreme, assuming the null is true.

An analogy: If you know that the probability of the streets being wet given that it is raining is 0.99, it does not mean that the probability that it's raining is 0.99.

\end{frame}

\begin{frame}[fragile]\frametitle{Type I error vs p-values}

Type I error is the probability of your incorrectly rejecting the null under repeated sampling. We can simulate this:

<<cache=TRUE>>=
nsim<-10000
n<-10
pvals<-rep(NA,nsim)
for(i in 1:nsim){
  x<-rnorm(n)
  pvals[i]<-t.test(x)$p.value
}
mean(pvals<0.05)
@

The single p-value you get from one experiment is just that, a single value. It will vary from experiment to experiment.

\end{frame}

\begin{frame}[fragile]\frametitle{Type II error (1-power) vs p-values}

\begin{enumerate}
\item Most studies in linguistics and psychology have very low power (maybe as low as 0.05).
\item This implies that if we get a so-called null result, i.e., fail to reject the null hypothesis (when $p>0.05$), we can't really conclude anything.
\item If power were high, then a null result could be more meaningful and we might be justified in accepting the null.
\end{enumerate}

But the situation with low power is not just that null results are inconclusive. Even ``statistically significant'' results are suspect with low power.

\end{frame}

\subsection{Type S and Type M errors}

\begin{frame}[fragile]\frametitle{Type S- and M-error}
\framesubtitle{Gelman and Carlin, Beyond Power Calculations:
Assessing Type S (Sign) and Type M (Magnitude) Errors, Perspectives on Psychological Science November 2014 vol.\ 9 no.\ 6 641-651}

If your true effect size is believed to be $D=15$, 
then we can compute (apart from statistical power) these error rates, which are defined as follows:

\begin{enumerate}
\item
\textbf{Type S error}: the probability that the sign of the effect is incorrect, given that (a) the result is statistically significant, or (b) the result is statistically non-significant.
\item 
\textbf{Type M error}: the expectation of the ratio of the absolute magnitude of the effect to the hypothesized true effect size (conditional on whether the result is significant or not). 
Gelman and Carlin also call this the exaggeration ratio, which is perhaps more descriptive than ``Type M error''.
\end{enumerate}


\end{frame}

\begin{frame}[fragile]\frametitle{Type S- and M-error}

Suppose a particular study has standard error 46, and sample size 37. And suppose that our estimated true D=15. Then, we can compute Type S error as follows:

<<typesandm,cache=TRUE,echo=TRUE>>=
## probable effect size derived from past studies:
D<-15
## SE from the study of interest:
se<-46
stddev<-se*sqrt(37)
nsim<-10000
drep<-rep(NA,nsim)
for(i in 1:nsim){
drep[i]<-mean(rnorm(37,mean=D,sd=stddev))
}
@
\end{frame}

\begin{frame}[fragile]\frametitle{Type S- and M-error}

<<typesandm2,cache=TRUE,echo=TRUE>>=
##power:
(pow<-mean(ifelse(abs(drep/se)>2,1,0)))
@

\end{frame}

\begin{frame}[fragile]\frametitle{Type S- and M-error}

<<typesandm3,cache=TRUE,echo=TRUE>>=
## which cells in drep are significant at alpha=0.05?
signif<-which(abs(drep/se)>2)

## Type S error rate | signif:
(types_sig<-mean(drep[signif]<0))
## Type S error rate | non-signif: 
(types_nonsig<-mean(drep[-signif]<0))

@

\end{frame}

\begin{frame}[fragile]\frametitle{Type S- and M-error}

<<typesandm4,cache=TRUE,echo=TRUE>>=
## Type M error rate | signif: 
(typem_sig<-mean(abs(drep[signif])/D))
## Type M error rate | not-signif: 
(typem_nonsig<-mean(abs(drep[-signif])/D))
@

\end{frame}

\begin{frame}[fragile]\frametitle{Type S- and M-error}
\framesubtitle{Funnel plot}

<<funnelplot,echo=FALSE,fig.height=4,cache=TRUE>>=
## funnel plot:
truemu<-15
se<-46
sampsize<-seq(10,10000,by=1)
n_expts<-length(sampsize)
means<-power<-rep(NA,n_expts)
SE<-46
for(i in 1:n_expts){
  x<-rnorm(sampsize[i],mean=truemu,
                       sd=se)
  means[i]<-mean(x)
  power[i]<-power.t.test(d=truemu,sd=250,
                      n=sampsize[i])$power
}

plot(means,power,
     main="Funnel plot",
     xlim=range(c(min(means),max(means))),
     xlab="effect",ylab="power",
     #ylim=c(0,0.003),
     cex.axis=0.8,cex.main=0.8)
abline(v=15)
abline(h=0.20,col="red")
@

\end{frame}

\begin{frame}[fragile]\frametitle{Type S- and M-errors}

\begin{enumerate}
\item
So, you can see that the Type S error and the exaggeration ratio, conditional on a result being significant, are pretty high. 
\item
The practical implication of this is that if most studies are low powered, then it doesn't matter much whether you got a significant result or not. You could be (and probably are) barking up the wrong tree. 
\item
The main point here is: run \textbf{high powered} studies, and \textbf{replicate} the results. There's really nothing that can match consistent replication.
\end{enumerate}

\end{frame}

\subsection{Stopping rules}

\begin{frame}[fragile]\frametitle{Stopping rules}

Psycholinguists and psychologists often adopt the following type of data-gathering procedure:

\begin{enumerate}
\item
The experimenter gathers $n$ data points, then checks for significance ($p<0.05$ or not). 
\item
If it's not significant, he gets more data ($n$ more data points). Since time and money are limited, he might decide to stop anyway at sample size, say, some multiple of $n$. 
\end{enumerate}

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

\begin{enumerate}
\item
One can play with different scenarios here. A typical $n$ might be $15$.
\item
This approach would give us a range of p-values under repeated sampling. 
\item
Theoretically, under the standard assumptions of frequentist methods, we expect a Type I error to be $0.05$. This is the case in standard analyses (I also track the t-statistic, in order to compare it with my stopping rule code below).
\end{enumerate}

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

<<stoppingrules>>=
##Standard:
pvals<-NULL
tstat_standard<-NULL
n<-10
nsim<-10000
## assume a standard dev of 1:
stddev<-1
mn<-0
for(i in 1:nsim){
  samp<-rnorm(n,mean=mn,sd=stddev)
  pvals[i]<-t.test(samp)$p.value
  tstat_standard[i]<-t.test(samp)$statistic
}

@
\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

<<>>=
## Type I error rate: about 5% as theory says:
table(pvals<0.05)[2]/nsim
@

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

But the situation quickly deteriorates as soon as adopt the strategy I outlined above. I will also track the distribution of the t-statistic below.

<<stoppingrules2>>=
pvals<-NULL
tstat<-NULL
## how many subjects can I run?
upper_bound<-n*6
@

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

\small
<<>>=
for(i in 1:nsim){
  significant<-FALSE 
  x<-rnorm(n,mean=mn,sd=stddev) ## take sample
while(!significant & length(x)<upper_bound){
  ## if not significant:
if(t.test(x)$p.value>0.05){
  x<-append(x,rnorm(n,mean=mn,sd=stddev)) ## get more data
} else {significant<-TRUE}   ## otherwise stop:
}
pvals[i]<-t.test(x)$p.value
tstat[i]<-t.test(x)$statistic
}
@

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

<<>>=
## Type I error rate: much higher than 5%:
table(pvals<0.05)[2]/nsim
@

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

Now let's compare the distribution of the t-statistic in the standard case vs with the above stopping rule (red):

<<stoppingrules3,echo=FALSE,fig.height=4>>=
plot(density(tstat_standard),main="",xlab="")
lines(density(tstat),col="red")
#hist(tstat_standard,
#main="The t-distributions for the standard case (white) \n
#     vs the stopping rule (gray)",freq=F)
#hist(tstat,add=T,col="#0000ff22",freq=F)
@

\end{frame}

\begin{frame}[fragile]\frametitle{Stopping rules}

\begin{enumerate}
\item
We get fatter tails with the above stopping rule---a higher Type I error than 0.05.
\item
The point is that one should fix one's sample size in advance based on a power analysis, not deploy a stopping rule like the one above; if we used such a stopping rule, we are much more likely to incorrectly declare a result as statistically significant.
\item
Of course, if your goal is to get an article published no matter what, such stopping rules are a great way to have a successful career!
\end{enumerate}

\end{frame}

\section{Summary so far}

\begin{frame}[fragile]\frametitle{Summary}
\begin{enumerate}
\item We learnt about the single sample t-test.
\item We learnt about Type I, II error (and power).
\item We learnt about Type M and Type S errors.
\end{enumerate}

\end{frame}

\section{Two sample and paired t-tests}

\subsection{Reminder about one-sample t-tests}

\begin{frame}[fragile]\frametitle{t-test}

These are the heights of students in one of my classes at Potsdam:

<<>>=
heights <- c(173,174,160,157,158,170,172,170,
            175,168,165,170,173,180,168,162,
            180,160,155,163,173,175,176,172,
            160,161,150,170,165,184,165)
@

We can do a t-test to evaluate the null hypothesis that $H_0: \mu = 170$ cm.

\end{frame}


\begin{frame}[fragile]\frametitle{The t-distribution}

The formal definition of the t-distribution is as follows:

Suppose we have a random sample of size $n$, say of heights, which come from a $Normal(\mu,\sigma)$ distribution. Then the quantity 

$T=\frac{\overline{X}-\mu}{S/\sqrt{n}}$

has a $\mathsf{t}(\mathtt{df}=n-1)$ sampling distribution. The distribution is defined as ($r$ is degrees of freedom):

$
f_{X}(x,r)=\frac{\Gamma[(r+1)/2]}{\sqrt{r\pi}\ \Gamma(r/2)}\left(1+\frac{x^{2}}{r}\right)^{-(r+1)/2},\quad -\infty < x < \infty.
$

[$\Gamma$ refers to the gamma function; in this course we can ignore what this is, but read Kerns if you are interested.]

\end{frame}


\begin{frame}[fragile]\frametitle{The t-test}

<<>>=
t.test(heights,mu=170)
@

\end{frame}

\begin{frame}[fragile]\frametitle{The t-test}

<<echo=FALSE,fig.width=4>>=
range <- seq(-4,4,.01)  
 
op<-par(mfrow=c(2,2),pty="s")

 for(i in c(2,5,15,20)){
   plot(range,dnorm(range),type="l",lty=1,
        xlab="",ylab="",
        cex.axis=1)
   lines(range,dt(range,df=i),lty=2,lwd=1)
   mtext(paste("df=",i),cex=1.2)
 }
@
\end{frame}



\begin{frame}[fragile]\frametitle{Computing the p-value by hand}

First, we compute $t=\frac{\bar{x}-\mu}{s/\sqrt{n}}$:

<<>>=
(observed_t<-(mean(heights)-170)/(sd(heights)/sqrt(31)))
@

Then we compute the probability of seeing that observed t or something more extreme, assuming the null is true:

<<>>=
2*pt(observed_t,df=30)
@

Notice that for n=31, we could have used the normal distribution:

<<>>=
2*(pnorm(mean(heights),mean=170,sd=sd(heights)/sqrt(30)))
@

\end{frame}

\subsection{The two-sample t-test}

\begin{frame}[fragile]\frametitle{Two-sample t-test}

This is a data-set from Keith Johnson's book (Quantitative Methods in Linguistics):

<<>>=
F1data<-read.table("data/F1_data.txt",header=TRUE)
head(F1data)
@

\end{frame}



\begin{frame}[fragile]\frametitle{Two-sample t-test}

Notice that the male and female values are paired in the sense that they are for the same vowel and language.

We can compare males and females' F1 frequencies, ignoring the fact that the data are paired.

Now, our null hypothesis is $H_0: \mu_m = \mu_f$ or $H_0: \mu_m - \mu_f = \delta = 0$.

\end{frame}

\begin{frame}[fragile]\frametitle{Two-sample t-test}
\framesubtitle{Assuming equal variance between men and women}

<<>>=
t.test(F1data$female,F1data$male,var.equal=TRUE)
@

\end{frame}


\begin{frame}[fragile]\frametitle{Two-sample t-test}

Doing this ``by hand'': The only new thing is the SE calculation, and the the df for t-distribution $(2\times n - 2)=36$.

$SE_\delta 
= \sqrt{\frac{\sigma_1^2}{n_1} + \frac{\sigma_2^2}{n_2}}$

<<>>=
d<-mean(F1data$female)-mean(F1data$male)
(SE<-sqrt(var(F1data$male)/19+var(F1data$female)/19))
observed_t <- (d-0)/SE
2*(1-pt(observed_t,df=36))
@

\end{frame}

\subsection{The paired t-test}

\begin{frame}[fragile]\frametitle{The paired t-test}

But this data analysis was incorrect.

This data is paired: each row has 
F1 measurements from a male and female for the same vowel and language.

For paired data,  $H_0: \delta=0$ as before. But since each row in the data-frame is paired (from the same vowel+language), we subtract row-wise, and get a new vector d with the pairwise differences.

\end{frame}


\begin{frame}[fragile]\frametitle{The paired t-test}

Then, we just do a one-sample test:

<<>>=
diff<-F1data$female-F1data$male
t.test(diff)
@

Notice that the p-value is now ``statistically significant''.
This is because the SE is much smaller:

<<>>=
(SE<-sqrt(var(diff)/19))
(observed_t <- (mean(d)-0)/SE)
2*(1-pt(observed_t,df=18))
@

\end{frame}

\begin{frame}[fragile]\frametitle{Summary so far}

We have worked through the

\begin{enumerate}
\item One sample t-test
\item Two sample t-test
\item Paired t-test
\end{enumerate}

\end{frame}


\subsection{An often-seen mistake in paired t-tests}

\begin{frame}[fragile]\frametitle{A note on paired t-tests}


Note that each row of the data frame cannot have more than one row for a particular pair. 

For example, doing a paired t-test on this frame would be incorrect:

\begin{table}[ht]
\centering
\begin{tabular}{rrrll}
  \hline
female & male & vowel & language \\ 
  \hline
391 & 339 & i & W.Apache \\ 
400 & 320 & i & W.Apache \\ 
\vdots & \vdots & \vdots & \vdots\\
   \hline
\end{tabular}
\end{table}

Why? Because the assumption is that each row is independent of the others. This assumption is violated here.

[In fact, it is arguable whether we can assume that rows containing the same vowel in different languages gives us independence.]

\end{frame}


\begin{frame}[fragile]\frametitle{A note on paired t-tests}

Note that each row of the data frame cannot have more than one row for a particular pair. 

Another example:

\begin{table}[ht]
\centering
\begin{tabular}{rrrll}
  \hline
cond\_a & cond\_b & subject & item \\ 
  \hline
391 & 339 & 1 & 1 \\ 
400 & 320 & 1 & 2 \\ 
\vdots & \vdots & \vdots & \vdots\\
   \hline
\end{tabular}
\end{table}

Here, we have repeated measures from subject 1. The independence assumption is violated.

\end{frame}

\begin{frame}[fragile]\frametitle{A note on paired t-tests}

\begin{enumerate}
\item
What to do when we have repeated measurements from each subject or each item?
\item
We aggregate the data so that each subject (or item) has only one value for each condition.
\item
This has a drawback: it pretends we have one measurement from each subject for each condition.
\item
Later on we will learn how to analyze unaggregated data.
\end{enumerate}


\end{frame}


\begin{frame}[fragile]\frametitle{Example of INCORRECT pair-wise t-test}

We have repeated measures data on noun pronunciation durations, in seconds:

<<>>=
dataN2<-read.table("data/dataN2.txt",header=T)
head(dataN2)
@

\end{frame}



\begin{frame}[fragile]\frametitle{Example of INCORRECT pair-wise t-test}

<<>>=
## significant effect:
with(dataN2,
t.test(N2_dur.2,N2_dur.1,paired=TRUE))
@

\end{frame}

\begin{frame}[fragile]\frametitle{Example of INCORRECT pair-wise t-test}

The above t-test was incorrect because we have multiple rows of (dependent) data from the same subject.

We need to aggregate the multiple measurements from each subject until we have one data point from each subject for each combination of vowel and language.

\end{frame}

\begin{frame}[fragile]\frametitle{CORRECT pair-wise t-test}

First, convert data to ``long'' form:

<<>>=
N2dur1data<-data.frame(item=dataN2$Sentence,
                          subj=dataN2$Speaker_id,
                          cond="a",
                          dur=dataN2$N2_dur.1)
N2dur2data<-data.frame(item=dataN2$Sentence,
                          subj=dataN2$Speaker_id,
                          cond="b",
                          dur=dataN2$N2_dur.2)

N2data<-rbind(N2dur1data,N2dur2data)
@

\end{frame}

\begin{frame}[fragile]\frametitle{CORRECT pair-wise t-test}

<<>>=
head(N2data)
@

\end{frame}

\begin{frame}[fragile]\frametitle{CORRECT pair-wise t-test}

Then aggregate so that we have only one data point per subject for each condition:

<<>>=
N2data_bysubj<-aggregate(dur~subj+cond,mean,
               data=N2data)
@

\end{frame}

\begin{frame}[fragile]\frametitle{Example of CORRECT pair-wise t-test (by items)}

Create a vector for each condition:

<<>>=
conda<-subset(N2data_bysubj,cond=="a")
condb<-subset(N2data_bysubj,cond=="b")
@

\end{frame}

\begin{frame}[fragile]\frametitle{Example of CORRECT pair-wise t-test (by items)}
\framesubtitle{Notice that the result is no longer significant}

<<>>=
## not significant:
t.test(condb$dur,conda$dur,paired=TRUE)
@

\end{frame}

\begin{frame}[fragile]\frametitle{Example of CORRECT pair-wise t-test (by items)}

Alternative syntax:

<<>>=
## alternative syntax:
t.test(dur~cond,paired=TRUE,N2data_bysubj)
@

\end{frame}

\begin{frame}[fragile]\frametitle{Exercise: Do a by items paired t-test}

<<>>=
head(N2data)
@

\end{frame}

\begin{frame}[fragile]\frametitle{Exercise: Do a by items paired t-test}

<<>>=
## STEP 1: Aggregate over items:
#N2data_byitem<-aggregate(dur~item+cond,mean,
#               data=N2data)
@

\end{frame}

\begin{frame}[fragile]\frametitle{Exercise: Do a by items paired t-test}

<<>>=
## STEP 2: Create a vector for condition a and b:
#conda<-subset(N2data_byitem,cond=="a")
#condb<-subset(N2data_byitem,cond=="b")
#conda<-...
#condb<-...
## Do a by subject paired t-test:
#t.test(condb$dur,conda$dur,paired=TRUE)
@
\end{frame}

\end{document}
\begin{frame}[fragile]\frametitle{t-test}


\end{frame}