Null Hypothesis Significance Testing, part 1

One-tailed vs Two-tailed Test

Our alternative hypothesis in this example was that the coin is unfair and \(p(\text{Heads}) > \frac{1}{2}\), which is also called a one-tailed significance test. In such case we are interested only in probabilities of obtaining results more extreme as 5 out of 5 Heads. Due to the experiment design there was no outcome that was more extreme as 5 out of 5 Heads (since we couldn’t get 6 Heads out of 5 coin flips).

Another type of test would be two-tailed significance test. In this way we would be interested in probabilities of obtaining results more or as extreme as obtained data. The alternative hypothesis would be:

• \(H_A\): the coin is unfair, \(p(\text{Heads}) \neq \frac{1}{2}\) (\(p(\text{Heads}) < \frac{1}{2}\) or \(p(\text{Heads}) > \frac{1}{2}\))

As you can see, under the two-tailed testing we are not concerned about the direction of alternative hypothesis.

In our example outcome of 0 Heads (5 Tails) out of 5 flips has the same probability as observed probability of 5 Heads out of 5 flips. Hence,

\[\scriptsize \text{p-value} = P(0) + P(5) = 0.0313 + 0.0313 = 0.0616\]

Since now \(\text{p-value} > \alpha\) we would fail to reject the null hypothesis saying there is not enough evidence to claim that the coin is not fair.

Example

Let’s assume that previous studies suggest that the average time to finish a memory test after 1mg of Alprazolam is 60 seconds. We have observed new data and we want to check if there is enough evidence to claim that the average time to finish a memory test is actually <60 seconds.

• \(H_0\): Average time to finish a memory test after 1mg of Alprazolam is 60 seconds, \(\mu = 60\)
• \(H_A\): Average time to finish a memory test after 1mg of Alprazolam is less than 60 seconds, \(\mu < 60\)
• \(\alpha = 0.05\)

temp_df <- islander_data %>% 
  filter(Drug == "A" & Dosage == 1)

hist(temp_df$Mem_Score_After, 
     main = "Distribution of Memory Test Results\nafter 1mg of Alprazolam",
     xlab = "Average time to finish a memory test")

Summary statistics:

temp_df %>% 
  summarise(mean = mean(Mem_Score_After),
            std = sd(Mem_Score_After),
            n = n()) %>% 
  kable()

Using formula from above we can calculate point estimate (average value of the set), null estimate (60 seconds) and standard error.

alpha <- 0.05
null_estimate <- 60
n <- length(temp_df$Mem_Score_After)
x_bar <- mean(temp_df$Mem_Score_After)
SE <- sd(temp_df$Mem_Score_After) / sqrt(n)

t_stat <- (x_bar - null_estimate) / SE
t_stat

## [1] -0.4255

After we calculated \(t\) statistic we can find the area under the \(t\) distribution curve from \(-\infty\) to \(t_{\text{calc}}\). \(t\) distribution has only one parameter - degrees of freedom, that are equal \(df = n-1\), where \(n\) is the sample size.

df <- n-1
x <- seq(-5,5,0.05)
t_dist <- dt(x, df)

ggplot() +
  geom_line(
    mapping = aes(x = x, y = t_dist),
    color = "black", size = 1.5) +   
  geom_vline(xintercept = t_stat) +
  geom_area(
    mapping = aes(x = x[x <= t_stat], y = t_dist[x <= t_stat]),
    fill="red", alpha=0.6) +
  labs(title = "t-Distribution",
       y = "Density") +
  annotate(
    geom = "curve", x = 2, y = 0.3, xend = t_stat, yend = 0.2, 
    curvature = .3, arrow = arrow(length = unit(2, "mm"))) +
  annotate(geom = "text", x = 2.05, y = 0.3, label = "t statistic", hjust = "left") +
  annotate(
    geom = "curve", x = -2.5, y = 0.2, xend = -1, yend = 0.1, 
    curvature = .3, arrow = arrow(length = unit(2, "mm"))) +
  annotate(geom = "text", x = -2.5, y = 0.22, label = "p-value", hjust = "top") +
  theme_classic()

