11 Non-parametric test

คือ การวิเคราะห์ทางสถิติที่เราไม่ทราบลักษณะการกระจายตัวของข้อมูลอย่างชัดเจน ข้อดีคือมีความยืดหยุ่นกว่า แต่มีความแม่นยำน้อยกว่า

11.1 Proportion test

คือ การทดสอบทางสถิติเพื่อทำการเปรียบเทียบอัตราส่วนของจำนวน โดยจะใช้กับข้อมูลประเภทจำนวนนับของตัวแปรจัดประเภท (Nominal variable) เช่น จำนวนคน จำนวนเซลล์ เป็นต้น

ข้อสังเกต การทดลองที่มีการทำ Technical replicate แล้วหาค่าเฉลี่ยนั้น ตัวข้อมูลยังเป็น จำนวนนับ การจะหาความต่างค่าเฉลี่ยโดยใช้ t.test นั้น จำนวน Sample ควรจะมากพอตาม Central limit theorem นอกเหนือจากนั้นควรใช้ Proportional test หรือ ควรใช้ Generalized linear model ประเภทอื่นมากกว่า

11.1.1 Chi-square test

คือ การทดสอบความต่างของอัตราส่วนโดยใช้การประมาณการของ ค่าที่คาดหวัง (Expected value) และ ค่าที่สังเกตได้จริง (Observed value)

\[ \chi^{2} = \sum_{i=1}^{n}\frac{(O-E)^{2}}{E} \]

11.1.1.1 Chi-square goodness of fit

เป็นการเปรียบเทียบว่า Observed value นั้นมาจากประชากรทางทฤษฎีหรือไม่ โดย Expected value นั้นคำนวนจาก

\[ E_{i} = CDF_{i}(Y_{u})-CDF_{i}(Y_{l}) \]

ยกตัวอย่างการทอดลูกเต๋า ซึ่งมีโอกาสการเกิดทุกหน้า = \(1/6\)

\(H_{0}\) การกระจายตัวของการทอดลูกเต๋านี้เท่ากับการกระจายทางทฤษฎี
\(H_{a}\) การกระจายตัวของการทอดลูกเต๋านี้ไม่เท่ากับการกระจายทางทฤษฎี

set.seed(123)
dice <- sample(6, 1000, replace = TRUE) # โยนลูกเต๋า 1000 ครั้ง
chisq_dice <- chisq.test(table(dice))
chisq_dice

## 
##  Chi-squared test for given probabilities
## 
## data:  table(dice)
## X-squared = 1.436, df = 5, p-value = 0.9203

ลองคำนวณเองตามสูตรขั้นต้น

chisq_gof_df <- data.frame(O = chisq_dice$observed, 
                           E.Freq = chisq_dice$expected) |> 
  mutate(X2 = ((`O.Freq` - `E.Freq`)^2)/`E.Freq`)

chisq_gof_df

chisq_gof <- sum(chisq_gof_df$X2)
chisq_gof

## [1] 1.436

pchisq(chisq_gof, df = 5, lower.tail = FALSE)

## [1] 0.9203351

จะเห็นว่าความแตกต่างระหว่าง Observed และ Expected นั้น อยู่ในพิสัยของของ \(H_{0}\)

แต่ถ้าลูกเต๋านั้นเป็นลูกเต๋าถ่วงน้ำหนัก เมื่อวิเคราะห์ Chi-square ตามการกระจายตัวของลูกเต๋าทั่วไป จะแตกต่างอย่างมีนัยสำคัญ

weight_dice <- sample(6,1000, replace = TRUE, prob = c(3,2,1,1,1,1)/9)
chisq_weight_norm <- chisq.test(table(weight_dice))
chisq_weight_norm

## 
##  Chi-squared test for given probabilities
## 
## data:  table(weight_dice)
## X-squared = 248.55, df = 5, p-value < 2.2e-16

