Statistical Inference: Khoảng tin cậy T (Confident Interval)

Confident interval

Confident interval

Trong bài viết này, ta sẽ thảo luận về một vài phương pháp thống kê trên tập dữ liệu nhỏ, cụ thể là phân phối T của Student/Gosset và khoảng tin cậy T.

Dẫn nhập

Ở bài viết trước (Asymptotics – tiệm cận), ta đã thảo luận về khoảng tin cậy sử dụng định lý giới hạn trung tâm (Central Limit Theorem (CLT)) và phân phối chuẩn. Cả hai đều cần kích thước mẫu dữ liệu lớn và công thức tính khoảng tin cậy bằng Est +/- qnorm *std error(Est). Trong đó, Est (estimated value) là giá trị ước lượng (ví dụ trung bình mẫu) với độ lỗi chuẩn (standard error). qnorm thể hiện giá trị phân vị (quantile) cụ thể từ phân phối chuẩn.

Ta cũng đề cập đến thống kê Z=(X’-mu)/(sigma/sqrt(n)) có phân phối chuẩn. Phân phối đã được chuẩn hóa (normalized) này trông rất “đẹp” vì ta có mean=0 và variance=1.

Như vậy mean và variance của phân phối chuẩn đã biết và không thay đổi. Bây giờ, ta sẽ định nghĩa thống kê t, trông khá giống Z. Ta định nghĩa t=(X’-mu)/(s/sqrt(n)). Như thống kê Z, t có trọng tâm tại 0. Điểm khác biệt duy nhất là độ lệch chuẩn của quần thể (population), sigma, trong Z được thay thế bởi độ lệch của của tập dữ liệu mẫu (sample) trong t. Do đó, phân phối của thống kê t độc lập với trung bình và phương sai quần thể. Thay vì vậy, t phụ thuộc vào kích thước tập dữ liệu mẫu n.

Kết quả là, đối với phân phối t, công thức tính khoảng tin cậy tương tự như những gì ta đã làm ở phần trước. Tuy nhiên, thay vì lấy phân vị cho phân phối chuẩn, ta sử dụng phân vị cho phân phối t. Cho nên, công thức sẽ là Est +/- t-quantile *std error(Est). Điểm lưu ý là ta sẽ sử dụng độ lệch chuẩn của tập dữ liệu mẫu khi ước lượng độ lỗi chuẩn của Est. Xem thêm

Khoảng tin cậy t rất tiện dụng, nếu bạn được chọn một trong hai với phân phối chuẩn, bạn nên chọn phân phối t. Ta sẽ thấy khi tập dữ liệu càng lớn, t xấp xĩ phân phối chuẩn.

Phân phối t

Phân phối t, được phát minh bởi William Gosset vào năm 1908, có phần đuôi “mập” hơn phân phối chuẩn. Thay vì có hai thông số mean và variance như phân phối chuẩn, phân phối t chỉ có duy nhất thông số bậc tự do (degrees of freedom (df)).

Khi df tăng, phân phối t trông giống phân phối chuẩn, trọng tâm tiến gần về 0. t cũng giải định dữ liệu là iid từ Gaussian, nên thống kê (X’ – mu)/(s/sqrt(n)) có n-1 bậc tự do. Xem thêm

Để minh họa điều trên, ta xây dựng hàm myplot với đối số là df. Hàm này sẽ xuất ra biểu đồ phân phối t với số bậc tự do là df.

k xvals myplot d x = xvals,
dist = factor(rep(c("Normal", "T"), c(k,k))))
g g print(g)
}

Ta thử hàm myplot với bậc tự do là 2

myplot(2)
T plot 2 df

T plot 2 df

Ta thấy phần “mông” của phân phối t (màu xanh) không cao như phân phối chuẩn. Theo đó, hai phần “đuôi” của phân phối t nhận nhiều giá trị hơn nên dày hơn phân phối chuẩn. Chú ý rằng với bậc tự do là 2, ta chỉ có 3 điểm dữ liệu. Tiếp theo, ta thử myplot với df là 20.

myplot(20)
T plot 20 df

T plot 20 df

