Overview of Today’s Lecture

• Non-Parametric Testing Procedures
• Dealing with Multiple Comparisons
• Review of Statistical Tests

Non-Parametric Testing Procedures

• Sign Test
• Komolgorov-Smirnov Test
• Mann-Whitney Rank Sum Test
  • a.k.a. Wilcoxen Rank Sum Test
• Resampling

An Unknown Population Distribution.

We’ve seen in the previous lectures that under certain conditions, the t-test can be used to run hypothesis tests even when the data is not normally distributed. The central limit theorem tells us that when a population is not normally distributed, when you take a large enough sample, the mean of the sample tends to be normal. The larger the sample, the more normal the sample mean will appear.

• When can you use a t-test?
  • when you know the data comes from a normal distribution
    • or
  • when the sample size is large (\(\gtrsim 30\))
• When should you not use a t-test?
  • When you have a small sample from a non normal distribution
• What can be done in these cases?
  • if the population distribution is known, a statistical test based on a sample point estimate can be derived mathematically.
  • if the population distribution is unknown, there are non parametric options.

One Sample Sign Test (Continuous data)

This is a non parametric test that lets us compare the sample median to a known population median. It is a simple alternative to the one sample t-test when you cannot assume normality of the data and is based on counting the number of observations that fall above/below the expected median. It’s easiest to understand with an example.

Q:
A manufacturer claims that the median lifespan of a certain brand of light bulbs is 1000 hours. A consumer protection group suspects that the actual median lifespan might be different. They collect a random sample of 10 bulbs from store shelves and record their lifespans (in hours):

# Sample data: lifespans of 10 light bulbs
lifespans <- c(980, 1020, 1010, 970, 1050, 1030, 1000, 950, 990, 1040)
lifespans

##  [1]  980 1020 1010  970 1050 1030 1000  950  990 1040

Determine whether there is evidence that the median lifespan differs from 1000 hours.

A:
Here we have a small sample of only 10 observations that clearly does not come from a normal distribution. Our goal is to see if the median lifespan is different from 1000.

• What’s the first thing we do?
  • plot the data?
  • see if the median is close to 1000

plot(hist(lifespans))

samp.med = median(lifespans)
samp.med

## [1] 1005

With such a small sample, the histogram doesn’t provide much useful information to us.

The median of 1005 is close to 1000, but also doesn’t give us much info either.

The trick to solving this problem is to think about the null hypothesis scenario, i.e. how the data would behave if the population had a median of 1000.

• If the pop median was 1000 and you sampled 10 observations
  • how many would you expect would be above 1000?
  • how many would you expect would be below 1000?

If the population median is \(\theta_0\) we expect half of the sample observations to fall below and half above. This is our basis of the statistical hypothesis test.

We can use the same framework we used for our categorical comparisons.

• If you flip \(n\) fair coins, what are the chances of seeing more than \(k\) heads or more than \(k\) tails (i.e. less than \(n-k\) heads)?
  • two sided binomial test with prob = 0.5

binom.test(k, n, p=0.5)

• If you sample \(n\) observations from a distribution with median \(\theta_0\) what are the chances that more than \(k\) have a value higher than \(\theta_0\) or more than \(k\) have a value lower than \(\theta_0\)?

# Hypothesized median under H0
median_h0 <- 1000

# Calculate the difference from hypothesized median
differences <- lifespans - median_h0

# Determine the signs of differences
# whether the observation falls above/below
signs <- sign(differences)

# Remove zeros (cases exactly equal to the hypothesized median)
signs_no_ties <- signs[signs != 0]

# Count positive signs (lifespans > 1000)
n_positive <- sum(signs_no_ties > 0)

# Total number of non-tied observations
n_total <- length(signs_no_ties)

# Perform a two-sided binomial test
test_result <- binom.test(n_positive, n_total, p = 0.5, alternative = "two.sided")

# Print results
print(test_result)

## 
##  Exact binomial test
## 
## data:  n_positive and n_total
## number of successes = 5, number of trials = 9, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2120085 0.8630043
## sample estimates:
## probability of success 
##              0.5555556

Notes on the Sign Test

• \(H_0\): median = \(M_0\)
• \(H_1\): median \(\ne M_0\)
• You can use it on small samples
• You don’t need to know the distribution of your data
• there should be no ties
  • take them out before running the test