data.frame(O = chisq_weight_norm$observed, 
           E.Freq = chisq_weight_norm$expected) |> 
     mutate(X2 = ((`O.Freq` - `E.Freq`)^2)/`E.Freq` )

แต่ถ้าวิเคราะห์เทียบกับการกระจายตัวเดียวกับลูกเต๋าถ่วงน้ำหนัก จะไม่แตกต่างกัน

chisq_weight_weight <- chisq.test(table(weight_dice), p = c(3,2,1,1,1,1)/9)
chisq_weight_weight

## 
##  Chi-squared test for given probabilities
## 
## data:  table(weight_dice)
## X-squared = 2.9135, df = 5, p-value = 0.7133

data.frame(O = chisq_weight_weight$observed, 
           E.Freq = chisq_weight_weight$expected) |> 
     mutate(X2 = ((`O.Freq` - `E.Freq`)^2)/`E.Freq` )

11.1.1.2 Chi-square test of independence

เป็นการเปรียบเทียบข้อมูลจำนวนสองกลุ่มขึ้นไปว่า มีความสัมพันธ์ที่ทำให้การกระจายตัวของข้อมูลเปลี่ยนไปจากปกติหรือไม่

\(H_{0}\): ข้อมูลทั้ง 2+ กลุ่มนั้นไม่มีความสัมพันธ์ต่อกัน

\(H_{a}\): ข้อมูลทั้ง 2+ กลุ่มนั้นมีความสัมพันธ์ต่อกัน

โดยการวิเคราะห์นั้นจะใช้กับข้อมูลความถี่แบบ \(m \times n\) โดย Expected event นั้นคิดจาก

Condition/Group	Group 1	Group 2	Total
Cond 1	A	B	A + B
Cond 2	C	D	C + D
Total	A + C	B + D	A + B + C + D

\[ E_{i,j} = P(Cond_{i,j}) \times P(Group_{i,j}) \times \text{total counts} \] \[ E = \frac{\text{Row total} \times \text{Column total}}{\text{Total sample size}} \]

อย่างเช่น Expected event สำหรับช่อง \(A\) คือ \[ \frac{(A + B) \times (A + C)}{(A + B + C + D)^{2}}(A+B+C+D) \]

ซึ่งอยู่ภายใต้ \(df = (n_{row} - 1)\times(n_{col} -1)\)

ยกตัวอย่างว่าอยากทราบว่าเพศมีผลต่ออัตราการตายในมะเร็งปอดหรือไม่

lung_ob

##         status
## sex      Alive Dead
##   Female    37   53
##   Male      26  112

lung_chisq <- chisq.test(lung_ob, correct = FALSE)
lung_chisq

## 
##  Pearson's Chi-squared test
## 
## data:  lung_ob
## X-squared = 13.511, df = 1, p-value = 0.0002371

\(p\) < 0.05 หมายความว่าเพศมีผลต่ออัตราการตายอย่างมีนัยสำคัญ

ลองคำนวณเองตามสูตรขั้นต้น

lung_exp <- expand.grid(rowSums(lung_ob), colSums(lung_ob)) |>
    apply( 1, prod) |> 
  matrix(nrow = 2)/sum(lung_ob) # สร้าง expected table

lung_exp

##          [,1]     [,2]
## [1,] 24.86842 65.13158
## [2,] 38.13158 99.86842

chi_value <- sum((lung_ob - lung_exp)^2/lung_exp)
chi_value # chi-squared

## [1] 13.51117

pchisq(chi_value, df = lung_chisq$parameter, lower.tail = FALSE) # p-value

## [1] 0.000237147

11.1.2 Fisher’s exact test

คือการทดสอบว่าข้อมูลนั้นมีความสัมพันธ์หรือไม่ โดยการเทียบกับความสัมพันธ์แบบสุ่ม การทดสอบนี้จะมีความแม่นยำกว่า Chi-square เนื่องจากเป็นการคำนวณความน่าจะเป็นโดยตรง เมื่อพิจารณาตาราง \(2 \times 2\) ในหลักของการหยิบสุ่ม จะได้ตารางดังนี้

