1. Kiểm định tỉ lệ
Giả sử ta muốn khảo sát tỉ lệ giới tính trong một tổng thể (population) là tất cả cư dân của một thành phố. Từ tổng thể, ta thu thập được một mẫu dữ liệu ngẫu nhiên (random sample) gồm 100 người trong đó có 60 nữ và 40 nam. Dựa vào mẫu dữ liệu này, ta kiểm định nghi vấn “tỉ lệ nữ của thành phố là cao hơn nam” với mức ý nghĩa (significance level) 5%.
Các bước tính toán có thể thực hiện trong R như sau:
1 2 3 4 5 6 7 8 9 10 11 |
n <- 100; p0 <- 0.5; p_hat <- 0.6; anpha <- 0.05 (se <- sqrt(p0*(1 - p0)/n)) 0.05 (z <- (p_hat - p0)/se) 2 (p_value <- 1 - pnorm(z)) 0.02275013 (crit_val <- qnorm(1 - anpha)) 1.644854 p_value < anpha; crit_val < z TRUE TRUE |
Hoặc có thể thực hiện nhanh bằng hàm prop.test như sau:
1 2 3 4 5 6 7 8 9 10 11 |
prop.test(60, 100, p = 0.5, conf.level = 1 - anpha, alternative = "greater") 1-sample proportions test with continuity correction data: 60 out of 100, null probability 0.5 X-squared = 3.61, df = 1, p-value = 0.02872 alternative hypothesis: true p is greater than 0.5 95 percent confidence interval: 0.5127842 1.0000000 sample estimates: p 0.6 |
Lưu ý, vì prop.test dùng công thức khác (có hiệu chỉnh liên tục với phân phối Chi-bình phương) nên kết quả tính có khác chút xíu so với cách dùng công thức dựa trên phân phối chuẩn ở trên.
Ta cũng có thể xây dựng khoảng tin cậy trên 1 − α cho e là [a, 1] với
a = ê − |zα|se = 0.6 − 1.6449 × 0.05 = 0.5178
Từ đó, khoảng tin cậy 95% trên cho e là [0.52, 1]. Vì eO = 0.5 không nằm trong khoảng tin cậy này nên ta bác bỏ HO với mức ý nghĩa 1 − 0.95 = 0.05.
Kiểm định bằng mô phỏng randomization:
Mô tả dữ liệu và thống kê tương ứng
1 2 3 4 5 6 7 8 9 |
n <- 100; p0 <- 0.5; p_hat <- 0.6; anpha <- 0.05 nullsample <- c(rep(1, 50), rep(0, 50)) #mẫu dữ liệu tương ứng với H0 stat <- function(data) #thống kê cần tính { return(sum(data)/length(data)) #tỉ lệ mẫu } |
Một mẫu randomization (randomization sample) là mẫu cùng kích thước như mẫu được cho, lấy theo phương pháp lấy mẫu tương ứng với giả sử giả thuyết mặc định đúng. Bằng cách lấy nhiều mẫu randomization và tính thống kê tương ứng, ta có tập kết quả thể hiện sự biến động của thống kê kiểm định theo mẫu (khi giả thuyết mặc định đúng), gọi là phân phối randomization (randomization distribution). Kĩ thuật mô phỏng này được gọi là randomization.
1 2 3 4 5 6 7 |
randomization <- function(B) { return(replicate(B, stat(sample(nullsample, n, replace = TRUE)))) } |
Chạy với số lượng mẫu randomization B đủ lớn thì phâp phối randomization có thể được dùng để “xấp xỉ” phân phối mẫu mặc định. Từ đó ta có thể xấp xỉ p-value cần tính bằng tỉ lệ điểm “khó xảy ra như hoặc hơn” giá trị tương ứng của thống kê trên mẫu được cho. Cũng lưu ý là mỗi lần chạy sẽ có kết quả khác nhau do ngẫu nhiên.
Chẳng hạn thử với B = 10000
1 2 3 4 5 6 7 |
rand_dist <- randomization(10000) (p_value <- sum(rand_dist >= p_hat)/length(rand_dist)) 0.0271 (crit_val <- quantile(rand_dist, 1 - anpha)) 0.58 p_value < anpha; crit_val < p_hat TRUE TRUE |
Vậy p-value xấp xỉ bằng randomization là 0.0271. Kết quả này rất giống với kết quả tính bằng công thức lý thuyết.
Ta cũng có thể xây dựng khoảng tin cậy bằng bootstrap và kiểm định như sau:
1 2 3 4 5 6 7 8 9 |
n <- 100; p0 <- 0.5; p_hat <- 0.6; anpha <- 0.05 sample1 <- c(rep(1, 60), rep(0, 40)) boot_dist <- replicate(10000, mean(sample(sample1, n, replace=TRUE))) (confint <- quantile(boot_dist, c(anpha, 1), names = FALSE)) 0.52 0.78 !(confint[1] <= p0 && p0 <= confint[2]) TRUE |
2. Kiểm định kì vọng
Giả sử ta muốn khảo sát chiều cao trung bình của một tổng thể là tất cả cư dân của một thành phố. Từ tổng thể, ta thu thập được một mẫu dữ liệu ngẫu nhiên gồm chiều cao của 40 người như bảng dưới. Dựa vào mẫu dữ liệu này, ta kiểm định nghi vấn “chiều cao trung bình của cư dân thành phố là 1.60 mét” với mức ý nghĩa 5%.
Bảng dữ liệu chiều cao 40 người (đơn vị mét)
1.56 | 1.47 | 1.59 | 1.65 | 1.62 | 1.78 | 1.69 | 1.49 | 1.92 | 1.55 |
1.65 | 1.52 | 1.65 | 1.60 | 1.71 | 1.48 | 1.69 | 1.65 | 1.59 | 1.74 |
1.70 | 1.61 | 1.58 | 1.65 | 1.75 | 1.65 | 1.46 | 1.53 | 1.59 | 1.62 |
1.60 | 1.55 | 1.57 | 1.46 | 1.57 | 1.63 | 1.46 | 1.68 | 1.53 | 1.48 |
Các bước tính toán có thể thực hiện trong R như sau:
sample1 <- c(1.56, | 1.47, | 1.59, | 1.65, | 1.62, | 1.78, | 1.69, | 1.49, | 1.92, | 1.55, |
1.65, | 1.52, | 1.65, | 1.60, | 1.71, | 1.48, | 1.69, | 1.65, | 1.59, | 1.74, |
1.70, | 1.61, | 1.58, | 1.65, | 1.75, | 1.65, | 1.46, | 1.53, | 1.59, | 1.62, |
1.60, | 1.55, | 1.57, | 1.46, | 1.57, | 1.63, | 1.46, | 1.68, | 1.53, | 1.48) |
n <- length(sample1); mu0 <- 1.6;
(x_bar <- mean(sample1)) |
anpha <- | 0.05 |
1.60675 | ||
(se <- sd(sample1)/sqrt(n)) 0.01557606
(t <- (x_bar – mu0)/se) |
||
0.4333572
(p_value <- 2*(1 – pt(t, df = n – |
1))) |
|
0.6671428
(crit_val <- qt(1 – anpha/2, df = |
n – 1)) |
|
2.022691
p_value < anpha; crit_val < t FALSE FALSE |
Hoặc có thể thực hiện nhanh bằng hàm t.test như sau:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
t.test(sample1, mu = mu0, conf.level = 1 - anpha) One Sample t-test data: sample1 t = 0.43336, df = 39, p-value = 0.6671 alternative hypothesis: true mean is not equal to 1.6 95 percent confidence interval: 1.575244 1.638256 sample estimates: mean of x 1.60675 |
Ta cũng có thể xây dựng khoảng tin cậy giữa 1 − α cho µ là [a, b] với
|
Từ đó, khoảng tin cậy 95% giữa cho µ là [0.52, 1]. Vì µO = 1.60 nằm trong khoảng tin cậy này nên ta không bác bỏ HO với mức ý nghĩa 1 − 0.95 = 0.05.
Kiểm định bằng mô phỏng randomization:
Để giả thuyết HO: µ = µO = 1.60 đúng, ta cần dịch (shift) các dữ liệu trong mẫu được cho một khoảng là µO − x̅ trước khi lấy mẫu randomization.
Với dữ liệu đã được nhập như trên, ta mô tả thêm thống kê tương ứng (ở đây là trung bình mẫu) và hàm thực hiện randomization như phần 1
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 |
nullsample <- sample1 - (x_bar - mu0) #mẫu dữ liệu tương ứng với H0 mean(nullsample) 1.6 stat <- function(data) #thống kê cần tính { return(mean(data)) #trung bình mẫu } randomization <- function(B) { return(replicate(B, stat(sample(nullsample, n, replace = TRUE)))) } |
Chạy thử với B = 10000
1 2 3 4 5 6 7 |
rand_dist <- randomization(10000) (p_value <- sum(abs(rand_dist-mu0) >= abs(x_bar-mu0))/length(rand_dist)) 0.6668 (crit_val <- quantile(rand_dist, 1 - anpha/2)) 1.631 p_value < anpha; abs(crit_val - mu0) < abs(x_bar - mu0) FALSE FALSE |
Vậy p-value xấp xỉ bằng randomization là 0.6668. Kết quả này rất giống với kết quả tính bằng công thức lý thuyết.
Từ mẫu dữ liệu được cho ở trên, ta cũng có thể xây dựng khoảng tin cậy bằng bootstrap và kiểm định như sau:
1 2 3 4 5 |
boot_dist <- replicate(10000, mean(sample(sample1, n, replace = TRUE))) (confint <- quantile(boot_dist, c(anpha/2, 1-anpha/2), names = FALSE)) 1.57700 1.63725 !(confint[1] <= mu0 && mu0 <= confint[2]) FALSE |
Lưu ý là code chạy mô phỏng trong phần 2 này rất giống code chạy mô phỏng trong phần 1, chỉ khác hàm tính thống kê (ở đây là trung bình mẫu so với phần 1 là tỉ lệ mẫu) và hiệu chỉnh cho kiểm định hai phía (so với phần 1 là kiểm định một phía).
3. Kiểm định trung vị
Giả sử cùng tổng thể và mẫu dữ liệu ở phần 2 nhưng ta muốn kiểm định trung vị (median) thay vì trung bình chiều cao. Mặc dù trung bình thường được sử dụng như là con số mô tả trọng tâm của phân phối, nó lại rất nhạy cảm với ngoại lệ (outlier). Nghĩa là, trung bình thay đổi rất nhiều khi mẫu dữ liệu có ngoại lệ. Chẳng hạn trong mẫu dữ liệu ở phần 2 ta có một dữ liệu là 1.92, lớn hơn rất nhiều so với các dữ liệu còn lại. (Trong trường hợp này, đây có thể là chiều cao của một người nước ngoài, đã bị lấy mẫu nhầm:)).
Trong những trường hợp như thế này ta có thể dùng trung vị là một thống kê ít bị ảnh hưởng bởi ngoại lệ. Việc xây dựng khoảng tin cậy cho trung vị phức tạp hơn cho trung bình. Tuy nhiên, với phương pháp tổng quát và đơn giản như randomization, ta dễ dàng làm được điều này.
Với dữ liệu và đoạn mã randomization như phần 2, ta chỉ cần sửa lại thống kê tương ứng để tính trung vị của mẫu thay vì trung bình
1 2 3 4 5 6 7 |
stat <- function(data) #thống kê cần tính { return(median(data)) #trung vị mẫu } |
Chạy thử với B = 10000
1 2 3 4 5 6 7 8 9 |
(x_bar <- mean(sample1)) 1.60675 rand_dist <- randomization(10000) (p_value <- sum(abs(rand_dist-mu0) >= abs(x_bar-mu0))/length(rand_dist)) 0.8533 (crit_val <- quantile(rand_dist, 1 - anpha/2)) 1.64325 p_value < anpha; abs(crit_val - mu0) < abs(x_bar - mu0) FALSE FALSE |
Vậy p-value xấp xỉ bằng randomization là 0.8533 > α = 0.05, do đó ta không bác bỏ giả thuyết HO: median= 1.60 (với đối thuyết H1: median≠ 1.60). Lưu ý là code chạy mô phỏng trong phần này hoàn toàn giống code chạy mô phỏng trong 2 phần trên, chỉ khác hàm tính thống kê (ở đây là trung vị mẫu).
Leave a Reply