pval <- pt(t_stat, df)
print(paste0("p-value is: ", round(pval, 3)))
## [1] "p-value is: 0.337"
alpha <- 0.05
nhst_result(pval, alpha)
## [1] "p-value is greater than alpha (0.05). Fail to reject the null hypothesis."

But would you still be able to get the conclusion of the test without software, if you did all the calculations by hand? The answer is “definitely yes”. When you have calculated \(t_\text{stat}\) (which was -0.4255 in our example) you could compare this value to \(t_\text{critical}\). \(t_\text{critical}\) is the \(\alpha\) (for one-tailed test) or \(\alpha/2\) (for two-tailed test) quantile of \(t\)-distribution, that can be easily found in any of t-tables.

We need to find a value of \(t_\text{critical}\) creates an area of 0.05 (our \(\alpha\) level) under the curve for a \(t\) distribution with 22 degrees of freedom:

(t_crit <- qt(alpha, df))

## [1] -1.717

df <- n-1
x <- seq(-5,5,0.05)
t_dist <- dt(x, df)

ggplot() +
  geom_line(
    mapping = aes(x = x, y = t_dist),
    color = "black", size = 1.5) +   
  geom_vline(xintercept = t_stat) +
  geom_vline(xintercept = t_crit, color = "red") +
  labs(title = "t-Distribution",
       y = "Density") +
  annotate(
    geom = "curve", x = 2, y = 0.3, xend = t_stat, yend = 0.2, 
    curvature = .3, arrow = arrow(length = unit(2, "mm"))) +
  annotate(geom = "text", x = 2.05, y = 0.3, label = "t statistic", hjust = "left") +
  annotate(
    geom = "curve", x = -4, y = 0.3, xend = t_crit, yend = 0.2, 
    curvature = .3, arrow = arrow(length = unit(2, "mm"))) +
  annotate(geom = "text", x = -4.05, y = 0.3, label = "t ctitical", hjust = "right") +
  geom_area(
    mapping = aes(x = x[x <= t_stat], y = t_dist[x <= t_stat]),
    fill="red", alpha=0.6) +
  geom_area(
    mapping = aes(x = x[x <= t_crit], y = t_dist[x <= t_crit]),
    fill="blue", alpha=0.6) +
  theme_classic() 

If out calculated \(t\) statistic (black line) would be on the left side from critical \(t\) value (red line) we would reject the null hypothesis, since we would be sure that the p-value (red area) is < 0.05. Here we see that our \(t\) statistic is greater than \(t_\text{critical}\), so we fail to reject the null hypothesis. It’s always a good idea to draw a test distribution.

Alternatively, we can let R and Python do all the calculations for us:

t.test(temp_df$Mem_Score_After, 
       mu = 60, 
       conf.level = 0.95,
       alternative = "less")

## 
##  One Sample t-test
## 
## data:  temp_df$Mem_Score_After
## t = -0.43, df = 22, p-value = 0.3
## alternative hypothesis: true mean is less than 60
## 95 percent confidence interval:
##   -Inf 64.38
## sample estimates:
## mean of x 
##     58.56

islander_data = pd.read_csv("memory-test-on-drugged-islanders-data/Islander_data.csv")
temp_df = islander_data[(islander_data["Drug"] == "A") & (islander_data["Dosage"] == 1)] 

t_stat, p_val = stats.ttest_1samp(a=temp_df['Mem_Score_After'], popmean=60)

print(f"Calculated test statistic: {t_stat: .4f}\np-value: {p_val/2: .4f}")

## Calculated test statistic: -0.4255
## p-value:  0.3373

Example

Let’s check if there is a difference in memory test results after 5 mg of Alprazolam and 3 tabs of sugar pills (placebo).