Condition/Group	Group 1	Group 2	Total
Cond 1	\(k\)	\(n-k\)	\(n\)
Cond 2	\(K-k\)	\(N-K-n-k\)	\(N-n\)
Total	\(K\)	\(N-K\)	\(N\)

ซึ่งความน่าจะเป็นที่จะหยิบสุ่มได้สถานการณ์ที่ต้องการจำนวน \(k\) ครั้ง คือ

\[ p_{k} = \frac{{K \choose k}{N-K \choose n-k}}{{N \choose n}} \]

ใน Fisher’s exact test \(p\) ที่คำนวณได้คือ ผลรวมของโอกาสที่จะเกิดสถานการณ์ \(p_{k}\) และโอกาสที่เกิดขึ้นสถานการณ์ที่ยากกว่า สถานการณ์ \(p_{k}\)

lung_fish <- fisher.test(lung_ob, alternative = "greater")
lung_fish

## 
##  Fisher's Exact Test for Count Data
## 
## data:  lung_ob
## p-value = 0.0002318
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
##  1.740038      Inf
## sample estimates:
## odds ratio 
##   2.991585

\(p\) < 0.05 หมายความว่าเพศมีผลต่ออัตราการตายอย่างมีนัยสำคัญ (โอกาสที่จะเกิดสถานการณ์นี้ ยากกว่าสถานการณ์หยิบสุ่มแบบไม่คืนโดยทั่วไป) สังเกตว่า \(p\) จะมากกว่า Chi-square เนื่องจากการทดสอบนี้มีความ Conservative กว่า

ลองคำนวณเองตามสูตรขั้นต้น หลักการคือ การหาผลรวมของโอกาสการสุ่มที่จะเกิดสถานการณ์ที่ต้องการ และ สถานการณ์ที่เกิดยากกว่านั้น

lung_ob_sum <- addmargins(lung_ob)
lung_ob_sum

##         status
## sex      Alive Dead Sum
##   Female    37   53  90
##   Male      26  112 138
##   Sum       63  165 228

## Margins
nm <- lung_ob_sum[1,3]
N_nm <- lung_ob_sum[2,3]
Km <- lung_ob_sum[3,1]
N_Km <- lung_ob_sum[3,2]

## Position in the matrix
c1 <- 0:nm 
c2 <- nm-c1
c3 <- 0:N_nm
c4 <- N_nm-c3

## Build a permutation grid with fixed margin (n, N-n, K, N-K)
filtered_permuted_grid <- expand_grid(c1,c2,c3,c4) |> 
  mutate(n = c1+c2) |> 
  mutate(N_n = c3+c4) |> 
  mutate(K = c1+c3) |> 
  mutate(N_K = c2+c4) |>
  filter(n == nm & N_n ==  N_nm & K == Km & N_K == N_Km)

## Compute p for each condition
permuted_grid_p <- filtered_permuted_grid |> 
  mutate(p = (choose(K, c1)*choose(N_K, c2)/(choose(n+N_n, n))) )
permuted_grid_p

## Fisher's exact pvalue
pval_greater <- permuted_grid_p |> filter(c1 >= 37 ) |> pull(p) |> sum()
pval_greater

## [1] 0.000231783

## or
phyper(lung_ob[1]-1, Km, N_Km, nm, lower.tail=FALSE)

## [1] 0.000231783

Note: การใช้ lung_ob[1] - 1 นั้นมีที่มาจากว่า เมื่อใช้ lower.tail = FALSE จะคำนวณโอกาสที่ ได้ \(P[X > x]\) ซึ่งเราต้องการ \(P[X \geq x]\) จึงต้องเริ่มคำนวณตั้งแต่สถานการณ์ก่อนหน้า 1 สถานการณ์