• data needs to be at least ordinal (have to be able to rank the data points)
• paired samples can be tested by looking at their matched differences.

Two sample Komolgorov-Smirnov Test (Continuous data)

This is a test to check if two independent samples have the same distribution. We’re not trying to figure out what that distribution is, we’re just checking if it’s reasonable to say the two samples were drawn from it.

• How do we approach the problem of comparing two independent samples when we don’t know their distributions?
  • with large sample size (\(n \gtrsim 30\))
    • use central limit theorem
      • t-test (compare means)
  • with small sample size
    • can’t use central limit theorem
• What can we test?
  • do the two samples come from the same distribution?
    • Is there a way to figure this out without needing to know the the mathematical form of the distributions?
• What’s an intuitive way to approach this?
  • plot the histograms and see if they align
    • not very good with small samples
• Instead of histograms, we can use a tool called the empirical cumulative distribution function (eCDF) to construct a CDF from a sample.
  • it returns a function that for any value x, it gives the proportion of the sample that appears less than or equal to x.
  • gives a idea of the distributional shape of a sample
    • a step plot that increase by 1/n at each sample point
    • flat regions are where there are no new data points
    • always starts at 0, ends at 1
    • doesn’t use binning like a histogram.
  • samples from similar distributions will have a similar eCDF shape.

The Komolgovorv Smirnov test computes the empirical cumulative distribution of each sample, and calculates the largest vertical distance \(D\) between the two. Under the null hypothesis that the two samples come from the same distribution, this distance should be small. Under the null, this statistic \(D\) follows the Komolgorov-Smirnov distribution and we can calculate p-value probabilities just like all the other hypothesis tests we’ve gone over so far.

set.seed(1)
# Sample A: from standard normal distribution
sample_A <- rnorm(100, mean = 0, sd = 1)

# Sample B: from uniform distribution
sample_B <- runif(100, min = -2, max = 2)


# Set up a two-panel plot: histogram and ECDF
par(mfrow = c(1, 2))  # 1 row, 2 columns

# --- Plot 1: Histograms ---
hist(sample_A, col = rgb(1, 0, 0, 0.5), breaks = 15,
     xlim = range(c(sample_A, sample_B)), main = "Histograms",
     xlab = "Value", freq = FALSE)
hist(sample_B, col = rgb(0, 0, 1, 0.5), breaks = 15,
     add = TRUE, freq = FALSE)
legend("topright", legend = c("Sample A: Normal(0,1)", "Sample B: Uniform(-2,2)"),
       fill = c(rgb(1, 0, 0, 0.5), rgb(0, 0, 1, 0.5)), border = NA)

# --- Plot 2: ECDFs ---
plot(ecdf(sample_A), verticals = TRUE, col = "red", lwd = 2,
     main = "ECDFs and KS Test", xlab = "Value", ylab = "F(x)")
lines(ecdf(sample_B), verticals = TRUE, col = "blue", lwd = 2)
legend("bottomright", legend = c("Sample A: Normal", "Sample B: Uniform"),
       col = c("red", "blue"), lwd = 2)

# Perform the two-sample Kolmogorov-Smirnov test
ks_result <- ks.test(sample_A, sample_B)

# Print the KS test result
print(ks_result)

## 
##  Asymptotic two-sample Kolmogorov-Smirnov test
## 
## data:  sample_A and sample_B
## D = 0.29, p-value = 0.0004453
## alternative hypothesis: two-sided

Notes on the KS Test

• two sample test
  • Tests whether two samples come from the same distribution.
    • \(H_0\): they are the same
    • \(H_1\): they are the different
  • Makes no assumptions about the specific form of the distribution.
  • Assumes both samples are from continuous distributions with no ties.
• one sample test
  • Tests whether a sample comes from a specific distribution, e.g., normal.
  • Assumes the parameters of the distribution are known and fixed.
    • you supply these parameters to the test.

What about data that doesn’t seem to fit any of these frameworks?

What about data that doesn’t seem to fit any of these frameworks?

idea #2 (resampling)

If enumerating every possibility is not feasible, we can still get an idea of our signal by enumerating what we can.

Let’s test out this out on our problem.