• \(H_0\): There is no difference in test results fare 5 mg of Alprazolam and 3 tabs of sugar pills, \(\mu_{\text{Alprozam}} - \mu_{\text{Placebo}} = 0\)
• \(H_A\): There is a difference in test results fare 5 mg of Alprazolam and 3 tabs of sugar pills, \(\mu_{\text{Alprozam}} - \mu_{\text{Placebo}} \neq 0\)
• \(\alpha = 0.05\)

islander_data %>% 
  filter(Drug %in% c("A", "S") & Dosage == 3) %>% 
  group_by(Drug) %>% 
  ggplot() +
    geom_boxplot(aes(x = Drug, y = Mem_Score_After, color = Drug), show.legend = FALSE) +
  labs(title = "Distribution of Memory Test Result after Drugs",
       y = "Memory Test Score") +
  theme_classic()

Summary statistics:

islander_data %>% 
  filter(Drug %in% c("A", "S") & Dosage == 3) %>% 
  group_by(Drug) %>% 
  summarise(mean = mean(Mem_Score_After),
            std = sd(Mem_Score_After),
            n = n()) %>% 
  kable()

Test statistic can be calculated in the same way:

\[\scriptsize t_{\text{calc}} = \frac{\text{Point Estimate} - \text{Null Estimate}}{SE}\]

• Point estimate is the difference of two sample means: \(\bar{x}_{\text{Alprozam}} - \bar{x}_{\text{Placebo}}\);
• Null estimate is the value under the null hypothesis: \(\mu_{\text{Alprozam}} - \mu_{\text{Placebo}} = 0\)
• “Pooled” standard error: \(SE = \sqrt{\frac{s^2_\text{Alprazolam}}{n_\text{Alprazolam}} + \frac{s^2_\text{Placebo}}{n_\text{Placebo}}}\)
• Degrees of freedom: \(df = n_{\text{Alprozam}} + n_{\text{Placebo}} - 2\)

alp_sample <- islander_data %>% 
  filter(Drug == "A" & Dosage == 3) %>% 
  select(Mem_Score_After) 
pla_sample <- islander_data %>% 
  filter(Drug == "S" & Dosage == 3) %>% 
  select(Mem_Score_After)

# sample means
x_alp <- mean(alp_sample$Mem_Score_After)
x_pla <- mean(pla_sample$Mem_Score_After)

# sample variance
var_alp <- var(alp_sample$Mem_Score_After)
var_pla <- var(pla_sample$Mem_Score_After)

# sample standard deviations
s_alp <- sd(alp_sample$Mem_Score_After)
s_pla <- sd(pla_sample$Mem_Score_After)

# sample size
n_alp <- length(alp_sample$Mem_Score_After)
n_pla <- length(pla_sample$Mem_Score_After)

# degrees of freedom
df <- n_alp + n_pla - 2 

Check that variances are equal:

For this purpose, we are going to use the proportion of two variances \(\frac{s_1^2}{s_2^2}\) to estimate the proportion of population variances \(\frac{\sigma_1^2}{\sigma_2^2}\) using \(F\)-distribution⁵.

• \(H_0\): population variances are equal \(\frac{\sigma_1^2}{\sigma_2^2} = 1\).
• \(H_A\): population variances are different \(\frac{\sigma_1^2}{\sigma_2^2} \neq 1\).

\(F\)-distribution has two parameters \(df_1 = n_1 - 1\) and \(df_2 = n_2 -1\), where \(n_1\) and \(n_2\) is the size of samples.

# check for equal variances
var_ratio <- var_alp / var_pla
pval <- pf(var_ratio, n_alp-1, n_pla-1, lower.tail = FALSE) * 2
print(paste0("p-value is: ", round(pval, 3)))

## [1] "p-value is: 0.896"

p-value is 0.9 so we failed to reject the null hypothesis that population variances are equal.

var.test(alp_sample$Mem_Score_After, pla_sample$Mem_Score_After)