สำหรับการทดสอบแบบ two.sided คือการหาสถานการณ์สุดโต่งทั้งสองด้าน โดยการคำนวณ \(p_{con}\) ของสถานการณ์ที่ต้องการก่อน และรวมผลลัพธ์ของสถานการณ์ทั้งหมดที่ \(p \leq p_{con}\)

pbase <- permuted_grid_p |> filter(c1 ==37) |> pull(p)
pval_twoside <- permuted_grid_p |> filter(p <= pbase) |> pull(p) |> sum()
pval_twoside

## [1] 0.0004349142

fisher.test(lung_ob, alternative = "two.sided")$p.value ## same

## [1] 0.0004349142

11.2 Rank test

เป็นการหาความต่างของลำดับ แทนที่จะหาความแตกต่างของ Mean/Variance ในภาวะที่ไม่ทราบการกระจายตัวของข้อมูลที่ชัดเจน

11.2.1 Test for normality

กลับมาที่ตัวอย่าง iris ดูการกระจายของข้อมูล

  ggplot(long_df, aes(x = cm, fill = Metrics)) + geom_histogram() +
  facet_grid(Metrics~Species, scales = "free") +
  theme_bw()

พิจารณาแล้ว Petal.Width ไม่น่าจะใช่ Normal distribution จะทำการทดสอบต่อโดย shapiro.test() ซึ่งเป็นการทดสอบการกระจายตัวของข้อมูลว่าหลุดออกจาก Normal distribution หรือไม่

long_df |> 
  group_by(Metrics, Species) |> 
  summarise(normality = round(shapiro.test(cm)$p.value, 4))

จะเห็นว่า Petal.Width มี \(p\) < 0.05 พอสมควร (ไม่เป็น Normal distribution) จึงสมควรใช้ Non-parametric test

11.2.2 Wilcoxon’s rank sum test (Mann-Whitney U test)

คือ การคำนวณว่ากลุ่มตัวอย่าง 2 กลุ่มนั้นมาจากประชากรกลุ่มเดียวกันหรือไม่ จากการพิจารณาผลรวมของลำดับข้อมูลทั้งสองกลุ่ม

ลักษณะการใช้คล้าย Independent t-test สำหรับ Non-normal distribution

\[ W_{j} = n_{1}n_{2} + \frac{n_{j}(n_{j}+1)}{2} - R_{n} \]

\[ W = min(W_{1}, W_{2}) \]

iris_pw <- df |> 
  select(Petal.Width, Species) |> 
  filter(Species != "virginica" )

iris_wx <- wilcox.test(Petal.Width~Species, data = iris_pw)

iris_wx

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Petal.Width by Species
## W = 0, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

iris_wx$p.value

## [1] 2.284669e-18

ลองคำนวณเอง หลักการคือ

จัดอันดับของข้อมูลโดยเรียงจาก น้อยไปมาก และให้อันดับเป็นตัวเลข (อันดับที่เท่ากัน ให้เป็นค่าเฉลี่ยของลำดับนั้น เช่น อันดับ 2, 3, 4 ที่มีค่าเท่ากัน ให้อันดับเป็น (2+3+4)/3 = 5 ทุกตัว)
คำนวณค่า \(W\) โดยแยกกลุ่ม 1 และ กลุ่ม 2 และเลือกค่า \(W\) ที่น้อยที่สุด

iris_pw_rank <- iris_pw |> 
  arrange(`Petal.Width`) |> 
  mutate(Rank = rank(`Petal.Width`, ties.method = "average")) # rank all values

setosa <- split(iris_pw_rank, iris_pw$Species)$setosa
versicolor <- split(iris_pw_rank, iris_pw$Species)$versicolor

all_combn <- nrow(setosa)*nrow(versicolor)
all_combn

## [1] 2500

setosa_rank_sum <- all_combn + (nrow(setosa)*(nrow(setosa)+1))/2 - sum(setosa$Rank)
setosa_rank_sum

## [1] 2500

