Comparing two populations through samples

Setting up example

Let's define two populations as following Normal (Gaussian) distributions and imagine taking samples from them. These represent the lengths of two whitefish populations in the Baltic Sea.

Population 1: $y_{i} \sim N (μ_{1}, σ^{2})$
Population 2: $z_{i} \sim N (μ_{2}, σ^{2})$
$μ_{1} = 20, μ_{2} = 30, σ = 10$
These two populations have different means but the same variance.

A plot of their probability distribution functions, true means, and values of $n = 15$ random sampls from each and their sample means is given in the lecture notes:
!lecture_notes1_week1, p.4
We end up with sample means of $\bar{y} = 25.1$ and $\bar{z} = 31.4$

Question

If we didn't know the true means, would the data given by our $n = 15$ random samples be enough to prove that the true means are different?

The sample means are called Estimators and they are a mathematical tool to estimate the true population parameters. Let's take a look at these values:

Distribution of sample means

If we repeat this random sampling 100 times, then take the average of all of these, the value will be closer to the true population mean. The lecture notes get into some other techniques in R for estimating the true mean from samples.

What's important from this is that estimators have their own variances.

Difference between populations

By analyzing the distribution (variance) of our samples, we can determine how confident we are that there is a difference between the true means of the populations (not just the sample means).

Back to our example. Let's consider $y_{i}$ and $y_{i^{'}}$ , two random independent samples from population 1. Since we know that the trait values in $y$ are Gaussian (normally) distributed, we can solve analytically that:

y_{i} - y_{i^{'}} \sim N (0, 2 σ^{2})

Why does it have this distribution?

Since $y_{i}$ and $y_{i^{'}}$ both come from the same population, they have the same mean, so subtracting one from the other should give a mean of zero.
- This is why the mean is 0 in the distribution
Variance is additive and has no positive or negative sign. So when you subtract one sample from the other, the variances will actually add on to each other as long as they are independent (which they are in this case).
- This is why the variance is $2 σ^{2}$ .
- See Probability - Variance#Sum of variances

Now let's consider the difference between a random sample from population 1 and population 2:

y_{i} - z_{i} \sim N (μ_{1} - μ_{2}, 2 σ^{2})

We are also assuming here that both populations are independent and have the same variance, which is why the variance is doubled when subtracting them.

Difference between sample means

The distribution of differences between the sample means ( $\bar{y}$ and $\bar{z}$ ) in our samples from populations 1 and 2 is:

\bar{y} - \bar{z} \sim N (μ_{1} - μ_{2}, \frac{2 σ^{2}}{n})

Why is the variance divided by

n

The Variance of the sample average is $V (\bar{ξ}) = \frac{1}{n} V (ξ)$ . If both our samples use the same $n$ and are independent, then we have:

V (\bar{y} - \bar{z}) = V (\bar{y}) + V (\bar{z})

V (\bar{y} - \bar{z}) = \frac{1}{n} V (y_{i}) + \frac{1}{n} V (z_{i})

V (\bar{y} - \bar{z}) = \frac{1}{n} (σ^{2} + σ^{2})

V (\bar{y} - \bar{z}) = \frac{2 σ^{2}}{n}

Since $σ = 10$ and $n = 15$ , the variance is 13.33 and standard deviation is 3.7

Cumulative Distribution Function

This kind of function tells us what the probability that a random variable is less than or equal to a certain value. See Continuous Distributions for more info.

Let's define $d = \bar{y} - \bar{z}$ and call the cumulative function $F_{d} (x)$ as the cumulative probability that $d$ is less than the value $x$ . This is defined as:

F_{d} (0) = Probability (d \leq 0) = \int_{- \infty}^{0} p (d) d d

Where $p (d)$ is the probability density function (distribution) of $d$ - in other words:

p (d) = N (μ_{1} - μ_{2}, \frac{2 σ^{2}}{n})

So, $F_{d} (0)$ is the integral (area under the curve) taken from minus infinity to zero of this function. Drawing:

This cumulative distribution function is useful for:

Frequentist inference

It lets us quantify how likely it is to observe a specific value in an estimator (in our case the difference between sample means)

Bayesian inference

It lets us calculate the probability that some parameter has a certain value. Not sure what this means yet, we'll look at this more later.

Frequentist hypothesis testing

Now we can figure out how likely it is that the two population trait means are different based only on our two samples.

Hypotheses:

H0: The two samples come from the same population distribution
- e.g. they have the same mean and standard deviation
- We can call this common distribution $N (μ, σ^{2})$ and assume the variance is known but the man is unknown.
- If this null hypothesis is true, then the difference between our two sample means would be normally distributed around zero, specifically according to the following equation:

\bar{y} - \bar{z} \sim N (0, \frac{2 σ^{2}}{n})

H1: The two samples come from different population distributions
- e.g. they might have different means or standard deviations

Assumption of distributions

Note that both hypotheses assume we have a normal (Gaussian) distribution of both populations and samples. We will return to this note later when we look at the Bayesian Approach to Statistics.

To evaluate how likely that the difference in means arises under the null hypothesis lets introduce two random vectors $Y$ and $Z$ that collect $n$ independent samples from the null distribution.

Then, we define another random variable $T$ as $T = \bar{Y} - \bar{Z}$ - equal to the difference between sample means in the null distribution.

Note

The difference between $y, z$ and $Y, Z$ is that $Y$ and $Z$ are known to come from the null distribution, but $y$ and $z$ may or may not have come from this, that's what we're trying to figure out.

(in reality, they came from different distributions, but we are pretending we don't know that for this exercise).

Now we can get the probability that the difference in population means ( $t = \bar{y} - \bar{z}$ ) is at least as extreme as the observed difference between our sample means under the null hypothesis:

Probability (| T | \geq | t |) = Pr (| \bar{Y} - \bar{Z} | \geq | \bar{y} - \bar{z} |)

This will give us a p-value, or probability that the test statistic $T$ is as least as extreme as test statistic $t$ . This is a two-tailed t-test. The lecture notes show how to evaluate this in R:
!lecture_notes1_week1, p.12

We can compare our p-value to a threshold value called $α$ or a significance level. Traditional statistical analysis holds that we reject our null hypothesis if the p value is smaller than the significance level, and accepted otherwise.

Remember:

p-value does not tell us the probability there is a difference between the two populations
p-value does not tell us how big the difference between populations is
The significance level $α$ is an arbitrary choice and should be considered carefully before just relying on $α = 0.05$ .
It's not recommended to test hypotheses and reject null hypotheses solely based on significance levels, but p-values are still useful for assessing statistical support for a phenomenon.

Student's t-test

Since variances are usually unknown, we need to estimate it using the sample variance. If the two populations are assumed to have different population variances, we need to use the Student's t-test (sample variance of population 1 = $s_{y}^{2}$ ):

\frac{\bar{y} - \bar{z}}{\sqrt{s_{y}^{2} / n + s_{z}^{2} / n}} \sim {Student-t}_{v = n - 2} (0, 1)

Confidence intervals

The output in R of the t-test mentions 95 percent confidence intervals. These are a measurement of uncertainty in an estimate.

Confidence level $γ$ is a value between zero and one that tells how confident we want to be about the statement of the true value of the parameter of interest (for example the difference in population means). It is often defined to be a number close to one – most typically 95%.

Confidence intervals will be denoted by $[l (D), u (D)]$ :

$D$ means the data (in our case, $D = {y, z}$ , e.g. all the population 1 and 2 data)
$l (D)$ and $u (D)$ are functions that determine the lower and upper ends of the confidence interval. Since they are functions of our data, they are also random variables (like the test statistics).

Confidence interval for a parameter $θ$ , with confidence level $γ$ is an interval $[l (D), u (D)]$ , such that

Pr (l (D) \leq θ \leq u (D)) = γ

under the statistical model used to describe data.

(The probability of $θ$ being between the two boundaries is equal to our Confidence level $γ$ , usually 95%.)

Implementing confidence intervals

Based on our probability density function $p (d)$ from earlier, we can calculate many iterations of getting the upper and lower bounds from $Y$ and $Z$ , random samples generated from our null model. We are skipping the details of how to actually solve the $l (D)$ and $u (D)$ functions for now, but this is how to do it in R:
!lecture_notes1_week1, p.15
Or:

l (D) = (\bar{y} - \bar{z}) - (quantile at 0.975) \frac{(s_{y} - s_{z})}{\sqrt{n}}

l (U) = (\bar{y} - \bar{z}) + (quantile at 0.975) \frac{(s_{y} - s_{z})}{\sqrt{n}}

Not exactly sure how the quantile function works.

Conclusion

This is one method of frequentist hypothesis testing to see if two populations have different means.

The qualitative part of the test comes from the t-test, since it doesn't tell us how big the difference is, just how likely it is to be different.

The Confidence interval gives us both a qualitative comparison of whether a difference exists, and a quantitative estimate of the size of the difference.

Next, we will look at using linear models to compare population means.

Linear models for comparing population means

Definitions

Linear models are another way to compare population means through regression.

First, let's define $y$ as a vector of all of the observations for populations 1 and 2:

y = [y_{1}, y_{2} . . . y_{n}, z_{1}, z_{2}, . . ., z_{n}]

Then we define a covariate (differenct from covariance) vector $x$ which is just $0$ for all the $y_{i}$ values and $1$ for all the $z_{i}$ values:

x = [0, 0, . . .0, 1, 1, . . ., 1]

(The first $n$ values of the vector $x$ are 0, the last $n$ elements are $1$ ).

The linear model for comparing our groups is:

Y_{i} = β_{0} + β_{1} x_{i} + ϵ_{i}

The intercept $β_{0}$ is $μ_{1}$ and the coefficient on $x$ $β_{1} = μ_{2} - μ_{1}$ (also known as $d$ ).

Implementation in R

# Create vector of observed outcome variables (=fish lenghts) 
Y = c(y,z) 
# create a vector of covariates (=0 or 1 to indcate sampling location) 
X = c(0*y,z/z)


# And then visualize the data in "regression manner" 
plot(X,Y) 
# overlay the data with sample means 
points(0,mean(Y[X==0]), pch=17, col="black", cex=2) points(1,mean(Y[X==1]), pch=17, col="blue", cex=2)

!lecture_notes1_week1, p.17
The plot created from this shows the distribution of the two vectors on the y-axis, and assigns them to zero and one based on their category on the x axis. You can start to see how a regression would work on this with the R function lm(). We just need to add 1 to the X-values since I guess this will make the math work better (x-values are 1 and 2 instead of 0 and 1).

model = lm(Y ~ 1+X)
summary(model)

!lecture_notes1_week1, p.18
Based on the summary of our model generated in R, we can see that they gave an estimate of the difference between the two groups as 6.307 (under the coefficients X1 section). Standard error is 3.429. You can also use confint(model) to get the confidence intervals for the model.

Now let's go back to the main lecture note and look at data transformations.
Statistics - Introduction#Data Transformations