## 
##  F test to compare two variances
## 
## data:  alp_sample$Mem_Score_After and pla_sample$Mem_Score_After
## F = 1.1, num df = 21, denom df = 21, p-value = 0.9
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4399 2.5518
## sample estimates:
## ratio of variances 
##              1.059

Calculate \(t\)-statistic:

se <- sqrt(s_alp^2/n_alp + s_pla^2/n_pla)
t_stat <- (x_alp - x_pla - 0) / se
t_stat

## [1] 3.542

x <- seq(-5,5,0.01)
t_dist <- dt(x, df)

ggplot() +
  geom_line(
    mapping = aes(x = x, y = t_dist),
    color = "black", size = 1.5) +   
  geom_vline(xintercept = t_stat) +
  geom_vline(xintercept = -t_stat) +
  geom_area(
    mapping = aes(x = x[x <= -t_stat], y = t_dist[x <= -t_stat]),
    fill="red", alpha=0.6) +
  geom_area(
    mapping = aes(x = x[x >= t_stat], y = t_dist[x >= t_stat]),
    fill="red", alpha=0.6) +
  labs(title = "t-Distribution",
       y = "Density") +
  theme_classic()

pval <- pt(-t_stat, df)*2
print(paste0("p-value is: ", round(pval, 3)))
## [1] "p-value is: 0.001"
alpha <- 0.05
nhst_result(pval, alpha)
## [1] "p-value is less than alpha (0.05). Reject the null hypothesis."

t.test(alp_sample$Mem_Score_After,
       pla_sample$Mem_Score_After,
       null_estimate = 0,
       conf.level = 0.95,
       alternative = "two.sided",
       var.equal = TRUE) 

## 
##  Two Sample t-test
## 
## data:  alp_sample$Mem_Score_After and pla_sample$Mem_Score_After
## t = 3.5, df = 42, p-value = 0.001
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   8.58 31.30
## sample estimates:
## mean of x mean of y 
##     77.54     57.60

alp_mask = (islander_data["Drug"] == "A") & (islander_data["Dosage"] == 3)
pla_mask = (islander_data["Drug"] == "S") & (islander_data["Dosage"] == 3)

alp_sample = islander_data["Mem_Score_After"][alp_mask]
pla_sample = islander_data["Mem_Score_After"][pla_mask]

t_stat, p_val = stats.ttest_ind(alp_sample, pla_sample, equal_var=True)

print(f"Calculated test statistic: {t_stat: .4f}\np-value: {p_val: .5f}")

## Calculated test statistic:  3.5422
## p-value:  0.00099

Example

Let’s now check if there is a significant difference (in any direction, so two-tailed test) in memory test results before and after 0.75mg of Triazolam. Now we are going to use the variable Diff which was calculated as Mem_Score_After - Mem_Score_Before.

• \(H_0\): There is no difference in test results before and after 0.75mg of Triazolam, \(\mu_{\text{diff}} = 0\)
• \(H_A\): There is a difference in test results before and after 0.75mg of Triazolam, \(\mu_{\text{diff}} \neq 0\)
• \(\alpha = 0.05\)

temp_df <- islander_data %>% 
  filter(Drug == "T" & Dosage == 3)

hist(temp_df$Diff, 
     main = "Distribution of Difference of Memory Test Results\nbefore and after 0.75mg of Triazolam",
     xlab = "Average time to finish a memory test")

There are two outliers (difference in test results is less than \(-20\)) that need to be removed for further analysis:

temp_df[temp_df$Diff < -20, ] %>% kable()

Summary statistics:

temp_df %>% 
  filter(Diff > -20) %>% 
  summarise(n = n(),
            mean = mean(Diff),
            std = sd(Diff)) %>% 
  kable()

results <- t.test(temp_df$Diff[temp_df$Diff > -20], 
       mu = 0, 
       conf.level = 0.95,
       alternative = "two.sided")

results