• What is our signal?
  • The contamination rank sum of the 2 family households with birth defects.

Under the null hypothesis there’s no association between contamination and defect so reordering the labels should not make a difference. We want to build a null distribution by simulation. What does our signal look like with the ranks are randomly assigned?

• repeat many times
  • shuffle the contamination levels across our samples.
  • compute the rank sum of the 2 patients with defects.

#repeat this many times
K = 10000
null.ranks = sapply(1:K, function(i){
    #since we're only interested in the 2 families with birth defects
    #we don't need to resample all 20 to find the rank-sum
    defect.ranks = sample(1:20, size=2, replace=F)
    sum(defect.ranks)
})

head(null.ranks, 30)

##  [1] 29 21 10 23 22 21 12 30  9 26 28 12 18  9 27 29 32 19 18 27 29 11 16  5 21
## [26] 24 19 16 28 14

hist(null.ranks, breaks=20)

#how many of these null ranks as less than our observed signal?
n.extreme = sum(null.ranks <= (1+4))

resampling.pval = n.extreme / K
resampling.pval

## [1] 0.02

• Notes
  • the idea is very flexible
    • if you know the distribution of the point estimate
      • you can calculate p-values mathematically
    • if you don’t know the distribution of the point estimate
      • you can build a null distribution, then count the number of extremes to get a p-value
  • not restricted to rank sum comparisons, any function of your sample will work
  • you should have a computationally efficient way to measure signal so the null distribution can be built in a timely manner.
  • doesn’t have any distributional assumptions.

back to idea #1

Enumerate all possible values, and see how many of them are more extreme than what we’ve observed in our sample.

This test already has a name

If you can enumerate all the extreme signal values, then you can calculate the pvalue directly.

• In our example there’s 190 ways of randomly choosing 2 families from 20, and only 4 cases more extreme than what we’ve observed in our sample.
  • (1,2), (1,3), (1,4), (2,3)

The p-value is:

4/choose(20,2)

## [1] 0.02105263

Or, you can work with the ranks and use the builtin procedures. If you know the ranks of one group, then you can take the set difference to get the ranks of the other group. With two sets of ranks, you can now run the Mann-Whitney Rank-Sum Test.

#ranks in group 1
ranks1 = c(1,4)
#ranks in group 2  (what's left over)
ranks2 = (1:20)[-ranks1] 

ranks1

## [1] 1 4

ranks2

##  [1]  2  3  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20

wilcox.test(ranks1, ranks2, alternative="less")

## 
##  Wilcoxon rank sum exact test
## 
## data:  ranks1 and ranks2
## W = 2, p-value = 0.02105
## alternative hypothesis: true location shift is less than 0

• Notes on the Mann-Whitney rank sum testc
  • it’s also known as the two sample Wilcoxon rank sum test.
    • in R, both the 1 and 2 sample versions are called with wilcox.test
  • instead of working with the samples values, you work with their ranks.
  • 1 sample: Assuming a symmetric distribution
    • \(H_0\): median = \(M_0\)
    • \(H_1\): median \(\ne\) \(M_0\)
  • 2 samples: Assuming the distributions have the same shape
    • \(H_0\): there is no difference in medians
    • \(H_1\): they have different medians
  • ties provide no useful information and will be removed from the analysis
  • When n is small, it uses the exact probabilities
  • When n is large, it uses a normal approximation

Dealing with multiple comparisons

Dealing with multiple comparisons

What is type I error?

• It’s when you falsely reject the null hypothesis (you see signal that isn’t there)
  • Testing positive for a disease when you’re healthy
  • Jury convicts a person when they’re innocent.

what is significance level (\(\alpha\))

• It’s the probability of getting a type I error in a hypothesis test.
  • a threshold on the amount of false positives you’re willing to put up with.

When you’re running many hypothesis tests, the probability of making a type I error in at least one of them is lot higher than if you just look at a single test.

Dealing with multiple comparisons

What is p-Hacking?

• P-hacking refers to manipulating the analysis until statistically significant results (typically p < 0.05) are found. It can include:

  • Running many tests and only reporting the significant ones

  • Cherry-picking data or variables

  • Trying different statistical models or transformations until p-values become significant