versicolor_rank_sum <- all_combn + (nrow(versicolor)*(nrow(versicolor)+1))/2 - sum(versicolor$Rank)
versicolor_rank_sum

## [1] 0

W <- min(setosa_rank_sum,versicolor_rank_sum)
W

## [1] 0

## Normal approximation
mW <- all_combn/2
sdW <- sqrt(all_combn*(nrow(setosa)+nrow(versicolor)+1)/12)
z = (W-mW)/sdW
pval <- 2*pnorm(-abs(z)) 
pval # not exactly equal due to tie adjustment in wilcox.test()

## [1] 6.856641e-18

Note: เมื่อจำนวนตัวอย่าง >50 ค่า \(W\) จะเข้าสู่ Normal distribution ซึ่งสามารถเปลี่ยนเป็นค่า \(Z\) ได้ ส่งผลให้การคำนวณ \(p\)-value จะแม่นยำขึ้น

\[ Z = \frac{W-m_{w}}{s_{w}} \]

\[ m_{w} = \frac{n_{1}n_{2}}{2} = \text{mean of W} \]

\[ s_{w} = \sqrt{\frac{n_{1}n_{2}(n_{1}+n_{2}+1)}{12}} \]

11.2.3 Wilcoxon’s signed-rank test

คือ การคำนวณว่ากลุ่มตัวอย่าง 2 กลุ่มนั้นมาจากประชากรกลุ่มเดียวกันหรือไม่ จากการคำนวณลำดับข้อมูลของทั้งสองกลุ่ม

ลักษณะการใช้คล้าย Paired t-test สำหรับ Non-normal distribution

\[ V = \sum_{i=1}^{N_{r}}[sgn(x_{2,i} - x_{1,i}) \cdot R_{i}] \]

\[ sgn(x) =\begin{cases} -1 \quad \text{if} \, x < 1, \\0 \quad \ \ \ \text{if} \, x = 0,\\1 \quad \ \ \ \text{if} \, x > 0\end{cases} \]

\[ V = min(V_{-}, V_{+}) \]

สมมติการวัด Wound healing assay (วัดการเคลื่อนที่ของเซลล์) ก่อนและหลัง Treat ยา

set.seed(123)

treat <- data.frame(sample = 1:20,
  Before = runif(20, min = 400, max = 500),
           After = detectnorm::rnonnorm(20, mean = 400, sd = 30, skew = 10, kurt = 5)$dat)

head(treat, 10)

wilcox.test(treat$Before, treat$After, paired = TRUE, exact = TRUE)

## 
##  Wilcoxon signed rank exact test
## 
## data:  treat$Before and treat$After
## V = 169, p-value = 0.01531
## alternative hypothesis: true location shift is not equal to 0

ลองคำนวณเอง หลักการคือ

คำนวณความต่างของก่อนและหลัง (หรือคู่เทียบ)
จัดอันดับของข้อมูลโดยเรียงจาก น้อยไปมาก และให้อันดับเป็นตัวเลข (อันดับที่เท่ากัน ให้เป็นค่าเฉลี่ยของลำดับนั้น เช่น อันดับ 2, 3, 4 เท่ากัน ให้อันดับเป็น (2+3+4)/3 = 5 ทุกตัวแปร)
คำนวณค่า \(V\) โดยแยกเครื่องหมาย + และ - และเลือกค่า \(V\) ที่น้อยที่สุด

treat_rank <- treat |> mutate(Diff = Before-After) |> 
                mutate(Absdiff = abs(Diff)) |> 
                mutate(Rank = rank(Absdiff, ties.method = "average")) 
treat_rank

sign_rank <- treat_rank |> group_by(sign(Diff)) |> 
              summarize(rank_sum = sum(Rank))

V <- min(sign_rank$rank_sum)
V

## [1] 41

pval <- psignrank(V, 20,20)*2
pval

## [1] 0.01531219

11.2.4 Kruskal-Wallis test