Cả hai phân phối gần như chồng lên nhau khi ta nâng df lên. Một cách khác để so sánh hai phân phối này là sử dụng phân vị (quantiles). Ta xây dựng hàm myplot2 xuất ra biểu đồ với đường màu xanh nhạt thể hiện phân vị của phân phối chuẩn và đường màu đen thể hiện phân vị của phân phối t. Cả hai phân vị đều bắt đầu ở phần trăm phân vị 50th đến 99th.

pvals myplot2 d p = pvals)
g g g g g print(g)
}

Ta thử dùng hàm myplot2 với 2 df

myplot2(2)
T quantile 2 df

T quantile 2 df

Khoảng cách giữa hai đường này thể hiện sai biệt về kích thước phân vị của chúng. Lưu ý đường kẻ mảnh hơn theo chiều ngang và dọc. Chúng thể hiện 0.975 phân vị của phân phối t và chuẩn. Dù gì thì bạn cũng nhận ra vị trí của đường kẻ dọc là 1.96 như bài viết về tiệm cận (Asymptotics).

Kiểm tra vị trí đường kẻ ngang bằng cách sử dụng hàm qt() của R với đối số phân vị là 0.975 và bậc tự do df là 2.

qt(.975,2)
[1] 4.302653

Ta thấy giá trị này trùng với đường kẻ ngang cắt qua trong biểu đồ trên. Bây giờ, ta chạy hàm myplot2 với 20 df.

myplot2(20)
T quantile 20 df

T quantile 20 df

Hai phân vị ngày càng gần nhau hơn khi bậc tự do càng cao. Tại phần trăm phân vị 97.5, phân vị của t vẫn lớn hơn phân phối chuẩn. Đây là luật Student.

Student’s Rules

Điều này có nghĩa là khoảng tin cậy của t luôn rộng hơn phân bố chuẩn. Do việc ước lượng độ lệch chuẩn phát sinh tính toán không chắc chắn nên kết quả tin cậy cần rộng hơn.

Vậy khoảng tin cậy t được định nghĩa bởi X’ +/- t_(n-1)*s/sqrt(n). Trong đó, t_(n-1) là phân vị. Khoảng tin cậy t giả định rằng tập dữ liệu là idd từ phân phối chuẩn. Mặc dù với giả định cứng nhắc, phép tính này hoạt động hiệu quả với bất kì phân bố dữ liệu nào có tính đối xứng và có dạng hình chuông.

Tuy nhiên khoảng tin cậy t không phải lúc nào cũng khả dụng. Với phân phối lệch, tinh thần giả định (trọng tâm nằm ở 0) của khoảng tin cậy t bị vi phạm. Ta có thể điều chỉnh lại bằng cách lấy logs hay sử dụng giá trị trung vị median.

Skewed Distributions

Skewed Distributions

Đối với dữ liệu rời rạc như phân phối nhị phân, ta cần áp dụng các khoảng tin cậy khác ngoài t. Tuy nhiên, các quan sát cặp (paired observations) thường được phân tích bởi khoảng tin cậy t bằng cách lấy sai biệt giữa các dữ liệu quan sát. Ta sẽ giải thích điều này ngay bây giờ.

Ví dụ về tập dữ liệu Sleep

Hy vọng các bạn không cảm thấy mệt bởi vì chúng ta sẽ làm việc với tập dữ liệu sleep. Đây là tập dữ liệu được phân tích trong bài báo Biometrika của Gosset, cho thấy việc tăng số giờ ngủ của 10 bệnh nhân trên hai loại thuốc ngủ.

sleep
extra group ID
1 0.7 1 1
2 -1.6 1 2
3 -0.2 1 3
4 -1.2 1 4
5 -0.1 1 5
6 3.4 1 6
7 3.7 1 7
8 0.8 1 8
9 0.0 1 9
10 2.0 1 10
11 1.9 2 1
12 0.8 2 2
13 1.1 2 3
14 0.1 2 4
15 -0.1 2 5
16 4.4 2 6
17 5.5 2 7
18 1.6 2 8
19 4.6 2 9
20 3.4 2 10

Sau khi load tập dữ liệu sleep lên, ta thấy R chia tập dữ liệu thành hai nhóm (group) thay vì theo cặp (paired). Ta có tổng cộng 20 dòng dữ liệu, 10 dòng đầu cho kết quả (extra) của nhóm thuốc đầu tiên (group 1) trên từng bệnh nhân (ID), và 10 dòng sau là kết quả nhóm thuốc thứ hai (group 2) trên từng bệnh nhân (ID).