Without multiple testing correction, p-hacking can easily produce misleading results that appear valid.

How do you address false positives when running multiple tests?

Multiple comparisons in R

The three methods above are already implemented in R and can be accessed through the ‘p.adjust’ function. You input the results of your tests as a list of p-values, specify the method you want to use, and it returns the new list of “adjusted p-values” that have been corrected for the multiple comparisons.

An Example

To illustrate the need to adjusting for multiple comparisons, we’re going to construct a dataset that simulates the results of many experimental runs we’re we’re trying to see if there’s a difference in the means between two groups. In most of the runs(n=1000), the two groups were generated from the same distribution, but in some of the runs(n=100) the two group come from different distributions.

We’ll see what happens when we don’t adjust for the multiple comparisons over all(n=1100) tests.

1. Set up the data

set.seed(1)
#create a fake dataset that has 1000 variables with 'no signal'.  
x.nosig <- as.data.frame(matrix(rnorm(100 * 1000), nrow = 100))  # 100 obs, 200 variables

#create another fake dataset that has 100 variables with 'signal'
x.sig = rbind(
  matrix(rnorm(100*50, mean=1, sd=1), nrow=50),
  matrix(rnorm(100*50, mean=1.5, sd=1), nrow=50)
)

#make a rectangular matrix for the data.  
#Variables are along the rows, subjects along the columns
x = t(cbind(x.nosig, x.sig))
dim(x)

## [1] 1100  100

#The condition to split the subjects by
group = rep(c("A", "B"), each = 50)

#make some fake names for the variables
rownames(x) = paste0("gene", 1:nrow(x))
#make some fake names for the subjects
colnames(x) = paste0(group, 1:50)


library("DT")
datatable(head(x))

datatable(tail(x))

2. Run the tests

# Compare means across variables
# loop over all the rows
#   run t-test split by group and return the p-value
pvals <- sapply(1:nrow(x), function(i){
  t.test(x[i,] ~ group)$p.value
})

3. Adjust for multiple comparisons

# Adjust p-values to account for the multiple comparisons over all the variables
padj.bonf <- p.adjust(pvals, method = "bonferroni")
padj.fdr <- p.adjust(pvals, method = "fdr")

#how many tests pass the (< 0.05) criteria?

#unadjusted
sum(pvals < 0.05)

## [1] 125

#bonferroni
sum(padj.bonf < 0.05)

## [1] 7

#FDR 
sum(padj.fdr < 0.05)

## [1] 27

#Which variables were significant at FDR < 0.05?
#get the indices of the rows that meet the threshold criteria 
ix = padj.fdr < 0.05

datatable(x[ix,])

Recap on the tests

Recap on the tests (One sample)

library("DT")
datatable(mtcars)

t.test(x=mtcars$mpg, mu=20, alternative="greater")

## 
##  One Sample t-test
## 
## data:  mtcars$mpg
## t = 0.08506, df = 31, p-value = 0.4664
## alternative hypothesis: true mean is greater than 20
## 95 percent confidence interval:
##  18.28418      Inf
## sample estimates:
## mean of x 
##  20.09062

wilcox.test(x=mtcars$mpg, mu=20, alternative="greater")

## Warning in wilcox.test.default(x = mtcars$mpg, mu = 20, alternative =
## "greater"): cannot compute exact p-value with ties

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  mtcars$mpg
## V = 249, p-value = 0.614
## alternative hypothesis: true location is greater than 20

#first make sure to remove the observations that
#provide no useful information.  These would be the 
#cases where the mpg is exactly equal to 20.
#there are no cases like this in this example.
binom.test(x=sum(mtcars$mpg>20), length(mtcars$mpg), alternative="greater")

## 
##  Exact binomial test
## 
## data:  sum(mtcars$mpg > 20) and length(mtcars$mpg)
## number of successes = 14, number of trials = 32, p-value = 0.8115
## alternative hypothesis: true probability of success is greater than 0.5
## 95 percent confidence interval:
##  0.2872749 1.0000000
## sample estimates:
## probability of success 
##                 0.4375

Recap on the tests (Paired sample)

library("DT")
datatable(sleep)

#first check that the order is correct
#If not correct, you have to rearrange the rows to align the paired observations
sleep$ID[sleep$group==1] ==  sleep$ID[sleep$group==2]

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