## 
##  One Sample t-test
## 
## data:  temp_df$Diff[temp_df$Diff > -20]
## t = 0.24, df = 18, p-value = 0.8
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.644  3.318
## sample estimates:
## mean of x 
##    0.3368

temp_df = islander_data[(islander_data["Drug"] == "T") & (islander_data["Dosage"] == 3)] 
temp_df = temp_df[temp_df["Diff"] > -20]

t_stat, p_val = stats.ttest_1samp(a=temp_df['Diff'], popmean=0)

print(f"Calculated test statistic: {t_stat: .4f}\np-value: {p_val: .4f}")

## Calculated test statistic:  0.2374
## p-value:  0.8150

Since it’s two-tailed test we are interested in area under the curve from two sides of the distribution:

n <- dim(temp_df)[1]
t_stat <- results$statistic
df <- n-1
x <- seq(-5,5,0.01)
t_dist <- dt(x, df)

ggplot() +
  geom_line(
    mapping = aes(x = x, y = t_dist),
    color = "black", size = 1.5) +   
  geom_vline(xintercept = t_stat) +
  geom_vline(xintercept = -t_stat) +
  geom_area(
    mapping = aes(x = x[x <= -t_stat], y = t_dist[x <= -t_stat]),
    fill="red", alpha=0.6) +
  geom_area(
    mapping = aes(x = x[x >= t_stat], y = t_dist[x >= t_stat]),
    fill="red", alpha=0.6) +
  labs(title = "t-Distribution",
       y = "Density") +
  theme_classic()

Example

Let’s look for simplicity at the small number of samples and check whether there is a significant difference in memory test results before and after 0.25mg (Dosage = 1), 0.5mg (Dosage = 2) and 0.75mg (Dosage = 3) of Triazolam (3 samples in total).

islander_data %>% 
  filter(Drug == "T") %>%
  ggplot() +
    geom_boxplot(aes(x = as.factor(Dosage), y = Diff, color = as.factor(Dosage)), 
                 show.legend = FALSE) +
  labs(title = "Distribution of Difference in Memory Test Result\nbefore and after Triazolam",
       y = "Memory Test Score",
       x = "Dosage") +
  theme_classic()

There are some outliers that need to be removed for further analysis:

islander_data %>% 
  filter(Drug == "T" & (Diff < -15 | Diff > 15)) %>% 
  kable()

Summary statistics:

islander_data %>% 
  filter(Drug == "T" & Diff > -15 & Diff < 15) %>%
  group_by(Dosage) %>% 
  summarise(mean = mean(Diff),
            std = sd(Diff),
            n = n()) %>% 
  kable()

## `summarise()` ungrouping output (override with `.groups` argument)

temp_df <- islander_data %>% 
  filter(Drug == "T" & Diff > -15 & Diff < 15)

oned_sample <- temp_df$Diff[temp_df$Dosage == 1]
twod_sample <- temp_df$Diff[temp_df$Dosage == 2]
threed_sample <- temp_df$Diff[temp_df$Dosage == 3]

# check for equal variances
var.test(oned_sample, twod_sample)
## 
##  F test to compare two variances
## 
## data:  oned_sample and twod_sample
## F = 1.5, num df = 21, denom df = 20, p-value = 0.3
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.6309 3.7454
## sample estimates:
## ratio of variances 
##              1.545
var.test(oned_sample, threed_sample)
## 
##  F test to compare two variances
## 
## data:  oned_sample and threed_sample
## F = 0.8, num df = 21, denom df = 18, p-value = 0.6
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.3143 1.9678
## sample estimates:
## ratio of variances 
##             0.7993
var.test(twod_sample, threed_sample)
## 
##  F test to compare two variances
## 
## data:  twod_sample and threed_sample
## F = 0.52, num df = 20, denom df = 18, p-value = 0.2
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2022 1.2944
## sample estimates:
## ratio of variances 
##             0.5175

Total variability (\(SST\), sum of squares total) in the memory test results can be divided into variability attributed to the drug (between group variability) and variability attributed to other factors (within group variability).

