```python
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import pearsonr
from scipy.stats import t as tdist
```

# DIY t-test

As a first exercise lets try to mimic a t-test using simulated data. As the t-test is just a theoretic representation of what we are doing here, our empirical result from the simulation should be the very close to the theoretical one. To generate a so called null-distribution for the correlation coefficient we repeatedly draw __uncorrelated__ samples from a normal distribution and correlate them with each other saving the result.


```python
# Total number of samples
nsample = 10000
# Number of observations in each sample
nobs = 100 

# Initialization of a variable for all our samples
r_sample= np.zeros(nsample)

# Draw samples and correlate
for i in range(nsample):
    x = np.random.randn(nobs)
    y = np.random.randn(nobs)
    r, p = pearsonr(x, y)
    r_sample[i] = r
```

We can take a look at the samples and how they are distributed using a histogram. Note that they have a mean value of 0. The width of the distribution changes with the number of observations we correlate in each of the samples. Go ahead and give that a try, by changing the value of `nobs` above.


```python
plt.hist(r_sample, histtype='step')
plt.xlabel('$r$')
```

For the t-test we use the value of $r$ and the number of observations to calculate the test statistic, the t-value using the following formula:

\begin{align}
    t = r \frac{\sqrt{n - 2}}{\sqrt{1 - r^2}}
\end{align}

below this formula is implemented in a function and than applied to the sample of $r$'s that we have produced above.


```python
def calc_t(r, n):
    """Calculate t statistic for Pearsons r"""
    t = r * np.sqrt(n - 2) / np.sqrt(1 - r**2)
    return t
```


```python
t_sample = calc_t(r_sample, nobs)
```

To be able to compare the sample with the theoretical distribution of this variable, the t-distribution we need to have a range of t-values that spans the range of our samples:


```python
t_plot = np.linspace(np.min(t_sample), np.max(t_sample), 100)
```

Now we have everything to compare our empirical distribution with the theoretical one:


```python
plt.hist(t_sample, histtype='step', density=True, label='empirical')
plt.plot(t_plot, tdist.pdf(t_plot, nobs - 2), label='theoretical')
plt.legend()
```

They should be a pretty close match as the theory is exactly describing what we have done before, just for infinitly many samples.

Lets now compare the $r$ value of the time series in the slides to both the theoretical and empirical ones. Recall. that $r=0.27$ and we had 99 observations in the correlation.


```python
r_data = 0.27
nobs = 99
```

First lets test the value using the classical t-test:


```python
t_data = calc_t(r_data, nobs)
print('t = %.2f' % t_data)

p_theoretical = 2 * (1 - tdist.cdf(t_data, nobs - 2))

print('Theoretical p-value: %.3f' % p_theoretical)
```

Now lets use the code from above to generate a sample for r and for t that we can test these values against.


```python
# Total number of samples
nsample = 10000
# Number of observations in each sample
# Initialization of a variable for all our samples
r_sample= np.zeros(nsample)
# Draw samples and correlate
for i in range(nsample):
    x = np.random.randn(nobs)
    y = np.random.randn(nobs)
    r, p = pearsonr(x, y)
    r_sample[i] = r
    
t_sample = calc_t(r_sample, nobs)
```

We can test both the fraction of simulated t-values that are larger than the t-value for the data as well as the r-value directly to obtain an empirical p-value for our data.


```python
np.mean(np.abs(t_sample) >= t_data)
```


```python
np.mean(np.abs(r_sample) >= r_data)
```

Again, these values should be very close to the theoretical one as they are essenitally generated the same way, just with a finite number of samples.

So why was this important, when the result is the same as in the t-test anyway?

The t-test as per the theory is not always applicable, depending on the data or the type of statistic that we are looking at. The same is true for any other statistical test. However, as long as you can generate an empirical null-distribution, you will have something to test against. You will not be limited to the special cases that the usual tests are usefull for! So this a very important tool to have.

0dcca953c1096e1018c8816bd486d9d29066e29d
ipynb
Jupyter Notebook
01-DIY_t-test.ipynb
terhardt/ICAT-stats
5e74aff69a4c1b65c7756cb4edcd0db1df8a60fc
MIT
01-DIY_t-test.ipynb
terhardt/ICAT-stats
5e74aff69a4c1b65c7756cb4edcd0db1df8a60fc
MIT
01-DIY_t-test.ipynb
terhardt/ICAT-stats
5e74aff69a4c1b65c7756cb4edcd0db1df8a60fc
MIT
Qwen/Qwen-72B
1. YES
2. YES