x1 = sleep$extra[sleep$group==1]
x2 = sleep$extra[sleep$group==2]
t.test(x1, x2, paired = TRUE)

## 
##  Paired t-test
## 
## data:  x1 and x2
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean difference 
##           -1.58

#another way to set it up is to look at the paired differences
#and treat it like a 1 sample comparison to a mean of 0.
t.test(x1-x2, paired=F)

## 
##  One Sample t-test
## 
## data:  x1 - x2
## t = -4.0621, df = 9, p-value = 0.002833
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.4598858 -0.7001142
## sample estimates:
## mean of x 
##     -1.58

wilcox.test(x1, x2, paired = TRUE)

## Warning in wilcox.test.default(x1, x2, paired = TRUE): cannot compute exact
## p-value with ties

## Warning in wilcox.test.default(x1, x2, paired = TRUE): cannot compute exact
## p-value with zeroes

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  x1 and x2
## V = 0, p-value = 0.009091
## alternative hypothesis: true location shift is not equal to 0

#similarly, you can look at the paired differences and treat it as 1 sample comparison
wilcox.test(x1-x2)

## Warning in wilcox.test.default(x1 - x2): cannot compute exact p-value with ties

## Warning in wilcox.test.default(x1 - x2): cannot compute exact p-value with
## zeroes

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  x1 - x2
## V = 0, p-value = 0.009091
## alternative hypothesis: true location is not equal to 0

#make a new variable that has the difference of the two drugs
extra.diff = x1 - x2


#be sure to take out non-informative rows 
extra.diff = extra.diff[extra.diff != 0]

#now check if the signs of the differences corresponds to an even coin flip
binom.test(x=sum(extra.diff > 0), n=length(extra.diff), alternative="two.sided")

## 
##  Exact binomial test
## 
## data:  sum(extra.diff > 0) and length(extra.diff)
## number of successes = 0, number of trials = 9, p-value = 0.003906
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.0000000 0.3362671
## sample estimates:
## probability of success 
##                      0

Recap on the tests (Two independant samples)

datatable(ToothGrowth)

#make a smaller dataset that only has the dose 1.0 samples
tooth.dose1 = ToothGrowth[ToothGrowth$dose == 1.0,]
datatable(tooth.dose1)

library("car")

## Loading required package: carData

leveneTest(len~supp, data = tooth.dose1)

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  1.6722 0.2123
##       18

t.test(len~supp, data = tooth.dose1, var.equal=T)

## 
##  Two Sample t-test
## 
## data:  len by supp
## t = 4.0328, df = 18, p-value = 0.0007807
## alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
## 95 percent confidence interval:
##  2.840692 9.019308
## sample estimates:
## mean in group OJ mean in group VC 
##            22.70            16.77

t.test(len~supp, data = tooth.dose1, var.equal=F)

## 
##  Welch Two Sample t-test
## 
## data:  len by supp
## t = 4.0328, df = 15.358, p-value = 0.001038
## alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0
## 95 percent confidence interval:
##  2.802148 9.057852
## sample estimates:
## mean in group OJ mean in group VC 
##            22.70            16.77

wilcox.test(len~supp, data = tooth.dose1)

## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  len by supp
## W = 88.5, p-value = 0.00403
## alternative hypothesis: true location shift is not equal to 0

ks.test(len~supp, data = tooth.dose1)

## 
##  Exact two-sample Kolmogorov-Smirnov test
## 
## data:  len by supp
## D = 0.8, p-value = 0.00158
## alternative hypothesis: two-sided

Notes

• In General, the stronger the assumptions, the more tailored the test is to the data and the more power you’ll have to detect signal.

  • Even though you can use a sign test for any one sample dataset, you’ll get better results from a t-test, so long as all it’s assumptions are met.

• Plan out a workflow of which tests to use before running the analysis. Failure to do so could be considered a form a p-hacking, where multiple testing procedures are tried, and only the most significant finding is reported.

• Make sure you’re plotting/summarizing your data, to ensure results are being interpreted correctly, especially the directionality of the signal. It’s very easy to make coding mistakes than can completely reverse the meaning of a result.

• Mistakes happen, so it’s very important to document the results as well as the code used to generate them.