\[SST = SSG + SSE\]

\[SST = \sum_{i=1}^n \left[ (y_i - \bar{y})^2 \right]\]

• \(\bar{y}\): total average value.

\(SSG\): the sum of squares for groups, measures the variability between groups. Explained variability.

\[SSG = \sum_{j=1}^k \left[ n_j (\bar{y_j} - \bar{y})^2 \right]\]

• \(n_j\): number of samples in \(j\) group.
• \(\bar{y_j}\): average value for the \(j\) group.

\(SSE\): the sum of squares for error, measures the variability within groups. Unexplained variability.

We can rewrite the formula for \(F\) statistic:

\[F = \frac{MSG}{MSE}\]

\(MSG\): Mean Squares for groups. Average variability between groups calculated as the total variability (sum of squares) scaled by degrees of freedom.

\[MSG = \frac{SSG}{df_{\text{group}}}\]

\(MSE\) Mean Squares for error. Average variability within groups calculated as the total variability (sum of squares) scaled by degrees of freedom.

\[MSE = \frac{SSE}{df_{\text{Error}}}\]

Degrees of freedom can be found as:

• \(df_{\text{total}} = n - 1\)
• \(df_{\text{group}} = k -1\)
• \(df_{\text{error}} = df_{\text{total}} - df_{\text{group}}\)

# calculations

# degrees of freedom
n <- dim(temp_df)[1]
k <- length(unique(temp_df$Dosage))
df_total <- n - 1
df_group <- k - 1
df_error <- df_total - df_group

# sum of squares
y_bar <- mean(temp_df$Diff)
(sst <- sum((temp_df$Diff - y_bar)^2))
## [1] 1757

(temp_df %>% 
  group_by(Dosage) %>% 
  summarise(n = n(),
            yj_bar = mean(Diff)) %>% 
  mutate(ssg = n * (yj_bar - y_bar)^2) %>% 
  summarise(sum(ssg)) %>% 
  as.numeric() -> ssg)
## `summarise()` ungrouping output (override with `.groups` argument)
## [1] 29.84

(sse <- sst - ssg)
## [1] 1727

# average variability
(msg <- ssg / df_group)
## [1] 14.92
(mse <- sse / df_error)
## [1] 29.27

F_stat <- msg / mse
print(paste0("Calculated F statistic: ", round(F_stat, 3)))
## [1] "Calculated F statistic: 0.51"

x <- seq(0,4,0.01)
f_dist <- df(x, df_group, df_error)

ggplot() +
  geom_line(
    mapping = aes(x = x, y = f_dist),
    color = "black", size = 1.5) +   
  geom_vline(xintercept = F_stat) +
  geom_area(
    mapping = aes(x = x[x >= F_stat], y = f_dist[x >= F_stat]),
    fill="red", alpha=0.6) +
  labs(title = "F-Distribution",
       y = "Density") +
  theme_classic()

pval <- pf(F_stat, df_group, df_error, lower.tail = FALSE)
print(paste0("p-value is: ", round(pval, 3)))
## [1] "p-value is: 0.603"
alpha <- 0.05
nhst_result(pval, alpha)
## [1] "p-value is greater than alpha (0.05). Fail to reject the null hypothesis."

As you can see there is a lot of calculations, so most of the time we rely on software:

aov(formula = Diff ~ as.factor(Dosage),
    data = temp_df) %>% 
  tidy() %>% 
  kable()

temp_df = islander_data[(islander_data["Drug"] == "T") & (islander_data["Diff"].between(-15, 15))]
F_stat, p_val = stats.f_oneway(
  temp_df["Diff"][temp_df["Dosage"] == 1], 
  temp_df["Diff"][temp_df["Dosage"] == 2],
  temp_df["Diff"][temp_df["Dosage"] == 3])
  
print(f"Calculated test statistic: {F_stat: .4f}\np-value: {p_val: .4f}")

## Calculated test statistic:  0.5098
## p-value:  0.6033