Sleep plot

Sleep plot

Biểu đồ trên biểu diễn dữ liệu theo cặp, các liên kết giữa 2 nhóm thuốc cho từng bệnh nhân được biểu diễn bằng một đường thẳng, nhóm 1 nằm bên tay trái và nhóm 2 nằm bên tay phải. Ta để ý đường màu tím có hệ số góc cao. Đó là bệnh nhân có ID 9, kết quả dịch chuyển từ 0 đến 4.6.

Nếu ta chỉ quan sát trên 20 điểm dữ liệu thì ta sẽ so sánh độ sai biệt (variations) của nhóm 1 và độ sai biệt của nhóm 2. Cả hai nhóm đều có khoảng giá trị lớn. Tuy nhiên, khi chúng ta quan sát theo cặp dữ liệu của từng bệnh nhân, ta thấy sai biệt này rất nhỏ và phụ thuộc vào từng đối tượng cụ thể.

Để làm rõ, ta định nghĩa hai biến g1 và g2. Đây là hai vector có độ dài là 10, ứng với các kết quả khảo sát được của từng bệnh nhân trên 2 loại thuốc. Ta hãy quan sát khoảng giá trị của g1 và g2.

range(g1)
[1] -1.6 3.7

range(g2)
[1] -0.1 5.5

Bây giờ, ta sẽ xem độ sai biệt giữa hai nhóm này.

difference [1] -1.2 -2.4 -1.3 -1.3 0.0 -1.0 -1.8 -0.8 -4.6 -1.4

Tiếp theo ta sử dụng hàm mean để tìm giá trị trung bình của độ sai biệt trên.

mn [1] -1.58

Ta thấy giá trị trung bình này nhỏ hơn giá trị sai biệt trong từng nhóm. Bây giờ, ta sử dụng hàm sd() để tìm độ lệch chuẩn của biến difference.

s [1] 1.229995

Nhắc lại công thức tính khoảng tin cậy T, X’ +/- t_(n-1)*s/sqrt(n). Ta cần tính khoảng tin cậy 95% cho giá trị trung bình của biến difference chúng ta vừa kiếm được. Lưu ý khi lập trình R, ta có thể sử dụng c(-1, 1) để dựng phép toán +/-, hàm qt() để tính phân vị 0.975 và bậc tự do bằng n-1.

mn + c(-1,1)*qt(.975,9)*s/sqrt(10)
[1] 0.7001142 2.4598858

Kết quả trên có nghĩa là với xác suất 0.95, sự khác biệt trung bình về hiệu quả giữa hai nhóm thuốc cho từng bệnh nhân là 0.7 đến 2.46 giờ.

Ta cũng có thể sử dụng hàm t.test với đối số là difference để nhận cùng giá trị trên.

t.test(difference)$conf.int
[1] 0.7001142 2.4598858
attr(,"conf.level")
[1] 0.95

Khoảng tin cậy với hai nhóm độc lập

Bây giờ chúng ta sẽ làm quen với phương pháp, sử dụng khoảng tin cậy t, để so sánh các nhóm độc lập.

Giả sử chúng ta muốn so sánh huyết áp trung bình giữa hai nhóm trong một thử nghiệm ngẫu nhiên. Ta sẽ so sánh những người được điều trị với những người được sử dụng giả dược (placebo). Không như nghiên cứu về giấc ngủ, ta không thể dùng cặp (paired) cho kiểm định t bởi vì các nhóm độc lập và có kích thước mẫu khác nhau.

Như vậy, mục tiêu của chúng ta là tìm khoảng tin cậy 95% cho độ sai biệt giữa hai trung bình quần thể này. Ta biểu diễn độ sai biệt này bằng mu_y – mu_x.

Nhắc lại công thức của chúng ta X’ +/- t_(n-1)*s/sqrt(n). Đầu tiên chúng ta cần trung bình mẫu, nhưng chúng ta có tới hai X’ và Y’, tương ứng với trung bình mẫu của từng nhóm. Cũng dễ hiểu khi chúng ta có thể lấy (Y’-X’) vì chúng ta đang tìm khoảng tin cậy có chứa độ sai biệt mu_y – mu_x. Bây giờ, ta chỉ định phân vị t. Giả sử hai nhóm có hai kích thước mẫu n_x và n_y.