คือ การเปรียบเทียบค่าเฉลี่ยของลำดับว่ามาจากประชากรเดียวกันหรือไม่ ลักษณะการใช้เช่นเดียวกับ ANOVA

\[ H = (N-1)\frac{\sum^{g}_{i=1}n_{i}(\bar{r_{i}}-\bar{r})^{2}}{\sum^{g}_{i=1}\sum^{n_{i}}_{j=1}(r_{ij}-\bar{r})^{2}} \]

iris_pw_all <- df |> 
               select(Petal.Width, Species) # 3 groups

iris_kw <- kruskal.test(Petal.Width ~ Species, data = iris_pw_all)
iris_kw$p.value

## [1] 3.261796e-29

ลองคำนวณเอง หลักการการจัดอันดับเหมือน Wilcoxon’s rank sum test

iris_pw_all_ranked <- iris_pw_all |> arrange(`Petal.Width`) |> 
                        mutate(Rank = rank(`Petal.Width`, ties.method = "average"))
iris_pw_all_ranked

average_rank <- mean(iris_pw_all_ranked$Rank) # also = (nrow(iris_pw_all_ranked)+1)/2
average_rank

## [1] 75.5

between_group_var <- iris_pw_all_ranked |>
                  group_by(Species) |> 
                   summarise(rank_var = ((mean(Rank) - average_rank)^2)*n())
between_group_var

within_group_var <- iris_pw_all_ranked |> 
                    mutate(rank_var = (Rank-average_rank)^2)
within_group_var

H <- ((nrow(iris_pw_all_ranked)-1))*sum(between_group_var$rank_var)/sum(within_group_var$rank_var)
H

## [1] 131.1854

pval <- pchisq(H, 2, lower.tail = FALSE)
pval

## [1] 3.261796e-29

11.3 Correlation test

11.3.1 Spearman’s correlation

เป็นการทดสอบความสัมพันธ์ในทิศทางของตัวแปรสองตัวแปรว่าเป็นไปในทางเดียวกันหรือไม่ (Monotonic relationship)

\[ \rho_{xy} =\frac{cov(x,y)}{s(x) s(y)} \]

\[ \rho_{xy} = \frac{\sum^{n}_{i=1}(x_{i}-\bar{x})(y_{i}-\bar{y})}{\sqrt{\sum^{n}_{i=1}(x_{i}-\bar{x})^{2}}\sqrt{\sum^{n}_{i=1}(y_{i}-\bar{y})^{2}}} \]

โดย \(x, y\) คือ ลำดับของข้อมูล (ไม่ใช่ตัวข้อมูลเอง) ส่งผลให้การทดสอบนี้ เป็นการทดสอบการเพิ่มขึ้นของลำดับ ไม่ใช่การเพิ่มขึ้นของข้อมูล ดังนั้น จึงไม่ใช่การทดสอบความสัมพันธ์เป็นเชิงเส้นเหมือน Pearson’s correlation แต่เป็นการทดสอบเพียงว่าข้อมูลไปในทิศทางเดียวันหรือไม่

poly_data <- data.frame(x = seq(-10, 10, length.out = 100)) |> 
  mutate(y = x^9+x+10+rnorm(100, mean = 0, sd =3))

ggplot(poly_data, aes(x = x,y = y)) + geom_point() + 
  geom_smooth(method = "lm", se = FALSE) +
  labs(y = quote(f(x) == x^9+x+10+epsilon)) +
  theme_bw()

cor(poly_data$x, poly_data$y, method = "pearson") # not quite linear

## [1] 0.6870804

cor(poly_data$x, poly_data$y, method = "spearman") # monotonic

## [1] 0.9993399

การทดสอบใน R ใช้ลักษณะ code เดียวกันกับ Pearson’s correlation แต่เปลี่ยน Argument เป็น method = "pearson" และมีข้อควรระวังที่เหมือนกัน

10 Parametric test

12 False discovery