Đối với một nhóm ta sử dụng phân vị t_(0.975,n-1). Trong trường hợp hai nhóm, ta có thể sử dụng phân vị t_(.975,n_x+n_y-2). Ta trừ 1 bậc tự do ở mỗi nhóm vì ta đang tính trung bình mẫu cho từng nhóm nên ta cộng n_x và n_y sau đó trừ đi 2.

Toán hạng còn lại là độ lỗi chuẩn đối với một nhóm được tính bằng s/sqrt(n). Ta giải quyết phần tử số trước. Giả sử phương sai của hai nhóm là s^2, ta lấy được giá trị này bằng cách sử dụng tổng có trọng số (weighted sum).

Ta gọi phương sai ước lượng là phương sai gộp (pooled variance). Công thức cho phương sai này cần hai phương sai ước lượng, S_x và S_y tương ứng cho từng nhóm. Ta nhân cho bậc tự do cho từng S_x và S_y, sau đó chia cho tổng của các bậc tự do này. Do vậy, mẫu dữ liệu càng lớn sẽ có trọng số càng cao.

Công thức ở tử số là: (n_x-1)(S_x)^2+(n_y-1)(S_y)^2
Công thức tổng bậc tự do: (n_x-1)+(n_y-1)
Công thức ở mẫu số: sqrt(1/n_x + 1/n_y)

Nhóm đầu tiên là 8 người dùng thuốc ngừa thai (oral contraceptive) và nhóm hai gồm 21 người (controls). Trung bình mẫu tương ứng là X’_{oc}=132.86 và X’_{c}=127.44. Độ lệch chuẩn từng mẫu là s_{oc}= 15.34 và s_{c}= 18.23.

Phương sai mẫu được tính gộp như sau

sp [1] 8293.867

Tổng bậc tự do được tính gộp như sau

ns [1] 27

Độ lỗi chuẩn được tính gộp như sau

sp [1] 17.52656

Như vậy, khoảng tin cậy 95% được tính như sau

132.86-127.44+c(-1,1)*qt(.975,ns)*sp*sqrt(1/8+1/21)
[1] -9.521097 20.361097

Ta để ý giá trị 0 nằm trong khoảng tin cậy 95%. Điều này có nghĩa là ta không thể loại trừ khả các trung bình mẫu của hai nhóm bằng nhau.

Các bạn đã mệt chưa? Hãy quay lại ví dụ về thuốc ngủ và thay vì quan sát dưới gốc độ cặp dữ liệu cho 10 đối tượng nghiên cứu ta quan sát dưới gốc độ là hai tập độc lập có kích thước 10. Nhắc lại ta đã lưu dữ liệu dưới dạng hai vector g1 và g2, giá trị sai biệt trung bình giữa hai nhóm là md.

Ta hãy tính phương sai gộp của hai nhóm và lưu vào biến sp. Nhắc lại đó là giá trị sqrt(tổng có trọng số của phương sai mẫu/bậc tự do). Trọng số tương ứng là kích thước mẫu – 1. Bậc tự do lúc này là 10+10-2 = 18.

sp [1] 1.898625

Toán hạng sau cùng trong công thức là độ lỗi chuẩn, đơn giản bằng cách lấy sp nhân cho căn bậc hai của 1/10 + 1/10. Ta tính khoảng tin cậy t 95% cho độ sai biệt trung bình mẫu giữa hai nhóm g1 và g2 như sau

md + c(-1,1)*qt(.975,18)*sp*sqrt(1/5)
[1] -0.203874 3.363874

Ta có thể kiểm chứng lại công thức bằng tay trên bằng hàm t.test của R. Do ta lấy g2 – g1, nên đối số đầu tiên phải là g2 sau đó là g1.

t.test(g2,g1,paired=FALSE,var.equal=TRUE)$conf
[1] -0.203874 3.363874
attr(,"conf.level")
[1] 0.95

Khá hay vì nó trùng khớp phải không nào? Để ý rằng giá trị 0 một lần nữa nằm trong khoảng độ tin cậy 95%, do đó ta không thể loại trừ phát biểu rằng hai nhóm này như nhau. Tuy nhiên điều này không như kết quả ta thấy trước đó khi khảo sát theo gốc độ cặp dữ liệu. Ta dùng t.test chạy lại kết quả trước đó với đối số pair=TRUE.

t.test(g2,g1,paired=TRUE)$conf
[1] 0.7001142 2.4598858
attr(,"conf.level")
[1] 0.95

Ta thấy giá trị 0 không nằm trong khoảng tin cậy này. Điều này có nghĩa rằng hai nhóm khi tạo cặp sẽ có độ sai biệt trung bình lớn.

Khoảng tin cậy với hai nhóm có phương sai khác nhau

Bây giờ, ta sẽ nói về cách tính khoảng tin cậy cho hai nhóm có phương sai khác nhau. Ta không gộp chúng như những gì ta đã làm trước đó.

Trong trường hợp này, công thức tính khoảng tin cậy tương tự như chúng ta thấy trước đó, Y’-X’ +/- t_df * SE. Tuy nhiên, độ lỗi chuẩn SE và phân vị t_df được tính theo cách khác. SE là căn bậc hai tổng của căn bậc hai độ lỗi chuẩn ứng với hai trung bình mẫu, (s_1)^2/n_1 + (s_2)^2/n_2.

Khi dữ liệu của X và Y là iid từ phân phối chuẩn và phương sai của chúng khác nhau, thống kê chuẩn ta học từ trước, (X’-mu)/(s/sqrt(n)), không tuân theo phân phối t. Tuy nhiên, nó có thể được xấp xĩ bởi phân phối t nếu ta xét bậc tự do một cách hợp lý.

Công thức tính bậc tự do trong trường hợp này khá phức tạp đến nỗi không ai nhớ. Trên tử số là bình phương (tổng của (bình phương độ lỗi chuẩn của hai nhóm – s^2/n))). Ở mẫu số là tổng của hai thành phần tương ứng với mỗi nhóm. Mỗi thành phần có dạng như nhau, là độ lỗi chuẩn của trung bình mẫu chia cho kích thước mẫu – 1, (s^4/n^2)/(n-1). Ta có thể dùng bậc tự do df này để tìm phân vị t.

t\verb|-|interval = \bar{Y} - \bar{X} \pm t_{df} \times (\frac{s_x^2}{n_x} + \frac{s_y^2}{n_y})^{1/2}

df = \frac{(s_x^2/n_x + s_y^2/n_y)^2}{(\frac{s_x^2}{n_x})^2/(n_x - 1) + (\frac{s_y^2}{n_y})^2/(n_y - 1)}

Ta thử gán các giá trị cho công thức trên trong ví dụ về huyết áp. Nhắc lại ta có hai nhóm, nhóm một có số lượng 8 và X’_{oc}=132.86, s_{oc}=15.34, nhóm hai có số lượng là 21 và X’_{c}=127.44, s_{c}=18.23.

Ta tính bậc tự do trước. Bắt đầu với tử số (numerator), là bình phương tổng của hai thành phần. Mỗi thành phần có dạng s^2/n.

num [1] 2046.642

Bây giờ đến mẫu số (denominator), là tổng của hai thành phần. Mỗi thành phần có dạng s^4/n^2/(n-1).

den [1] 136.1235

Bây giờ chia tử với mẫu ta sẽ có mydf.

mydf [1] 15.03518

Bây giờ sử dụng hàm qt(.975,mydf) để tính khoảng tinh cậy t 95%. Nhắc lại công thức, X’_{oc}-X’_{c} +/- t_df * SE. Trong đó, SE là (s_1)^2/n_1 + (s_2)^2/n_2.

132.86-127.44 +c(-1,1)*qt(.975,mydf)*sqrt(15.34^2/8 + 18.23^2/21)
[1] -8.913327 19.753327

Tương tự, ta có thể dùng t.test với var.equal=FALSE, khi đó R sẽ tính bậc tự do cho bạn. Bạn không cần phải nhớ công thức.

Xin chúc mừng, bạn đã hoàn tất bài viết về t-dious (buồn chán) này.

Nguồn tham khảo: http://swirlstats.com/

Advertisements

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Đăng xuất / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Đăng xuất / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Đăng xuất / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Đăng xuất / Thay đổi )

Connecting to %s